texts: How to access a Feature service hosted on arcgis.com which is shared only with particular groups or organization

texts:
In my C# application, I am trying to access a Feature service which is hosted on arcgis.com. The feature service is shared only with a group or organization and is NOT public. How can I access this service in my client code. Here is my client code which throws an exception saying the Authentication Token required. I do have the token but dont know how to use that token to be authenticated to use the service. This same code runs fine for services that are shared as PUBLIC.

   IPropertySet sipPS = new PropertySet();
    sipPS.SetProperty("DATABASE", layerUrl);
    IWorkspaceFactory sipWSF = (IWorkspaceFactory)new ESRI.ArcGIS.Carto.FeatureServiceWorkspaceFactory();
    IWorkspace sipWS = sipWSF.Open(sipPS, 0);
    IFeatureWorkspace sipFWS = (IFeatureWorkspace)sipWS;
    IGroupLayer grpLayer = new GroupLayerClass();
    grpLayer.Name = Constants.ArcGISOnlineGroupLayerName;
    IFeatureClass sipFC = sipFWS.OpenFeatureClass("0");
    IFeatureLayer fl = new FeatureLayerClass();
    fl.Name = featureLayerName;
    fl.FeatureClass = sipFC;
    grpLayer.Add(fl as ILayer);


texts: How do I rename an out db POSTGIS raster

texts:
I used the following command to import a bunch of raster to postgis out-db style (only pointers to the rasters are stored in the database):

raster2pgsql -s {srid} -c -R -I -C -F -t auto {file_name} public.{table} | psql -h {host} -p {port} -d {database} -U {user} --no-password

Now I want to rename all those files that I previously imported. I already renamed them on the filesystem. Now when I try to retrieve them I get:

psycopg2.InternalError: rt_band_load_offline_data: Cannot open offline raster: /myfilepath/file_name.tif

Where in the Postgis db to I rename the files (edit the pointer path) so that I can access them again? I'm hoping I can do this without re-importing everything again using raster2psgsql.


texts:
I need a open-source server-side component for a web-application that can export basemaps + layers to an image format (saved on a db).
Being slightly new to GIS, I've managed to find a few tools that do this on the client-side.

leaflet's plugin Leaflet.print does this. I'm using leaflet to display CloudMade tiles on the browser.  it has npm packages, so I'm guessing it can run server-side as well.
OpenLayers 3 (beta) has this as an example but I've not tested it out.
MapFish Print
GeoServer which uses MapFish Print again.

Is there any component that could take in map bounds & a features as geoJson and output image file on the server-side ?
It would not be possible to use ArcGIS or proprietary tools


texts: How to keep a parameter blank in a script tool?

texts:
I've developed a sort of find and replace script tool that functions. But I'm not able to figure out how to find empty cells in fields and replace them.For the "Find" parameter, I have tried entering "", '', None, null, and < null >.

An example that might make this clearer. Say you had a field in a table with every row empty, and you wanted to make them say "Not Available" or something instead, which you would enter in the "replace" parameter.

I feel like this is a tool/syntax issue and not really a script one, no?

import arcpy

FC = arcpy.GetParameterAsText(0)
oldText = arcpy.GetParameterAsText(1)
replaceText = arcpy.GetParameterAsText(2)
targetField = arcpy.GetParameterAsText(3)

queryString = targetField + ' = ' + "'" + oldText + "'"
try:

    with arcpy.da.UpdateCursor(FC, (targetField), queryString) as cursor:
        for row in cursor:
            if row[0] == oldText:
                row[0] = replaceText
                cursor.updateRow(row)

                arcpy.AddMessage("String where " + queryString + " found and      replaced with " + replaceText)

    except:
        arcpy.AddMessage(arcpy.GetMessages())
        arcpy.AddError("Failed to find and replace where " + queryString)


texts:
I am building a GIS system. All my data is stuffed into GeoServer. I need to have a front end able to do the following:

Select one or several layers from GeoServer and show it to the end user on a map. It is important to say that there are many layers (speaking 1000s or so) so I want to make a search and pile up some layers the user has chosen and then generate a map with those layers, being able the check boxes for display
Create some group of layers to simplify the work as in point 1.
To be nice, i.e. a bit like extjs (but extjs is out for licensing issues) having the possibilities to overlay windows on top of the map and collapse them if needed.
Preferably based on Java, JavaScript, but not mandatory

So far I have checked out different things but it is always one piece of the stack. OpenLayers, extjs, mapfish....

Anything out there?

If there is nothing, how many people would be ready to team up to come up with open source stuff. I was thinking about a grails plugin.


texts: Polygon Neighbors tool produces incorrect results

texts:
BACKGROUND:
At 10.8.1, I have used the Polygon Neighbors tool with a countywide parcel shapefile (approximately 47,000 polygons).  The tool runs to completion (taking < 1 minute) without reporting any errors.  I used the shapefile's PARCEL_NO field as the "Report by Field", and I unchecked "Include both sides..." option.  All other settings were left at their default.
THE PROBLEM:
When I manually check the results, I find that Polygon Neighbors has failed to identify many of the neighboring polygons.  For example, one of the ~47k polygons has 15 neighbors, but the tool only found 4 neighbors.
In several cases, none of the surrounding polygon neighbors were identified.
Out of approximately 20 manual checks, only one parcel's neighbors were correctly identified, and that was a simple situation with only two neighbors.
ATTEMPTED SOLUTIONS:

Ran the Geometry Check tool on the parcel shapefile.  It reports ok.
Ran Polygon Neighbors from ArcCatalog instead of ArcMap.  The same errors resulted.
Selected a subset (169) of the ~47k polygons and Data > Export them to a new shapefile.  This subset included the polygon with the 15 neighbors mentioned above.  In this case, Polygon Neighbors correctly(!) identified all neighbors.
Used QGIS to perform the same polygon neighbor concept with the ~47k polygon shapefile, using the methods described in https://spatialthoughts.com/2019/05/23/neighbor-polygons-aggregate-qgis/.  This method correctly identified all neighbors.  However, it took about 20 hours to complete, so I'd like to find a faster option.

QUESTION:
I've run out of ideas.  Why does Polygon Neighbors fail to identify surrounding polygons in the ~47k polygon shapefile, yet correctly identify them in a much smaller shapefile subset?


texts: How to find 2 parallel lines with given line?

texts:
I have a line that is drawn on a map.
I want to find 2 parallel lines with this one, one at 50m above this and one at 50m below this.
Something like in the picture below

In the end I want to transform the line in a rectangle of width 100m.

I have tried some implementations, but I can't get it right.

How can I achieve this, keeping in mind that the line can have a length of hundreds of kilometers?

Edit:

  function getBearing(a, b) {
    const TWOPI = 6.2831853071795865;
    const RAD2DEG = 57.2957795130823209;
    var theta = Math.atan2(b.lat - a.lat, a.lng - b.lng);
    if (theta < 0.0)
      theta += TWOPI;
    console.log(a, b);
    return RAD2DEG * theta;
  }

  function getPolygonForLine(start, end, distance) {
    var bearing = getBearing(start, end);
    distance /= 2;
    var dx = distance * Math.cos(bearing);
    var dy = distance * Math.sin(bearing);

    var delta = {
      lat: dx / (111320 * Math.cos(start.lat)),
      lng: dy / 110540
    };

    var bounds = [
      L.polygon(getBoundsForPoint(start, delta)).getBounds(),
    ];

    var points = [];
    _.forEach(bounds, function(val) {
      points.push(
        val.getNorthEast(),
        val.getNorthWest(),
        val.getSouthWest(),
        val.getSouthEast()
      );
    });
  }

  function getBoundsForPoint(point, delta) {
    var points = [
      [
        point.lat + delta.lat,
        point.lng + delta.lng
      ],
      [
        point.lat - delta.lat,
        point.lng - delta.lng
      ],
      [
        point.lat + delta.lat,
        point.lng - delta.lng
      ],
      [
        point.lat - delta.lat,
        point.lng + delta.lng
      ],
      [
        point.lat - delta.lat,
        point.lng
      ],
      [
        point.lat + delta.lat,
        point.lng
      ],
      [
        point.lat,
        point.lng - delta.lng
      ],
      [
        point.lat,
        point.lng + delta.lng
      ]
    ];
    return points;
  }

This is the code that I have at the moment(I left only the code that calculates one end of the line). The problem is that the rectangle seems to not have the desired width and it also does not form a 90 degree angle with the line (it is parallel with the edge of the screen).


texts: A list of world projection that make "beautiful" map in QGIS

texts:
I was making a globe world map with the "on the fly" function in the project Properties, which changing all the project CRS simultaneous
The "World_Robinson" projection make a beautiful oval map. I was wondering what other projection can be used for a beautiful 3D map?


texts: How do I password protect a CartoDB Visualization

texts:
We have some data and a visualization that we would like to make available to a small user community.

How do I password protect a CartoDB visualization without making it public.


texts: cartodb: how to display data from a related table in the info window?

texts: I have a table of buildings. I have a table of all the things in each building. It is a one to many relationship. How can i make info windows which display a list of the things in the building?


texts: What does pgRouting's architecture look like?

texts:
I know that pgRouting is a pretty nice tool, but I would like to use something of my own in my application, mainly because I don't know C and I don't know how to tweak the performance of pgRouting. Thus I am interested in how pgRouting's internals are organized.

How is the actual shortest path calculation done? Is it anything like the following:

Load the whole ways/edges table into memory with one 'select * from ways';
Build a graph from that;
Apply a C-based (based on the Boost library) shortest path algorithm to the graph;
Return a result into psql server?

The other interesting thing is - why does pgRouting store vertices coordinates (x1, y1, x2, y2) in the ways table? Is it the only table from which the graph is built?


texts: How to creat a feature type using rest/ php on geoserver?

texts:
I'm trying to create a feature type for a ship in geoserver but when I post the request to the server I get this response:

[geoserver.rest] - Trying to create new feature type inside the store, but no attributes were specified

This is xml data that I'm sending to the server :

<featureType>
        <name>'.$layerName.'</name>
        <enabled>true</enabled>
        <srs>RGF93_Lambert_93</srs>
        <projectionPolicy>FORCE_DECLARED</projectionPolicy>
        <nativeName>'.$layerName.'</nativeName>
        <description>'.htmlentities($description, ENT_COMPAT).'</description>
        <store class="dataStore"><name>'.htmlentities($datastoreName, ENT_COMPAT).'</name></store>
</featureType>'


texts: Randomly selecting single point from within multiple fishnet grid-squares?

texts:
Using ArcGIS, I have a large data set of occurrence data many of which are in similar places. I have created a fishnet around the total area and want to select a single random point from within each grid-square of the fishnet.

How do I go about this?

I want to run this multiple times so that I can then bootstrap it. I know I could subset each individual grid-square but there are approximately 100 of them so this would take a very long time.


texts: QGIS: CRS, ellipsoids and area calculation

texts:
I have a doubt about dynamic CRS, ellipsoids and area calculations. I am using QGIS 2.14.0 (still waiting for IT to install 2.14.11...)

When I load a polygon table into a new, empty project, it comes in as EPSG:32723. Looking at the Project Properties > General tab, no ellipsoid is defined. When I go to the CRS tab and activate dynamic transformation, I see that the ellipsoid has now been defined in the General tab (corresponding to the ellipsoid used by EPSG:32723).

Does this mean that until I tell QGIS to activate dynamic transformation, it is ignoring ellipsoid parameters for area and distance calculations? Because when I calculate area values for these polygons using EPSG:32723 and EPSG:32723 (Dynamic), I get different results.

EDIT: I saw a question similar to this from 2015, but it (and the replies) make no mention of the ellipsoid parameters being used or not used for area calculations.


texts: Get pixel count according to a shapefile's identity

texts:
I have a raster that has 5 classes and a polygon shapefile with 10 rows.

This is the raster:

This is the shapefile with various shapes

I want to get the pixel count

by classes - 5
by identity (there are 10 in the shapefile).


texts:
I am trying to upload a geoTIFF file into my PostGIS database. Using raster2psql, it works.

raster2pgsql -I -C -M *.tif -F -t auto mytablename > mytable.sql
psql -d mydatabase -f mytable.sql

However, my datatable is not what I expect it to be.

I get:

rid | rast | filename
----+------+----------
 1    01000001000F1111111111813F0F111111111181BF00000000008066C0FAFFFFFFFF7F564000000000000000000000000000000000E61000001E001E006AFDFF7FFFFDFF7... filename.tif

What I want is:

rid | geometry | value | filename
 1    0100AF...   12.3   filename.tif

Where geometry is the hex code associated, and the value is the band value from the raster.

How would I go about doing that?


texts:
How to use QGIS batch-processing dissolve:

Using QGIS 2.16 under Win 10 (see screenshot). Questions:

Why must I have at least two rows. When I open the window I have 3 rows. Using the red minus-icon, I can delete one row, but not two rows. I think I only need one row, because I want to dissolve only one file (layer).
The last right-most cell of each row is for the output of the dissolve-process, right? Should I simple enter a filename like output.shp here?
Anybody got this working under Windows?


texts: Iterate into Calculate field to populate Field with unique condition rating

texts:
I am new to ModelBuilder and I've been learning as I go. Almost have this model working but I've hit a snag. I have a simple Calculation that I iterate to spit out condition ratings of pipelines based on date and durability etc. This I have working fine. Problem is writing to the new field I created. I am trying to use Calculate Field inline with the iterated calculation so as it calculates it writes to my field. Problem is it updates the entire field every time it iterates over the calculation.

How can I make calculate field only write one value at a time and skip to the next cell to populate my field with satisfying pipeline condition ratings?


texts: Why am I getting small float numbers when aggregating points containing large integer values in QGIS?

texts:
I have a layer of points containing the centroids of European NUTS regions at NUTS 3, with population data in the atttribute table. I want to aggregate this population data to NUTS level 1, which is contained in a polygon layer.

As such, I use the QGIS tool Merge attributes by location. The operation finishes without error, but the resulting summed population values in the created polygons are all very small float numbers. See the picture below for examples on input and output data:

I have double checked that the source data is integers, and are formatted correctly. I also tried running the tool both from the vector menu as well as the processing toolbox.

Am I doing something wrong, or is this a bug in the QGIS spatial join?


texts: How to display a control (zoom bar) outside map container in openlayers?

texts:
Since I'm new to OpenLayers this question might be the simplest of them all but after a quick search on the subject I couldn't find anything useful, so:

Is it possible to put a OpenLayers zoom bar/slider outside the map container in another div? Like for example OpenLayers.Control.Zoom().


texts: R: distance from one point to many polygons

texts:
OK so I have a spatial polygons data frame, and I have a point of interest. I want to add a row to the SPDF that contains the distance from the centroid of each polygon to this given point. How do I do this ?

Here is a toy dataset:

    library(sp)
    library(spdep)
    library(rgeos)

    # Create toy SPDF
    fn <- system.file("etc/shapes/eire.shp", package = "spdep")[1]
    prj <- CRS("+proj=utm +zone=30 +units=km")
    eire <- readShapeSpatial(fn, ID = "names", proj4string = prj)

    # Create a point of interest
    pt.coords <-readWKT("POINT(260 5900)")
    eire.proj<-CRS(proj4string(eire))
    proj4string(pt.coords) <- eire.proj

    # check point is on the map
    plot(eire)
    plot(pt.coords, col="Red", add=TRUE)

Ok so the SPDF data has following columns:

    > names(eire)
   [1] "A"       "towns"   "pale"    "size"    "ROADACC" "OWNCONS" "POPCHG"  "RETSALE" "INCOME"  "names"

I want to add a extra column showing the distance from each polygon centroid to the point we defined. How do I do this ?


texts: Wildcard for field in arcpy.GetParameterAsText

texts:
I have this parameters for my python script:

import arcpy

Layer1 = arcpy.GetParameterAsText(0)
Tabela1 = arcpy.GetParameterAsText(1)
Estado = arcpy.GetParameterAsText(2)
Ano = arcpy.GetParameterAsText(3).split(';')
mxd = arcpy.mapping.MapDocument("CURRENT")
Saida = arcpy.GetParameterAsText(4)
stats = []

Parameter 'Ano' is a multivalue field selection and this parameter obtained the value list from the feature class 'Tabela1'. There is many fields in this FC and all of them are listed in 'Ano' but I only want to show the fields with a prefix ACH*. The script works well, But don`t want to get the user confusing with too many fields in the 'Ano' field.


texts: Select raster tiles within vector polygons

texts:
I have a shapefile with polygons of all the countries of the world. I also have a raster .tiff (250m square pixels) of average wind speed per pixel.

I want to select all the tiles which lie within each country polygon so I end up with lots of individual raster files which correspond to every single country.

I know there will be border/overlap issues but it's supposed to be general. How can I do this in QGIS?


texts: Rasterizing point layer with multiple features falling inside raster cell using GDAL?

texts:
How to make a raster with values that are the sum of all points inside each raster cell?

The burn value is a specific attribute.

Using gdal_rasterize or gdal.RasterizeLayer(), the resulting raster only has the value of a random (or probably the first or last) point inside a raster cell. I’ve already tried gdal_rasterize option ALL_TOUCHED, but it does not change the behaviour.


texts: Are there any examples of the Esri GeoPortal REST API in use?

texts:
Are they any public examples of sites built on top of the Esri GeoPortal's REST API?

I would like to use the back-end functionality of the GeoPortal to manage the metadata, but serve it up using a custom front-end.

Before I undertake this task I'm interested to see whether anyone has successfully achieved this.

Thanks.


texts: Setting final orientation of polygons in CityGML model

texts: I am converting shapefile model of buildings in LoD2 to CityGML. In CityGML really matter about poylgon orientation(Countercloskwise can be seen from outside, Clockwise from inside). I have problem with orientation of some polygons.  Althought I've set almost every polygon to right direction, there still remain some of them to orientate. I know that in 2D it's pretty much easy task, but it changes by going to 3D (polygonZ). Is there some kind of tool, which can solve this task straight forward?


texts: Are ESRI Shapefiles and KML files two different ways of representing data?

texts:
I have a really basic newbie question. I have following data with me:

Esri Shapefiles (.shp, .shx, .dbf, .prj, .sbn and .sbx)
KML/KMZ files

Are Esri Shapefiles and KML files two different ways of representing data? I am suspecting that the data provider has given me the same data in two formats: Shapefiles and KML files.

Is there a way I can verify that both of them are representing same data and that there is no loss of information when using one over the other?


texts: Moving polyline to another location in ArcGIS Desktop?

texts:
In desktop ArcGIS 10.3, I have this polyline with xy center coordinate of it is:

18157037.6 / 65108911.9

When I zoom out (1:5000000) the map, in order to grab the polyline into the right location,  I can't see the right located layers and the problematic polyline.

I understand the polyline has a different coordinate system- but when i see the the properties of the layer- it has no coordinate system:

and I try to move this polyline (center of it) into 181586/651079 (the right coordinate system).

I tried to use the tool "Move" in the editor- but the polyline just located somewhere in the space - but not in the wanted xy coordinate.

Also, I tried to convert the coordinate system of the polyline into Israel TM grid- but with no success.


texts: Difficulty using variables as input parameters for arcpy.Intersect_analysis

texts:
I have 2 variables created from a previous geoprocesses:

fc = arcpy.JSONToFeatures_conversion("cslf.json", os.path.join("cslf.gdb", "cslf"))
selectFeatures = arcpy.CopyFeatures_management(select, "C:\Workspace\Sandbox\MapChangeProject\selects.gdb\select_{}".format(item))

Both of these layers have viable feature layers that can be added in arcMap. I want to use these 2 variables as the input feature parameters (in_features) for a new feature created with arcpy.Intersect_analysis() like this:

 intersect = arcpy.Intersect_analysis([fc, selectFeatures],"intersects", "ALL", "", "")

However, when I run this, I get the following error: ERROR 000628: Cannot set input into parameter cluster_tolerance. Any suggestions as to how to fix my input parameters? I've tried a few variations (adding ""'s) but I get the same thing.


texts: Adding multiple OpenLayers Vector Layer markers with data parsed from JSON file?

texts:
I would like to parse a local JSON file and add markers to a map based on the latitude and longitude parsed from the JSON file.

Here is a sample from the local JSON file:

"LOCATION_A": {
              "latitude": 10.2070,
              "longitude": 25.3215,
              "altitude": 1.2,
              "status": "CLOSED",
              "serverAddress": "192.168.0.1"
 }

I am unsure how to proceed beyond the following block of code:

var markers = new OpenLayers.Layer.Vector("Sites", {
  projection: map.displayProjection,
  strategies: [new OpenLayers.Strategy.Fixed()],
  protocol: new OpenLayers.Protocol.HTTP({
    url: "filename.json",
    format: new OpenLayers.Format.JSON()
    }),
  });

I am unsure how to use the parser to do this.  I would like to use jQuery for this exercise now.


texts:
I've attempted to write a script that will export all the coded value domains in a geodatabase to a table. Before I dig in the to script and what's wrong, I have looked at the similar question on this topic, but that script would not run for me. Here is my code:

import arcpy

arcpy.env.workspace = r"F:\ITree"
gdb = r"ITree.gdb"
arcpy.env.overwriteOutput = True

domains = arcpy.da.ListDomains(gdb)

for domain in domains:
     if domain.domainType == 'CodedValue':
         domain_name = domain.name
         print 'Exporting {0} coded values to table in {1}'.format(domain_name, gdb)
         coded_value_list = domain.codedValues
         print "The coded values / descriptions are"
         for value, descrip in coded_value_list.iteritems():
             print "{0} : {1}".format(value, descrip)
         out_table_name = domain_name.lower()
         #arcpy.DomainToTable_management(gdb, domain_name, out_table_name, "item", "descrip")
     else:
         print "{0} not a coded value domain. Passing it up.".format(domain.name)

When I comment out the Domain to Table, the script runs and produces, this:

Exporting LANDUSE coded values to table in ITree.gdb
The coded values / descriptions are
1 : Single family residential
0 : Not entered
3 : Small commercial
2 : Multi-family residential
5 : Park/vacant/other
4 : Industrial/Large commercial

However, when I uncomment the Domain to Table, it produces the following error:

Exporting LANDUSE coded values to table in ITree.gdb
The coded values / descriptions are
1 : Single family residential
0 : Not entered
3 : Small commercial
2 : Multi-family residential
5 : Park/vacant/other
4 : Industrial/Large commercial
Traceback (most recent call last):
  File "E:/Work/MDC/GISData/MOTools/Scripts/export_gdb_domains.py", line 29, in <module>
    arcpy.DomainToTable_management(gdb, domain_name, out_table_name, "item", "descrip")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 1491, in DomainToTable
    raise e
arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Failed to execute (DomainToTable).

It also produces the first table from the domain, but does not copy the coded values and creates an empty table. I've looked through the help documents for Domain to Table (Data Management) and List Domains, but I cannot figure out why it's throwing an error. Clearly something is off in the Domain to Table function, but I can't figure out what, especially since it prints all the parameters correctly.

Using ArcGIS 10.3 and PyCharm Community 4.0.5 with Python 2.7.9


texts:
I recently used a GDAL-command working with the OSGeo4W command line, but it does not feel that comfortable. Is there a way to somehow execute a GDAL-command directly in QGIS?

Right now the gdal-extension is installed in QGIS, but i can not find a command line.


texts: Compare multiple lists of field names thought to be in shapefile with multiple list of fields found in shapefile

texts:
In the code below, I am able to create and print multiple lists that contain field names for each shapefile. This is done by creating a list of file paths, pointing to a shapefile, and looping through the list passing each shapefile into the arcpy.ListFields function.

import arcpy

#Define the Workspace
arcpy.env.workspace = r"C:\TestFolder"

dbf1 = "APExistOnlyEdit.dbf"
myfield = "FI_PATH"

filepathlist =[row[0] for row in arcpy.da.SearchCursor(dbf1,myfield)]
#print(filepathlist)

possiblefieldlist =['AP_FIPS','AP_BLDGCOM','AP_TYPE','AP_STATUS','AP_ID','AP_HOUSENU','AP_HALFADD','AP_PREDIR','AP_PRETYPE','AP_STNAME',
                    'AP_SUFTYPE','AP_SUFDIR','AP_UNITTYP','AP_UNIT','AP_BUILDIN','AP_FULLADD','AP_CITY','AP_STATE','AP_ZIP']

#Loop through shapefiles in fielpathlist
#Get field names for each shapefile and pass into lists
for thisFile in filepathlist:
    ContainedFieldNames = [f.name.upper() for f in arcpy.ListFields(thisFile)]
    print(ContainedFieldNames)

#Fields with Field Names thought to be in Shapefile
#Search though Fields and rows and make lists of field names though to be in shapefile

with arcpy.da.SearchCursor(dbf1,possiblefieldlist) as cursor:
    for row in cursor:
        fields = [a.upper() for a in row if a.strip() != '']
        print(fields)

Below is the output when I print the lists of field names found in each shapefile. There are three separate lists, one for each shapefile.

['FID', 'SHAPE', 'FULL_ADDRE', 'STREET_ADD', 'CITY_ST_ZI', 'PREFIX_DIR', 'HALF_ADDRE', 'HOUSE_NUMB', 'STREET_NAM', 'STREET_TYP', 'SUFFIX_DIR', 'UNIT_TYPE', 'UNIT_ID', 'CITY', 'STATE', 'ZIPCODE', 'PARCEL_NUM', 'PARCEL_N_1', 'HOUSE_NU_1', 'HOUSE_NU_2', 'STATUS', 'ASSEMBLY_S', 'SENATE', 'DISTRICT', 'COUNCIL', 'GRIDNUM', 'PHYSICAL_A', 'PUBDATE']

['FID', 'SHAPE', 'ADD_ID', 'HSE_NUMB', 'SNAM_PREMO', 'SNAM_PREDI', 'SNAME', 'SNAM_POSTY', 'SNAM_POSDI', 'SNAM_POSMO', 'FULL_STREE', 'SUB_ADD_TY', 'SUB_ADD_ID', 'FULL_ADDRE', 'PLACE_NAME', 'FNSB_COMM', 'FECC_COMM', 'ZIPCODE', 'ADD_TYPE', 'NOTES', 'STATUS', 'MS_EXCEPTI', 'NO_MSAG', 'GLOBALID', 'CREATED_US', 'CREATED_DA', 'LAST_EDITE', 'LAST_EDI_1']

['FID', 'SHAPE', 'OBJECTID', 'P_ID', 'ACCOUNT', 'TAXID_LOKI', 'ADRSNUM', 'P_ROADNME', 'ROADNME', 'S_ROADNME', 'PS_ROADNME', 'ADRSNUM_S', 'ZIP', 'LAT', 'LONG', 'ADDRESS', 'COMMUNITY', 'GLOBALID']

In addition, I use the arcpy.da.SearchCursor to loop through multiple fields, in a dbf file that contain names of fields, that are thought to be within each shapefile. I stripped out rows that do not contain any text. Below is the output when I print the lists of field names thought to be within the shapefile.

['HOUSE_NUMB', 'HALF_ADDRE', 'PREFIX_DIR', 'STREET_NAM', 'STREET_TYP', 'SUFFIX_DIR', 'UNIT_TYPE', 'UNIT_ID', ' BUILDINGA', 'STREET_ADD', 'CITY', 'STATE', 'ZIPCODE']

['HSE_NUMB', ' PRETYPEA', 'FULL_STREE', 'SUB_ADD_TY', 'SUB_ADD_ID', 'ZIPCODE']

[' TYPEA', 'OBJECTID', 'ADRSNUM', 'P_ROADNME', 'ROADNME', 'S_ROADNME', 'ADRSNUM_S', 'ADDRESS', 'COMMUNITY', 'ZIP']

My goal with this script is to compare the lists of field names though to be within the shapefile (possible fields) to the lists of field names actually found within the shapefile (actual fields). For each comparison of lists, I want to capture the text of field names that are found and field names that are not found. I would like to use arcpy.da.UpdateCursor to fill two new fields, one field will hold the field names that are found and one will hold the field names that are not found. I am using a dbf table to control processing and it holds all field names thought to be within the shapefile and file paths to shapefiles and I want to write results to this dbf table as well. All of the examples that I have seen compare a single shapefile that is defined against one defined list. I am wondering if it's possible to compare multiple lists with multiple lists using ArcPy functions.

Below is a pic of my dbf table and I'm using ArcGIS Pro.


texts: Json not loading from GeoServer into Leaflet?

texts:
I'm creating an web map application using Geoserver and Leaflet. I'm trying to load a few Json feature through an WFS service.The problem is that nothing is being displayed on the final map, the Geoserver console emits the following meassage: [info.json] - about to encode JSON

What makes me to believe that there is nothing wrong with the service.

Here is the code.

var map = L.map('map').setView([39.752171, -104.998817], 17);

L.tileLayer('https://{s}.tiles.mapbox.com/v3/{id}/{z}/{x}/{y}.png', {
        maxZoom: 18,
        attribution: 'Map data &copy; OpenStreetMap contributors, ' +
            'CC-BY-SA, ' +
            'Imagery © Mapbox',
        id: 'examples.map-20v6611k'
    }).addTo(map);

var geojsonLayer = new L.GeoJSON();

function handleJson(data) {
console.log(data)
geojsonLayer.addData(data);
}

var geoJsonUrl = "http://localhost:8080/geoserver/PastDenver/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=PastDenver:dataset1&maxFeatures=50&outputFormat=text/javascript&format_options=callback:getJson"

$.ajax({
url : geoJsonUrl,
dataType : 'json',
jsonpCallback: 'getJson',
success : handleJson
});

map.addLayer(geojsonLayer);

Most of this code I took from examples from here, I have no idea what may be happening.

Does the version of JQuery matters from this AJAX? I'm using this link https://code.jquery.com/jquery-2.1.3.js


texts: How to convert a symbol size in millimeter into that in meter by PyQGIS?

texts:
Recently I am trying to use PyQGIS API to help auto-mapping in QGIS, and some symbols' overlapping on a map that is to be printed at a certain scale is not in accordance with cartography principles.
To solve this problem in a simplest way that one of two features will be "hidden" by setting its symbol size to 0 if the distance between the two features is less than the sum of the two features' symbol size in meter converted from the size in millimeter, I directly use Diameter=MaxSymbolSize*CanvasScale/1000 to do the conversion, where Diameter represents the size in meter, MaxSymbolSize=TestLayer.renderer().symbol().size(), and CanvasScale=iface.mapCanvas().scale(). This can really handle the symbol overlap problem correctly, but I guess I think there should be some errors if the mapping region is considerable large.
I try to find the related API in PyQGIS and find some functions such as convertFromMapUnits(), convertMetersToMapUnits(), and convertToMapUnits() in QgsRenderContext() class. But the definition of MapUnits and so on in QGIS is a little confusing to me and their respective return values seem the same in different map scales.
I wonder how to understand correctly about the QgsRendererContext() class and how to convert a symbol size in millimeter into that in meter by PyQGIS.


texts: Touched lines but not connected in ArcGIS 10.x

texts: Ι convert a contour file DWG to SHP file. The lines are touched but disconnected and they are broken in many places. Is there an automated way with a particular tolerance to automatically connect the touched vertex and merge the lines?


texts: Return FeatureLayer in geoprocessing REST service

texts:
I have created a Python Script Tool that returns a FeatureLayer using the following two lines of code:

arcpy.MakeFeatureLayer_management(scratch_fc, "GeocodedZIPsLayer")
arcpy.SetParameterAsText(1, "GeocodedZIPsLayer")

This works as expected in ArcMap, and resulting POINT features are added to the map without an issues.  When I publish this as a Geoprocessing Service, the service fails, and I see the following error in the server logs:

Invalid return value: GeocodedZIPsLayer

What is the best practice for returning a FeatureLayer through a REST service?


texts: How to append a geometry column to a data.frame using the R package sf

texts:
I'm working with Debian 8 usin R 3.4 and the spatial package sf 0.9.4. I try to append a  geometry column but cannot find the right iteration mode. The data table was generated from a SQL query and contains two columns for UTMX and UTMY coordinates representing the center of a cell.
head(
   tab.ds[tab.ds$rgn !='NN' & ! is.na(tab.ds$dens),
          c('irow', 'icol','ux','uy','ndist','dens', 'rgn')])

##    irow icol     ux      uy     ndist     dens rgn
## 50   10    3 316114 6028314 11134.732 1.783611   D
## 56    4    3 316958 6050121  8137.872 3.267525   D
## 57    5    3 316817 6046494  4747.973 1.747866   C
## 58    6    3 316677 6042864  2281.063 1.329282   B
## 59    7    3 316536 6039231  3694.376 3.088088   B
## 60    8    3 316395 6035595  6161.736 1.869873   C

Where the columns repersenting:

icol - Raster column index,
irow - Raster row index,
ux   - UTM-X component of the cell center,
uy   - UTM-Y component of the cell center,
ndis - the minimal distance to an obstacle,
dens - the density of an spieces and
rgn  - the regin class within the obstacle lies.

But for a test procedure the table:
 # Cell width
 RST.WIDTH = 123

 tab.ds <- data.frame( ux = c( 1,  2, -3),
                       uy = c(-1, -2, -3))

should work. I want to append a geometry column for each table row represented by the cell polygon.
To generate the polygon I've written a function which looks like this:
sf.rect.xywh <- function(x, y, w, h) {
    # using the half cellsize
    dw <- w/2;  dh <- h/2
    # create the corner points
    nw.x <- x - dw; nw.y <- y + dh;
    ne.x <- x + dw; ne.y <- y + dh;
    se.x <- x + dw; se.y <- y - dh;
    sw.x <- x - dw; sw.y <- y - dh;
    # use vectors to arrange them in a matrix
    cx <- c(nw.x, ne.x, se.x, sw.x, nw.x);
    cy <- c(nw.y, ne.y, se.y, sw.y, nw.y);
    # create the coordinate matrix
    coord <- matrix(c(cx, cy), ncol=2, nrow = 5 )
    # create the polygon
    poly <- st_polygon(list(coord))
    # give back the results
    return(poly)
}

#Test
sf.rect.xywh(10,-10, 4, 4)

## POLYGON ((8 -8, 12 -8, 12 -12, 8 -12, 8 -8))

Which works fine and gives me a closed polygon. As I understand, I've to convert the data frame into a sf object using but don't know, how to do it. If I try to append a create polygon column
  # Number of rows
  N<-nrow(tab.ds)

  # add column geom
  tab.ds$geom<-sf.rect.xywh(tab.ds$ux,
                            tab.ds$uy,
                            w=rep(RST.WIDTH,N),
                            h=rep(RST.WIDTH,N))

I get:
MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) polygons not (all) closed.

I suppose, that I need to address a row wise iteration. What is the right code to do that?


texts: Orthographic projection in Python Basemap

texts:
I'm trying to plot a GeoTIFF in Python using Basemap. I can do so correctly with cyl projection like this:

m = Basemap(projection='cyl',
            llcrnrlat=LL[1], urcrnrlat=UR[1],
            llcrnrlon=LL[0], urcrnrlon=UR[0])

But if I try to do so in the ortho projection, it does not work. I'm using this function call:

m = Basemap(projection='ortho',
            lon_0 = origin[0],
            lat_0 = origin[1],
            llcrnrlat=LL[1], urcrnrlat=UR[1],
            llcrnrlon=LL[0], urcrnrlon=UR[0])

Basemap / matplotlib crashes with the error:

AttributeError: 'Basemap' object has no attribute '_height'

If I remove the llcrnX arguments, the ortho projection works, but I see the whole globe. I guess this is OK, if I can then somehow plot my image using the corner coordinate. Basemap documentation states that the corners may be specified.


texts:
I am very new to all this coding/scripting thing and now don't know how to use https://github.com/mapbox/mbutil/blob/master/patch. It just says

$ patch [source] [dest]

but where do I put this? Python, I guess, but it says that the syntax is invalid when I put

patch C:/Temp/1.mbtiles C:/Temp/2.mbtiles

Or do I have to use another command interpreter? Execute it from a certain folder?


texts: OpenLayers 2 active select control on layer prevents drag and zoom?

texts:
I've build this web-map http://gimoya.bplaced.net/Openlayers/Trails/map.html with bike-trails (retrieved from a KML on my webserver), where several trails from a region are covered by a selectable region's polygon. Selecting a region's polygon will trigger a zoom to the polygons extent. Now, when zoomed inside the selected polygon the active select control prevents dragging and zooming.

Where in my script would I have to deactivate/activate the select control to enable dragging and zooming on top of a selected polygon?


texts: Exporting QGIS styles to SLD for GeoServer?

texts:
I´m starting to dive into sld creation for Geoserver 2.7.1

Therefore I´ve also looked at Qgis 2.8.2 as it can export styles to SLD.

Where I´m clueless is the big amount of used prefixes/namespaces.

<se:Name …
<ogc:Name …
<sld:Name …

QGIS SLD exports result in files with se: (symbology encoding).

Is it important for GeoServer which prefix/namespace is used in SLDs? If so, where can I find information and what should be noticed? Can you give a use case for when to choose which prefix/namespace?


texts: GDAL RasterIO performance in C++ vs ReadAsArray in python/numpy

texts:
I'm currently writing optimisations for a raster viewer program which uses gdal as it's base. It's currently written purely in python, and has some major speed issues which cause problems when we are reading many files at a time.
After making some optimisations in the python, and getting quite a minimal speed increase, I proceeded to profile the program quite heavily and found that our getImage method was our slowest call. I had already performed some optimisations on this function so decided to write a C-Extension so that we could get some speed increases through a lower level language.

This has worked for the most part, however there is still one issue, we have found a speed increase of ~2s for some of our larger files in the bulk of the code. But this is negated by the GDALRasterIO call, which is actually about 3s slower than the python ReadAsArray.

This doesn't make any sense to me as ReadAsArray is a wrapper around a C++ call to GDALRasterIO, and thus should be slower than having a call straight to GDALRasterIO.

I was hoping someone here might know of a way to read the rasters more efficiently. I have tried to implement a method using ReadBlock rather than RasterIO, but due to the replication that RasterIO does it didn't work at all. (I'm currently trying to figure out a way to do that replication without losing too much speed).

The RasterIO call i'm using is

band->RasterIO(GF_Read, this->ovleft, this->ovtop, this->ovxsize,
               this->ovysize, dataTmp, this->ovxsize, this->ovysize,
               band->GetRasterDataType(), 0, 0);

The old python call was:

dataTmp = band.ReadAsArray(ovleft, ovtop,
    ovxsize, ovysize,
    dspRastXSize, dspRastYSize)


texts:
I have a table of WKB HexString values, of the form:

01010000000812368FA5D106C01C25B4DAE5AE4A40

in a column called 'Centroid' in spatialite

I want to create a geometry from these values. When I run:

select GeomFromWKB(x'01010000000812368FA5D106C01C25B4DAE5AE4A40');

I get a valid geometry. However, if I run:

select GeomFromWKB('x'||"'"||Centroid||"'") from table;

I get invalid geometries. It seems that it is impossible to insert these hex values as binary values. I think may be due to Sqlite thinking that the bracked value in the second statement is a string literal and not being able to transform it into a blob/binary literal.

Does anyone have any solutions either within spatialite itself or using a third party library in python or something?

Note that I have also tried:

select GeomFromEWKB(Centroid) from table;

I seem to recall doing this successfully before in PostGIS, but I'd rather use spatialite if possible.


texts: QGIS Raster Calculator: Need values to be set to 0 instead of no data

texts:
I have a landcover map with different values for the landuse options e.g. 190 for commercially used areas. Now I want to extract only these areas using the Raster Calculator with the expression

  "Raster@1" = 190

This works fine, the resulting Raster contains value 1 for all cells which where 190 before (as it is supposed to) but no data values (high negative values) for all other cells, which I need to be zero instead. When I did the same operation to other rasters, QGIS would set the other cells to 0, which in this case is necessary for me to add the result raster to another raster later on.

Any ideas how I can avoid that QGIS sets the cells to no data cells? And what could be a reason that it deals differently with rasters of the same format when applying the same operation with the raster calculator? Both rasters are .tif but from different sources.


texts: IF (there are three or more developed cells in the Moore Neighbourhood of a cell) THEN (the cell is developed)

texts:
I need to create an expression for that principle, so far I have figured out

CON ("FocalStatisticsSum" >= 3,"FocalStatisticsSum"="FocalStatisticsSum"+1)

Although I pretty sure I am missing things, and I am not including the requirement of More neighbourhood in the expression (I don't know how to).

Also I am not sure that  >= 3 is Ok for this purpose.

Finally, how would I translate "developed"

Could somebody tell me what would be a good expression for this?

Note: Moore neigbourhood is


texts: Changing features attributes in OpenLayers

texts:
I am using OpenLayers 5.3.3 and I want to change attributes of the specific feature in my map in OpenLayers. After changing the features I want to the changes to be updated in my map.

this.vectorSource = new VectorSource({
     features: this.iconArr,
     wrapX: false,
});
this.markerVectorLayer = new VectorLayer({
     source: this.vectorSource,
});

this.map.addLayer(this.markerVectorLayer);


texts: Dynamically drawing a polygon around another shapefile, ArcGIS 10.2

texts:
Let's say I have a vegetation raster that indicates the density of forest all over the world. I also have a shapefile that contains the admin boundaries of each country in the world. The easy part is to make sure the grid cells are within the country boundary, but what I want to do afterwards is to create a separate polygon boundary that will surround the region in question. It is best to do it dynamically based on a certain density threshold, that could vary from country to country. I haven't exactly figured out how to handle that yet, but once I do, how would I go about encasing that data within a polygon that can be saved as another shapefile. I have attached an image of what I want the result to look like.

So after a threshold is determined, a polygon will be drawn around the area as tight to the boundary as possible, excuse the freehand above, and then output it as a new shapefile. Is this possible and any ideas on how it can be done? Would using Model Builder be feasible to create a tool that would take the density shapefile, boundary shapefile, and threshold value as inputs to produce a polygon?


texts: Remove certain points from point shapefile

texts:
I have a point shapefile with 9674 points.  The attribute column has three fields: X, Y and value. Out of these 9674 points, I want to remove 1700 points based on X and Y values (I know the X and Y of these 1700 points) and create a new point shapefile consisting of the remaining 7974 points.

I could easily do Sort Ascending and select the 1700 points but the problem is these 1700 points are scattered across the original shapefile and cannot not be arranged in an order.

Any suggestions?


texts:
i have an administrative unit feature dataset (polygons) with 3 levels: Admin 1 - Province, Admin 2 - District and Admin 3 - Commune.

Now i want to make a geocoding process of an address table which consists of 3 columns:Admin1, Admin2, Admin3 as same as 3 levels above. In case any record of the table have full address, it would return centroid of Admin3, if the column Admin3 is empty, it should return centroid of Admin 2 and if both of Admin3 and Admin2 are empty, the result should be centroid of Admin1.

Any ideas to do this task in ArcGIS Desktop?


texts: Converting .kml file to .tab (Map info) format

texts:
If I have saved polygons as a (.kml) Google Earth File and I want to draw it in MapInfo, is there a tool to do this ?

If not, how can be it be made using MapBasic?


texts: Calculating Distances between coordinates

texts:
I have 75 coordinates in WGS84 that I brought in from an Excel table. They are now on my imagery basemap spanning across about a half mile on either side of a dirt road in Prescott National Forest in Arizona.

My advisor wants me to get the distances between all 75 points and I first tried to use the measure tool and found that it would take too long to get all the distances that way.
I then tried Point Distance and ArcGIS just kept quitting unexpectedly.
I imagine there are at least a few more tools in ArcGIS that can do this but am not sure where to go from here.
I tried to do this in Excel as well, but I'm not familiar enough with coding to figure that out.

I am just now getting into this program and don't know very much about it.

If it helps, I'm using ArcGIS 10.2.2 on an ASUS desktop at my university. This seems to be the only PC that has ArcGIS installed on it for my use so I hope the problem is not with my computer.


texts: Landscape metrics using Google Earth Engine

texts: Has anyone developed scripts/package that calculates fragmentation metrics similar to those of FRAGSTATS using Google Earth Engine?


texts: Finding points closest to latitude longitude point using PostGIS

texts:
Right now I have a data base with PostGIS. Several of the data sets contain different locations at latitude longitude points. I will be doing several queries based at different latitude longitude points, I need to find the closest locations at each point as well as all the points within a certain number of meters. Right now I was just sorting all the data by closest to the lat/lon point that I was querying, so:

SELECT
     rest.station_id
FROM
     rest
ORDER BY rest.geom <-> ST_SetSRID(ST_MakePoint('%(lat)s',' (lon)s')::geometry, 4326)
LIMIT 1

where the lat lon variables are what I am querying.

Is there a way to do this better so that I don't have to sort each time I do a query?


texts: Import processing module from QGIS into Python Window of ArcCatalog

texts:
It is possible to import processing module from QGIS and run it in Python Window of ArcCatalog?

I tried in QGIS

 import sys
 output = sys.path

and copy output to ArcGIS Python console:

import sys
for x in output :
     sys.path.append(x)

then I tried to import:

import PyQt4.QtGui
import PyQt4.QtCore
from qgis.gui import *
from qgis.core import *
import processing

For example import processing return error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python/plugins\processing   \__init__.py", line 28, in <module>
 from processing.tools.dataobjects import *
File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python/plugins\processing\tools\dataobjects.py", line 29, in <module>
  from qgis.core import *
 File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\__init__.py", line 36, in <module>
from qgis.core import QgsFeature, QgsGeometry
  File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\core\__init__.py", line 1, in <module>
from qgis._core import *
ImportError: DLL load failed:Can`t find the specified module.


texts: Where can I find spatial data on polar fish populations? I am looking specifically for polar cod, but any fish will do at this point in my project. I also need arctic ocean salinity, temperature, and sea ice coverage


texts:
I have an ArcPy script that applies the slope function to raster files in ArcGIS. I tried running the script with two raster files in the table of contents, and the script runs, but only saves one layer and loads it into the table of contents. I think that the script is overwriting, and only saving the last file.
Here is my script:


texts: Visualizing BIG datasets with leaflet-heat

texts:
I am trying to find a fast method to visualize very BIG datasets (up to 1,6 Mio points) as a heatmap with javaScript. I tried to use leaflet-heat.js, but it is extremely slow with such a large number of points. I also tried leaflet-clustering, but on the beginning it is 10x slower as a heatmap and sometimes crashes (once the clusters are built, it is fast).

Any ideas how can a make a heatmap-layer that draws million of points fast?


texts: Export KML from points plotted in javascript Google Maps API?

texts:
I want to export a KML file from points plotted in Javascript using the Google Maps API. Is there a way to do this? The user logs in and checks Google Maps, then plot points on the map and finally export those points in a KML format.

NOTE:I am not using PHP.


texts: Keep NoData values when using gdal_merge

texts:
I'm using gdal_merge to join several BioClim rasters, which I will then feed into a Java program that "snaps" points to the nearest BioClim raster tile of non-NoData value. However, when I run gdal_merge, it assigns the NoData value of the BioClim file (NoData Value=-3.39999999999999996e+38) to each of the tiles as if it's a proper value (I think).

Running gdal_info on the merged raster shows that value as the minimum, with no indication of a NoData value; whereas running gdal_info on one of the pre-merged raster tiles shows it assigned to the NoData values.

I think this is creating an IndexOutOfBoundsException error when I then run the Java program... preventing me from snapping the points successfully.

I have two questions:

1) How do I use gdal_merge such that the NoData classification of cells are preserved in the merged output?

2) Could something else be causing the IndexOutOfBoundsException? I'm a novice with Java programs.


texts:
I am a bit new to GRASS, but I have managed to get a lot of good work done thus far. However, the following error is really throwing me off, because I can't find any help anywhere for it. I am using GRASS 7 that I downloaded today (9.23.13) on Windows 7.

Basically, I have a set of bands from a new Landsat 8 image acquired in May 2013. I want to convert these (GeoTiff format) to Top-of-Atmosphere using i.landsat.toar. According to what I've read, this should be the tool to use. I have the metadata file for the imagery and have made sure that i.landsat.toar is accessing it and that the computational region is set. However, I'm getting the following TWO errors.

Note this is the output from the tool's dialog box.

Error 1: The tool appears not to be picking up the Earth-sun distance or the Solar elevation angle from the metadata file

i.landsat.toar --overwrite --verbose input_prefix=LC80630462013143LGN01_B output_prefix=TOAR2_ metfile=C:\Users\Andrew\Desktop\All Projects\PTA\Mapping\Landsat\LC80630462013143LGN01_MTL.txt sensor=ot8

    Metada file is MTL file: new format
    RADIANCE & QUANTIZE from data of the metadata file
    ESUN evaluate from REFLECTANCE_ADDITIVE_FACTOR_BAND of the metadata file
     LANDSAT: 8 SENSOR: OLI/TIRS
     ACQUISITION DATE 2013-05-23 [production date 2013-05-28]
       Earth-sun distance    = 0.00000000
       Solar elevation angle = 0.00000000
       Atmospheric correction = UNCORRECTED

Error 2: Because of this (?) it thinks that BAND 1 (and all the rest of the bands) have calibrated DN's in the range 0 to 0.

     BAND 1  (code 1)
       calibrated digital number (DN): 0.0 to 0.0
       calibration constants (L): -0.000 to 0.000
       at-sensor radiance = -0.00000000 * DN + -0.000
       mean solar exoatmospheric irradiance (ESUN): -0.000
       at-sensor reflectance = radiance / -0.00000
    -------------------
     BAND 2  (code 2)
       calibrated digital number (DN): 0.0 to 0.0
       calibration constants (L): 330444430637036151551722492294
    587268563145980464374645729863366094385568172627644951713179
    378481117952893338397577640994585275959535639861949498028426
    262590359799651950264727525610431453748910545956564791713413
    700827628649921419448208128081385936515331098210632350384577
    822089022182894321784571544979905732238361760724622021578607
    979232195464053310102313472742663525983586004143091065701111
    514715113697333585943053423533112929011860884077135607428377
    220404697113400382954102158210112807905333778313832237152143
    633839279456924366120960750591866213481411548079987319784626
    444589368787016457003827891606484279320248476371919970747319
    49817417372212139642585088.000 to 0.000
       at-sensor radiance = -0.00000000 * DN + -0.000
       mean solar exoatmospheric irradiance (ESUN): 0.000
       at-sensor reflectance = radiance / -0.00000


texts:
I have been working on old aerial photos on Agisoft Metashape and there is a button that company added this called the "show" point residuals.
When you turn it on, it shows where the line shows.
I downloaded their software manual version from their website and have search for more information about it, but it is not on their manual

I am trying to understand this and how it works…

For example, if you have the control points on  several ungeoreferenced imagery and then you plot spread out the control points. You put down around 20 of them and you turn it on.
If you found one that has 3 different lines on three imagery and should I need to fix it?

You can see what I am drawing here is that if you see imagery 1, I move my control points from A to B. Same thing with imagery 2 and 3

Is that what am I suppose to do? or do I need to leave them alone.

Here is the attachment

if you have some links that might be very helpful.


texts: Loading in GCS shapefile into State Plane database - will it change to State plane or still be GCS?

texts: I have GCS shapefiles that I loaded back into an SDE state plane database. The database shows that they are in State Plane, not GCS. Are they actually in the correct state plane or are they still in GCS and will be distorted or shifted in some small way?


texts:
I received data saved as geodatabase from Esri (10.x) and want to open or work with it in Qgis 2.0.1. The files are like following

I read a lot in the site about this subject but I still have no idea how it could work.

Can somebody give me step by step advice?


texts: Selecting limited TOP n from specific column in ArcGIS query builder?

texts:
I'd like to know what is the best way to select a limited top 'n' of rows from a specific column of a geodatabase file. Currently I use a function in the query builder of the select tool (arcgis10) in model builder which look a lot like the one posted on another item on this site (Select top 5 most items in a geodatabase). It looks like this in the expression of the query builder:

[Join_Count] in(SELECT TOP 9  [Join_Count]  FROM spatial_join_geodat ORDER BY [Join_Count] DESC)

With this function I am able get a top n, but if the 9,10,11,12..... rows have the same value also these rows are selected. And sometimes the 9th row has a 0 as value and then all the rows are selected which I don't want. So, how do I limit the number of rows to 9 if this is the case?

Has anyone a solution for this?


texts: Iterating through values in Attribute field and performing selection using ArcPy?

texts:
I am using ArcGis Pro 2.4 python window.
I have seen suggestions of using "SearchCursor" but am unclear how to utilize this tool in the Where-Clause in my SQL selection statement to select by CPAName_X and AlcInv = 'T'.
I need help with creating a statement for my selection that iterates through each CPA Name so that I do not have to write the selection 12 times per crash type.

Here is my code so far:

import arcpy

#Set Workspace:
arcpy.env.workspace = 'N:\AGISFILE\PROJECTS\CompPlanCPAMetrics\CompPlanMetrics\CrashData_ID0612041\CrashData_ID0612041_2019_0626\CrashData_ID0612041_2019_0626.gdb'

#Input Feature Class containing the data for selection and field calculation
CrashData = "JOIN_CrashDataByCPA_2010_2016"
CPALayerCalculate = "CPAs_2016CrashesTotaled"

#Lists to Use to iterate through in selection statement
CPAName = ["Central ABQ","East Gateway", "Foothills", "Mesa del Sol", "Mid Heights", "Near Heights", "Near North Valley", "North Albuquerque", "North I-25", "Northwest Mesa", "Southwest Mesa", "West Mesa"]
CrashType = ["ALCInv", "DRUGInv","PEDInv","PECInv"]

#Select crash data By CPA from CrashData
SelectField = "CPAName" #groups all records together for selection by CPA

#Make Table View
CrashView = arcpy.MakeTableView_management(CrashData,"CrashTableView")

#SQL and Selection: Total Crashes
ListTotal = []
for i in CPAName:
    SelectStatement = "{0} = '{1}'".format(arcpy.AddFieldDelimiters(CrashData, SelectField),CPAName[i])
    SelectCPA = arcpy.SelectLayerByAttribute_management(CrashView, "New Selection", SelectStatement)
    TotalSelected = int(str(arcpy.GetCount_management(SelectCPA)))
    ListTotal.append(TotalSelected)

ListTotal1 = []
for i in CPAName:
#SQL and selection: Total Crashes By Type per CPA: add and part of sql
    SelectStatement1 = "{0} = '{1}' And {2}= 'T'".format(arcpy.AddFieldDelimiters(CrashData, SelectField),CPAName[i],CrashType[0])
    SelectCPA = arcpy.SelectLayerByAttribute_management(CrashView, "New Selection", SelectStatement1)
    TotalSelected2 = int(str(arcpy.GetCount_management(SelectCPA)))
    ListTotal1.append(TotalSelected2)
    print(ListTotal1)

Alternatively, this code iterates but does not select so I get a list with [1831,0,0,0,0,0,0,0,0,0,0,0] as a result.

ListTotalAlcInv = []
for i in CPAName:
    SQL1 = "CPAName = 'i' and ALCInv = 'T'"
    SelectCPASQL1 = arcpy.SelectLayerByAttribute_management(CrashView, "New Selection", SQL1)
    NumSelected = int(str(arcpy.GetCount_management(SelectCPASQL1)))
    print(NumSelected)
    ListTotalAlcInv.append(NumSelected)
    print(ListTotalAlcInv)


texts:
I’ve followed these instructions in the past to get this NASA hdf file (subset “PrecipitationCal”) to render in ArcMap with no problems, but now it no longer seems to work as something has changed.
Specifically, the document says to rotate the raster by -90 then warp with these coordinates:

but the image appears to be flipped and cut in half.  I'm really stumped at this point and have rearranged the order of the coordinates many different ways with no luck.


texts:
I'm super new on here and have also not much experience with Python but I really like to learn to work with this language. I'm tying to write a script for this:

I think I need a for loop instead the 'Iterate Datasets' from the ModelBuilder but how can I program it with using arcpy.ASCIIToRaster_conversion, arcpy.DefineProjection_management, arcpy.AddField_management, arcpy.CalculateField_management.
This is what I could export from the ModelBulider to a Python script:

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("Modellfunktionen")
arcpy.ImportToolbox("D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx")

# Local variables:
testdatenB = "D:\\Julia_T\\projekt\\MODIS\\karadarya\\2000_neu\\testdatenB"
File_asc =        "D:\\Julia_T\\projekt\\MODIS\\karadarya\\2000_neu\\testdatenB\\2000112_cloud_fr=e.asc"
v_Name__output_img = "D:\\Julia_T\\projekt\\erste_schritte\\zweiter_tag_a1\\ergebnisse\\test_c\\ergebnisse\\%Name%_output.img"
v_Name__output_img__2_ = "D:\\Julia_T\\projekt\\erste_schritte\\zweiter_tag_a1\\ergebnisse\\test_c\\ergebnisse\\%Name%_output.img"
v_Name__output_img__3_ = "D:\\Julia_T\\projekt\\erste_schritte\\zweiter_tag_a1\\ergebnisse\\test_c\\ergebnisse\\%Name%_output.img"

# Process: Dateien iterieren
arcpy.IterateFiles_mb(testdatenB, "*", "asc", "NOT_RECURSIVE")

# Process: ASCII in Raster
arcpy.gp.toolbox = "D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx";
# Warning: the toolbox D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx   DOES NOT have an alias.
# Please assign this toolbox an alias to avoid tool name collisions
# And replace arcpy.gp.ASCIIToRaster(...) with arcpy.ASCIIToRaster_ALIAS(...)
arcpy.gp.ASCIIToRaster(File_asc, v_Name__output_img, "INTEGER")

# Process: Projektion definieren
arcpy.gp.toolbox = "D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx";
# Warning: the toolbox D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx   DOES NOT have an alias.
# Please assign this toolbox an alias to avoid tool name collisions
# And replace arcpy.gp.DefineProjection(...) with arcpy.DefineProjection_ALIAS(...)
arcpy.gp.DefineProjection(v_Name__output_img,   "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRI    MEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")

# Process: Feld hinzufügen
arcpy.gp.toolbox = "D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx";
# Warning: the toolbox D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx     DOES NOT have an alias.
# Please assign this toolbox an alias to avoid tool name collisions
# And replace arcpy.gp.AddField(...) with arcpy.AddField_ALIAS(...)
arcpy.gp.AddField(v_Name__output_img__2_, "Prozent", "DOUBLE", "", "", "", "",     "NULLABLE", "NON_REQUIRED", "")

# Process: Feld berechnen
arcpy.gp.toolbox = "D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx";
# Warning: the toolbox D:/Julia_T/projekt/erste_schritte/zweiter_tag_a1/a1_prozent.tbx    DOES NOT have an alias.
# Please assign this toolbox an alias to avoid tool name collisions
# And replace arcpy.gp.CalculateField(...) with arcpy.CalculateField_ALIAS(...)
arcpy.gp.CalculateField(v_Name__output_img__3_, "Prozent", "([Count]/53959)*100", "VB",    "")

I thing my main problem is how to program the loop by using the last created raster in the next tool.
I have folders (15) with around 350 asc files which should turn through this loop. The result should be written in a new folder with the same name plus the new type (img). Is there something I need to do with the workspace or is it like this?

import arcpy
from arcpy import env

env.workspace = r"D:\Julia_T\projekt\MODIS\karadarya\2000_neu\testdatenB"


texts:
I'm trying to install mcc-lidar on Mac but get into problems when compiling.

[  3%] Building CXX object CMakeFiles/mcc-lidar.dir/LasDataFormat.cpp.o
/Users/Oskar/Documents/mcclidar-code/src/LasDataFormat.cpp:59:29: error: no
      member named 'str' in namespace 'boost'
        std::string scale = boost::str(boost::format("%d") % scaleDoma...
                            ~~~~~~~^
/Users/Oskar/Documents/mcclidar-code/src/LasDataFormat.cpp:59:40: error: no
      member named 'format' in namespace 'boost'
        std::string scale = boost::str(boost::format("%d") % scaleDoma...
                                       ~~~~~~~^
/Users/Oskar/Documents/mcclidar-code/src/LasDataFormat.cpp:60:33: error: no
      member named 'str' in namespace 'boost'
        std::string threshold = boost::str(boost::format("%d") % curva...
                                ~~~~~~~^
/Users/Oskar/Documents/mcclidar-code/src/LasDataFormat.cpp:60:44: error: no
      member named 'format' in namespace 'boost'
        std::string threshold = boost::str(boost::format("%d") % curva...
                                           ~~~~~~~^
4 errors generated.
make[2]: *** [CMakeFiles/mcc-lidar.dir/LasDataFormat.cpp.o] Error 1
make[1]: *** [CMakeFiles/mcc-lidar.dir/all] Error 2
make: *** [all] Error 2

Does anyone know why I get this problem? Has anyone at all managed to compile and run this app on Mac?


texts: Using ArcPy SeachCursor in Point SHP to define new variable made of point coordinates?

texts:
I have a polyline, and I am trying to extract the endpoint x,y coordinates to define them as variables.

For example, if I want to define as PtA = 551799.478099, 3762853.77587; and PtB =561687.918363, 3764343.32249 (of course these change as I loop over multiple lines.I'm having trouble understanding/utilizing the SearchCursor concept.

I was trying to dump "PtA" and "PtB" into "testAB"

import arcpy
arcpy.env.overwriteOutput = True

centerLine = r'C:\Users\X\Documents\_temp working\center_line.shp'
outpath = r'C:\Users\X\Documents\_temp working'

# Set local variables
inFeatures = centerLine
outFeatureClass = outpath + '\\' + 'centerLine_vertices.shp'
pointLocation = "BOTH_ENDS"

# Execute FeatureVerticesToPoints
arcpy.FeatureVerticesToPoints_management(inFeatures, outFeatureClass, pointLocation)

# this prints the xy coords, but how to save as variable?
for row in arcpy.da.SearchCursor(outFeatureClass, ["SHAPE@XY"]):
        # Print x,y coordinates of each point feature
        x, y = row[0]
        print("{}, {}".format(x, y))

# this doesn't work, something wrong with how I have set up "test"?
testAB = []
counter = 0
for row in arcpy.da.SearchCursor(outFeatureClass, ["SHAPE@XY"]):
    testAB[counter] = row[0]
    counter += 1

print test


texts: How can I programatically create and add features to a memory layer in QGIS 1.9?

texts:
I had a working plugin in QGIS 1.8 which read data from an MSAccess database and added it to a series of memory layers.  There is some processing involved in the interim so I don't think simply using QGIS to read directly from the database is an option.

I would like to move from QGIS 1.8 to 1.9 (primarily due to the improved quality of the print composer).  The plugin doesn't work with the new API.

I've tried a variety of methods that came up in google searches.  One, modifying the below code - from http://www.qgis.org/pyqgis-cookbook/vector.html#memory-provider, i.e. adding geometry and attributes to the dataprovider and then updating the layer - to suit the new API worked a little but the attributes were not visible until I manully entered editing mode (similar to http://hub.qgis.org/issues/3713).  An alternative approach, detailed in reply #1 of the above link, added the layer and attributes correctly but I was unable to add features to the layer.

Given that it should be a pretty simple task I'm hopeful that someone on here can offer a working example of how this should be done.
(P.S. I am not a professional programmer and most of my coding is fairly crude - I welcome any guidance but ask that you excuse some ignorance on my part)

# Receivers = a list of lists returned from a database query

# create layer
vl = QgsVectorLayer("Point", item, "memory")
pr = vl.dataProvider()

# add fields
pr.addAttributes( [ QgsField("Rec_No", QVariant.Int), QgsField("Include",  QVariant.String), QgsField("Label",  QVariant.String), QgsField("X", QVariant.Double),
                    QgsField("Y", QVariant.Double), QgsField("Z", QVariant.Double), QgsField("Height", QVariant.Double),
                    QgsField("Project_Re", QVariant.String), QgsField("NCA", QVariant.String),
                    QgsField("DayCrit", QVariant.Int), QgsField("EveCrit", QVariant.Int), QgsField("NightCrit", QVariant.Int) ] )

for i in range(len(Receivers)):
  # add a feature
  fet = QgsFeature()
  X = Receivers[i][3]
  Y = Receivers[i][4]
  fet.setGeometry( QgsGeometry.fromPoint(QgsPoint(X,Y)) )

  # Details = a list of results returned from a database query specific to each result in 'Receivers'

  if Receivers[i][3] != 0:
    Include = 'Yes'
  else:
    Include = 'No'

  fet.setAttributeMap( { 0 : QVariant(Receivers[i][0]), 1 : QVariant(Include), 2 : QVariant(Receivers[i][2]),
                         3 : QVariant(Receivers[i][3]), 4 : QVariant(Receivers[i][4]), 5 : QVariant(Receivers[i][5]), 6 : QVariant(Receivers[i][6]),
                         7 : QVariant(Details[0]), 8 : QVariant(Details[1]), 9 : QVariant(Details[2]), 10 : QVariant(Details[3]), 11 : QVariant(Details[4]) } )
  pr.addFeatures( [ fet ] )

# add a style
vl.loadNamedStyle('C:/OSGeo4W/apps/qgis/python/plugins/Gopher2QGIS/styles/Receiver_Style.qml')

# update layer's extent when new features have been added
# because change of extent in provider is not propagated to the layer
vl.commitChanges()
vl.updateExtents()
vl.updateFieldMap()

QgsMapLayerRegistry.instance().addMapLayer(vl)


texts: Finding number of categories that fall within buffer in ArcGIS Desktop?

texts:
I have two categorical datasets (polygons) and I am trying to find out how many of the categories fall within a buffer for each (one has 8 categories, the other 33).

I have used 'tabulate intersection' to get the percent of each buffer covered by each category so could manually go through and work out how many categories are in each one but as I have more than 600 points I would rather try to avoid this!


texts:
I have been trying to set up PyScripter so that I can make use of code completion for QGIS python code.  So far, I have been unable to get things to work.

PyScripter when set up as per instructions in the python cookbook is not starting.   I have set up a batch file with the following contents to launch PyScripter.

SET PYTHONHOME=C:\OSGeo4W\dev\bin
PyScripter --python27 --pythonpathdll C:\OSGeo4W\dev\bin %1 %2 %3 %4 %5
PAUSE

The next issue is registering QGIS core in PyScripter. In the PyScripter help file this is done by just adding the names of the libraries into the Tools\Options\IDE Options..Special packages area.  Some special packages such as wx and scipy are in there already but I suspect that qgis needs path specified somewhere so that correct qgis can be found.  Just adding qgis to the list of special packages and then placing from qgis import core into a py file does not result in qgis classes and method becoming available in the code explorer.

Regarding point 1, PyScripter starts and displays a list of valid command line arguments and then quits without loading.  PyScripter was installed via an installer (32 bit) and I have also tried the roaming installation based on the zip file version of PyScripter.  My goal is to get PyScripter using the same python that QGIS dev is using.  I have no idea why the command line options fail.  I am using Win7 x64 and my QGIS install is in C:\OSGeo4W\dev.

I suspect that version creep has made the existing instructions obsolete.  Any help in getting PyScripter able use qgis' python and to recognise qgis' classes would be very much appreciated.  Any set up instructions for an alternative IDE would also be welcome.

AndrewM


texts: How to create new columns in QGIS by writng Conditional query on other columns of attributes table

texts:
First of all, sorry, I couldn't phrase my question better.

So I want to be able to add a new column in my attributes table, say B, and it's value should be 1 if Column A's value is 1 and 0 if Column A's value is 2. in other words:

If A==1 : B=1

elif A==2 : B=2
I couldn't figure out a way to do this using Field Calculator.

Thanks


texts: Add layers (satellite, night mode, ...) to map server (Ubuntu)

texts: I used https://www.linuxbabe.com/ubuntu/openstreetmap-tile-server-ubuntu-20-04-osm instructions and built my offline map of my country. Now I want to add satellite map and night mode map layers to my map server. How can I do this?


texts: Run function over different raster layers depending on the grid cell values

texts:
I have a Rasterstack X of 92 layers (all days of June, July, August). Which maps the difference between the maximum temperature on that day and a threshold which I have defined.
class      : RasterStack
dimensions : 201, 464, 93264, 92  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -40.5, 75.5, 25.25, 75.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs

To give you an idea of what these layers look like I have plotted the first layer at the bottom of this post: the difference T_max - T_threshold for June 1st. What I want to do now is to take the sum over the layers, but only where X > 0 for at least 3 consecutive days. The result would be just 1 RasterLayer.
I have created the following function to do this. The script does exactly what I want it to when I run it on a vector of numbers.
total <- stack()
sum <- stack()
consec <- stack()

test <- function(X) {
  for (day in 1:nlayers(X)) {
    if (X[[day]] > 0) {
      consec = consec + 1
      sum = sum + X[[day]]
    }
    else if (consec >= 3) { total = total + sum
    consec = 0
    sum = 0}
    else {consec = 0
    sum = 0 }
  }
  return(total)
}
calc(x = X, fun = test)

This gives me the following error:
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function

I think I understand why it gives me this error: calc() usually runs a function on all layers, but I want it to run on different layers for different pixel values.
I have been breaking my head over this for a long time. Does someone know how I can best solve this?


texts: Generate non-overlapping shapes from a buffer

texts:
I have a collection of points, and I need to generate shapes around these points.

Let's assume that these points are locations of offices and that the shapes are the 'area of responsibility' around the offices. The area of responsibility is defined by airline distance. This can be easily solved by a buffer function.

But now: It is possible that offices are very close to each other and that the areas of responsibility of two or more offices overlap. In this case, only one (the closest) office is responsible.

Is it possible to generate such shapes without iterating over the positions in the overlapped area?


texts: Automation of two related datasets using ArcPy with ArcGIS Pro

texts:
I am using ArcGIS Pro.

I cannot share details or data for this question as it is for a client.

Here is a simplified version of what I am working with. I am working with a dataset that was generated by the client to be a network design for telecommunications. One feature class represents spans which are the aerial pathway between poles or underground pathway between pedestals. The other feature class represents the actual telecom cables that are to be placed that follow the pathways found in the spans table. These cables have different types, can continue across multiple spans, and more than one can be placed in a single span. In between each span is a point called pits with its own feature class, however that tables has nothing more then a Pit_ID. I don't think that this table or information will be useful for what I am looking to do, but maybe it will be.

See the example tables below:

Spans Table

Cable Table

Now what I am trying to do is fill out the fields (that will have to be added) to the Spans Table that are highlighted. There number of cables in a span don't have a predetermined limit and I have seen 10+ near the site of the cable feed source in a single span, hence the "etc." Because these cables are not directly buried but will be placed in conduits, the number of cables and types/sizes are relevant to know to determine the number of conduits needed. It is also relevant for aerial spans as only so many are allowed in a aerial span (at least without a more specialized engineering and permit.)

There maybe more than one direction this can go to arrive at the needed solution. One is more of a spatial solution since all of the cable lines overlap the span lines and share line segments. A simple Spatial join will not work since it is not a one-to-one situation. So I think that leaves some type of select by location but would have to be automated some how since each feature would have to be selected individually and then ran through the process to relate it to the other feature class. I am currently incapable of sorting this process out in an automated fashion. There maybe some other spatial route that could be done as well.

The other direction would be using just the information in the tables and I figure would be the more likely solution to be used. I don't see how a simple table join could be preformed to come to the solution either. This probably just leaves some python script to churn through the data using the two tables to fill in the needed information. With my current python skills, this too is also beyond me at the moment. I do alright using basic python in the field calculator, checking against multiple other fields in the row for a variable output using if/else and other python functions. When it comes to running arcpy functions, I am inexperienced and but looking at other examples, I can usually follow.

An added bonus for any direction would be to automatically add fields for the Cable[n]_ID and Cable[n]_type if there isn't enough blank fields for the number of cables in the individual span. I can get around this Manually adding a fields by large enough amount that I should have a few extra in the expect highest cable count in a single span.

I am also willing to export the tables to excel, working the process and then rejoining the finished excel product back with the feature class.

Given that I will have to go through this process many a times with new data-sets, I really need a semi/automatic approach to this.


texts: How do I create a geoJSON layer with a div icon instead of the leaflet standard icon?

texts:
I'm trying to load a geoJSON layer with div icons, since they are supposed to be more lightweight.

My problem is that I cannot get the divIcons to display.  The geoJSON object is stored as individualPoints.  I've tried styling the icon directly:

var myIcon = L.divIcon({className: 'my-div-icon'});

var geoJsonLayer = L.geoJson(individualPoints, {icon: myIcon});
map.addLayer(geoJsonLayer); //shows standard icon

and styling the icon on each layer:

var geoJsonLayer = L.geoJson(individualPoints, {
    onEachFeature: function (feature, layer) {
        layer.setIcon(myIcon);
    }
});
map.addLayer(geoJsonLayer); // nothing displays

and even creating an additional marker, and adding that to the map:

// create object to save markers
var markersID = [];

// exchanges latitude & longitude
function switchLatLong(oldLocation) {
    return new L.LatLng(oldLocation[1], oldLocation[0]);
}

var geoJsonLayer = L.geoJson(individualPoints, {
    onEachFeature: function (feature, layer) {
        // make a marker and add to map individually
        markersID[feature.properties.id] = L.marker(switchLatLong(feature.geometry.coordinates),{icon: myIcon});
        map.addLayer(markersID[feature.properties.id]); //nothing displays on map
    }
});

The first code snippet shows the standard leaflet icon; the second and third snippets do not display anything.

Question: How does one load a geoJson using divIcon?

Edit:
Using the new keyword on the second and third snippets displays the divIcon, but this doesn't work in the first code snippet.  Also, I would still like to have the divIcon appear without iterating through the layers to generate a new marker.  Also, why must the new keyword be used?

var geoJsonLayer = L.geoJson(individualPoints, {
    onEachFeature: function (feature, layer) {
        // set icon to a new DivIcon
        layer.setIcon(new L.DivIcon());
    }
});
map.addLayer(geoJsonLayer); // divIcon displays on map


texts: How to let the user resize a Leaflet map

texts:
I use Leaflet to display some map.
I would like the user to be able to resize the map by dragging its borders as she likes.
I would like to do this without additional Javascript libraries, instead by using the CSS resize property: https://developer.mozilla.org/en/docs/Web/CSS/resize .

For this, I added resize:both;overflow:auto; to my div#map. It now has some ugly scrollbars (a problem for later) and I can resize it by dragging the corners.

How can I make Leaflet recognise the changed size of the div now and make it update the map contents accordingly (fetch tiles, fire some event so I can request some POIs via AJAX)?

https://stackoverflow.com/questions/24412325/resizing-a-leaflet-map-on-container-resize has some hints but I am too dumb with Javascript to apply them.


texts: I want to display the production quantities for lentil by country, i.e. countries that produced more would be highlighted in a darker color. A similar example can be found here: http://faostat3.fao.org/browse/Q/QC/E


texts: Seeking shapefile of city boundaries in US?

texts:
Where can I find a shapefile that contains boundaries for the major cities in the United States?

I have tried the census website but could not find anything.


texts: Display DrillHole locations in QGIS from SQL Server database source

texts:
I'm using QGIS 3.10.9, and I want to be able to display some drill hole data from my source database, which is on a SQL server. (I've been able to export the data into CSV format and display the data correctly in the project however as I'm working in a team, it would be preferable to use the database directly rather than use an export/manage broken links to the file etc..)
To keep it simple, I've created a view that has Hole_ID, X,Y I wish to display the holeID in the project area using the x,y.
I've tried Layer > Add Layer > Add MSSQL Spatial Layer... In the pop up I have added a connection to the SQL Database.
It will only display tables when the Also list tables with no geometry check box is selected.
I then select the view that has the data (HoleID, X, Y) and click ADD.
It adds it as a layer in QGIS, but under properties I cant manipulate it to display in the project window
I've tried creating a virtual layer from this to no avail.
Can anyone point me in the right direction? It feels like I'm close but just going around in circles.


texts: Label overlaps when zooming out in QGIS3

texts:
I would like to have my label at a fixed distance to the object despite the zoom level.
The problem is, that I cannot really set the right place for it, as it's being constantly overlapped when zooming out and too far from my object when zooming in.

Looking on the E11 label above, you can see, that it's too far at the biggest zoom level and overlays the point, when zoomed out.
I found some solution here:
Making labels stay at some defined position in QGIS?
which states, that I must set XY coordinates for it, but I've already got them.

Is there any way to keep them at some defined position against the object?


texts: How can I create my own minifed openlayers js file?

texts:
I want to create my own minified openlayers js file with only the relevant classes my app used (including some custom class).

Is there a way to know by using the full openlayers.js file which are the ones I need to include in my custom build (instead doing trial and error and start deleting files and see which breaks my app) ?


texts: Uncaught TypeError: Map.addLayer is not a function earth engine javascript api

texts:
I am new to the Google Earth Engine javascript api. Sorry if my question is silly.
I have included ee_api_js.js file from the git library of Google Earth Engine javascript api. I have validated the user with the ee.data.authenticate method successfully. Now when I am trying to add layer to map it gives me "Uncaught TypeError: Map.addLayer is not a function"

var map, ee;

function initMap() {
      var latlng = {lat: 31.046051, lng: 34.851612};
      map = new google.maps.Map(document.getElementById('map'), {zoom: 8, center: latlng, disableDefaultUI: true});
    }

var runAnalysis = function() {
            ee.initialize();
var image = ee.Image('COPERNICUS/S2/20180521T081559_20180521T082244_T36SYB');
            var visParams_ndvi = {min: 0.1, max: 0.9, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
                '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};
            var image_ndvi = image.normalizedDifference(['B8','B4']);
            map.addLayer(image_ndvi,visParams_ndvi,'NDVI');
    };

    $(document).ready(function(){
        //$('#date').flatpickr();

        // Shows a button prompting the user to log in.
        var onImmediateFailed = function() {
            ee.data.authenticateViaPopup(function() {
                // If the login succeeds, hide the login button and run the analysis.
                //$('.g-sign-in').addClass('hidden');
                runAnalysis();
              });
        };

        // Attempt to authenticate using existing credentials.
        ee.data.authenticate(CLIENT_ID, runAnalysis, null, null, onImmediateFailed);
    });


texts: The Difference of clipping Images based on Gdalwarp and QGIS software?

texts:
I clipped same MODIS image by using QGIS and GDAL based on same shapefile.

Fig.1 and Fig.2 are clipped using QGIS software and GDAL coding, Fig.3 and Fig.4 are image geolocation of original image(upper red box) and clipped image(bottom red box) in QGIS and GDAL.

We can see, from Fig.3, the upper left geolocation of clipped image(bottom red box) and original image(upper red box) in QGIS are same.However,geolocation of  GDAL-clipped result(Fig.4) is on the contrary. QGIS-clipped image include some background pixels outside of shapefile boundary(Fig. 1 compared to Fig. 2).

I want to get the result as Fig.2 using GDAL or QGIS.

Fig.5 is clipping command in GDAL,Fig.6 is clipping command in QGIS.
Why exits result difference?

              Fig.1

              Fig.2

              Fig.3

              Fig.4

              Fig.5

              Fig.6


texts: How to improve performance of Oracle Spatial function sdo_aggr_mbr

texts:
I experienced bad performance when using the sdo_aggr_mbr function on a data set of more than 1 million rows, I had to terminate the operation after several minutes.

The statement was:

select sdo_aggr_mbr(MY_GEOMETRY) from MY_TABLE

When I execute a similar operation, to create the bounding box in a manual way, the performance is very good:

select min(t.MY_GEOMETRY.sdo_point.x),
max( t.MY_GEOMETRY.sdo_point.x),
min(t.MY_GEOMETRY.sdo_point.x),
max(t.MY_GEOMETRY.sdo_point.y)
from MY_TABLE t;

That indicates that the database and the spatial index of the geometry column is configured in a correct way.

Why is the aggregation function of oracle is so slow and what can be done to improve the performance?


texts: ArcGIS Server 10.1 Map Service Not Using Maximum Instances

texts:
I'm running a stress test tool that's simulating 10 concurrent requests for an ArcGIS Server 10.1 (Windows) Map Service.  I have the Map Server set to use a maximum pool of 4 instances.

I'm not using caching - these are full screen map images based on random coordinates within a bounding area of interest.

Regardless of whether I'm using File Geodatabase, PostGIS or SQL Server as the source of the data - the maximum number of instances that gets used is 2.

FYI - the requests are coming from separate threads and are staggered by 200ms with an in-between request delay of 1000ms. So I would expect it to fire on all cylinders during the test.

Any ideas why it never uses the maximum number of instances?


texts: Remove ArcSDE repository and associated files from PostGresql DB

texts:
A bit of background first.

I have a PostGres DB with PostGIS, I also have ArcSDE installed in this DB and all feature classes created through ArcSDE are using the PG_Geometry. Everything works well!

Over the next twelve months I will be looking at getting all of my desktop and web clients away from esri products to open source equivalents which for me means i'll be free of ArcSDE at last (I'm not an esri basher, I enjoy working with esri products but I can't justify the existance of ArcSDE doing something that should be free). Anyway I want to know if there is a way which I can "uninstall" ArcSDE from PostGres? I know that if the ArcSDE tables are left there users wont notice a difference but from an admin perspective it would be much nicer and tidier.

Has anyone done this in the past? Any tips or experiences?

Thanks
Ando


texts: Trying to find the manual for a Sokkia GIR1000

texts: I'm trying to track down a manual for a Sokkia GIR1000.  My department has one but over the years the manual seems to have disappeared.  Does anyone have one or know where I can get one?


texts: Where is the field mapping for Attribute Transfer tool stored?

texts:
I can't seem to find any patterns in when the source-destination field mapping is preserved. Is this stored in the database, map document, or elsewhere? I'd like to know so I don't have to create the mapping every time I begin work.

Dialog for reference below.


texts: Does QGIS handle AREA_OR_POINT metadata tag?

texts:
I used to think that Sentinel-2 L1C images stored on AWS are AREA_OR_POINT=Point, i.e. the origin lies in the centre of the pixel. Please, correct me if I'm wrong. For convenience I convert jpeg2000 to geotiff with gdal_translate, but it sets AREA_OR_POINT=Area by default. So I change the tag with gdal_edit. When I open the images as different layers the origin of the second image is not shifted by a half of the pixel size. Does QGIS handle AREA_OR_POINT tag?
Here are the steps to reproduce:

    # download an image
    wget --quiet http://sentinel-s2-l1c.s3.amazonaws.com/tiles/15/S/TD/2016/8/20/0/B03.jp2
    # convert to geotiff
    gdal_translate -of GTiff B03.jp2 B03.tif
    # make a copy
    cp B03.tif B03_point.tif
    # change the tag
    gdal_edit.py -mo AREA_OR_POINT=Point B03_point.tif
    # now open B03.tif and B03_point.tif in QGIS, the images are displayed as identical

I'm using QGIS 2.18.12 compiled against GDAL 2.1.0


texts: Placing value on top of bar chart thematic map using MapInfo?

texts:
Is it possible to place the value on top of each bar in a bar chart thematic map in Mapinfo?

This is how I did it so far: I wanted to show a value per zone (the amount of trips generated in one area) with a bar: the areas which generate more trips have a larger bar. I also wanted to place the value of such bar on top, but as I couldn't find how, I placed the label at the base of the bar by choosing the following label display options:

This is what the thematic map looks like:

I would like the numbers to be on top of the bars, if it is possible, especially because for my next thematic map I will have to represent several fields at a time (origins and destinations, that is: two bars), so I will no longer be able to use labels to show the value.


texts: Trigger geocoding without having to visit page

texts:
I'm using Google's geocoder to store results in database then display them on a map on another page. The problem is that each company profile page needs to be visited before it's geocoded.

Is there a way to run the geocoder if a company address's lat/lng is 0 (i.e it hasn't been geocoded)?

<?php

require("database.php");
// Opens a connection to a MySQL server
$con = mysql_connect("localhost", $username, $password);

if (!$con)
{
    die('Could not connect: ' . mysql_error());
}

mysql_select_db("medicom_wp", $con);

    $company = get_field('company_name');
    $address = get_field('address');
    $city = get_field('city');
    $post_code = get_field('post_code');

    $sql = sprintf("select count('x') as cnt from markers where `name` = '%s'", mysql_real_escape_string($company));
    $row_dup = mysql_fetch_assoc(mysql_query($sql,$con));
    if ($row_dup['cnt'] == 0) {
        mysql_query("INSERT INTO markers (`name`, `address`, `lat`, `lng`, `type`) VALUES ('".$company."', '".$address.", ".$city.", ".$post_code."', '0.0', '0.0', '')");
}
wp_reset_query();

define("MAPS_HOST", "maps.google.com");
define("KEY", "");

// Opens a connection to a MySQL server
$connection = mysql_connect("localhost", $username, $password);
if (!$connection) {
  die("Not connected : " . mysql_error());
}

// Set the active MySQL database
$db_selected = mysql_select_db($database, $connection);
if (!$db_selected) {
  die("Can\'t use db : " . mysql_error());
}

// Select all the rows in the markers table
$query = "SELECT * FROM markers WHERE 1";
$result = mysql_query($query);
if (!$result) {
  die("Invalid query: " . mysql_error());
}

// Initialize delay in geocode speed
$delay = 0;
$base_url = "http://" . MAPS_HOST . "/maps/geo?output=xml" . "&key=" . KEY;

// Iterate through the rows, geocoding each address
while ($row = @mysql_fetch_assoc($result)) {
  $geocode_pending = true;

  while ($geocode_pending) {
    $address = $row["address"];
    $id = $row["id"];
    $request_url = $base_url . "&q=" . urlencode($address);
    $xml = simplexml_load_file($request_url) or die("url not loading");

    $status = $xml->Response->Status->code;
    if (strcmp($status, "200") == 0) {
      // Successful geocode
      $geocode_pending = false;
      $coordinates = $xml->Response->Placemark->Point->coordinates;
      $coordinatesSplit = split(",", $coordinates);
      // Format: Longitude, Latitude, Altitude
      $lat = $coordinatesSplit[1];
      $lng = $coordinatesSplit[0];

      $query = sprintf("UPDATE markers " .
             " SET lat = '%s', lng = '%s' " .
             " WHERE id = '%s' LIMIT 1;",
             mysql_real_escape_string($lat),
             mysql_real_escape_string($lng),
             mysql_real_escape_string($id));
      $update_result = mysql_query($query);
      if (!$update_result) {
        die("Invalid query: " . mysql_error());
      }
    } else if (strcmp($status, "620") == 0) {
      // sent geocodes too fast
      $delay += 1000;
    } else {
      // failure to geocode
      $geocode_pending = false;
      echo "Address " . $address . " failed to geocoded. ";
      echo "Received status " . $status . "
\n";
    }
    usleep($delay);
  }
}

?>

Before anyone says, I actually followed a tutorial from Google themselves saying how to store geocode results in a database, so in my eyes it doesn't go against their TOS because they're telling you to do it!


texts: How to resolved ORA-12154: TNS:could not resolve the connect identifier specified..?

texts: When application try to connect to an Oracle Service via Oracle Networking SQL*Net, it looks up the connect data for the service. When this lookup fails, Oracle client show an error message ORA-12154.


texts:
I have these lines (blues lines) that extract from a raster and need to unite them automatically (red lines).
Currently I'm doing it manually using a plugin (Join multiple lines) but it's very tiring. Is there any way to grow and join the lines automatically? I thought of joining based on the angle of inclination but I do not know how to proceed.
I use QGIS and Ecognition.

What I have:

What I want:


texts: Adding support to Esri's FileGDB driver when fiona/gdal/osgeo were installed through Anaconda

texts:
Is there a way to tweak the fiona/gdal/osgeo packages installed through Anaconda so that they can be used to read/write GDB files using the ESRI FileGDB driver? The main thing I'm actually trying to achieve is writing to GDB a file.
I have found multiple answers online, but all of them involve re-installing OSGeo or GDAL from the python wheel files. These methods do not work for me because I installed fiona, OSGeo and GDAL using Anaconda. So none the solutions I found online (listed below) seem to work.
I have already downloaded the DLL file from ESRI, placed it in my C:\...\Lib\site-packages\osgeo folder and added a new environment variable called GDAL_DRIVER_PATH indicating the appropriate path. However, contrary to the suggestions in the posts listed below, the __init__.py inside the osgeo folder does not have a commented line that reads os.environ['GDAL_DRIVER_PATH'] = os.path.join(os.path.dirname(__file__), 'gdalplugins'). Therefore, I am not sure where to include this line in my own file system.
Similar (unresolved) discussions can be found on the GeoPandas and on the Fiona Github pages.
Also, here is a list of related posts/questions/sites with relevant info but which ultimately did not solve my problem.

How to add support for FileGDB (Esri file gdb API) driver in
fiona?
How to edit file geodatabase in Python using FileGDB Driver?
Geopandas: Write layer back into GeoDataBase
Exporting PostGIS Data
Setting up GDAL/OGR with FileGDB Driver for Python on Windows
GDAL FileGDB Driver with Anaconda / ArcGIS Pro (Unanswered)
GDAL Documentation on using ESRI File Geodatabase (FileGDB)
How to access feature classes in file geodatabases with Python and GDAL?
ESRI FileGDB API on GitHub

How I installed Fiona and other packages:
I'm currently using Windows 10. After installing Anaconda, I created a new environment called myenv and installed the main libraries I was going to use:
conda create -n myenv pandas geopandas fiona spyder numpy scipy

GDAL was automatically installed as one of the associated/required packages for fiona and geopandas.
In my command prompt, when I activate myenv and type gdalinfo --version, I get the following response:
GDAL 3.0.2, released 2019/10/28

When I use the python-based route, I get the same answer.
import osgeo.gdal
print(osgeo.gdal.__version__)
>3.0.2

Additional Notes:
If the answer involves something like "you need to compile/build GDAL from source files", please include a step-by-step on how to do so, including what to do after "building" (for example, where to put the files and how to get Python to use the newly-built version of GDAL within a specific conda environment). I know this is a lot to ask, and I apologize for doing so. It's just that I have tried looking into this process and it is not very newbie-friendly.