texts:
I am having trouble deleting the features I have selected using 'arcpy.SelectLayerByLocation_management' tool.  To use the select layer by location tool, you much pass in a .lyr file.  My script SHOULD create a .lyr file of the input data, select by location on the .lyr file, make a copy of the selected features and save to a .shp file, and then delete the selected features in my input data.  Instead, it creates a .lyr file, selects features in the layer file and makes a copy of the layer file (correct up until this point).  Then it deletes the .lyr file and im stuck at with an unchanged input data set.  Should i make a copy of the .lyr file after i have used the arcpy.DeleteFeatures_management tool?

I guess the question in fewer words is:

When i call the arcpy.DeleteFeatures_management tool on a .lyr file that has selected features, does it delete the entire .lyr file or just the selected features?

Here is my code:

def dBaseComparison(self):
        '''This function compairs the input data to the database and selects
            the features that already exist'''
        try:
            # Create a feature layer to select the points that already exist in      the database.
            arcpy.MakeFeatureLayer_management(self.MakeXYeventLayer(), 'PointInFile.lyr')
            arcpy.SelectLayerByLocation_management('PointInFile.lyr', "ARE_IDENTICAL_TO", self.dBase)
            # Make a copy of the existing features and delete them from the new file.
            arcpy.CopyFeatures_management('PointInFile.lyr', self.MakeXYeventLayer()[:-4] + "_selected.shp")
            arcpy.DeleteFeatures_management('PointInFile.lyr')
            arcpy.CopyFeatures_management('PointInFile.lyr', self.MakeXYeventLayer()[:-4] + "_NewData.shp")
            return self.MakeXYeventLayer()

        except arcpy.ExecuteError:
            return arcpy.GetMessages( )

Here is a stand alone script that works great! what is the difference?:

import arcpy
arcpy.env.workspace = 'C:\Temp\Project_Output'
arcpy.env.overwriteOutput = 1

dataBase = "C:\Temp\Project_Output1\dBase.shp"
infile = "C:\Temp\Project_Output1\NC_CSV_TEST_Point.shp"

# Create a feature layer to select the points that already exist in the database.
arcpy.MakeFeatureLayer_management(infile, 'PointInFile.lyr')
arcpy.SelectLayerByLocation_management('PointInFile.lyr', "ARE_IDENTICAL_TO", dataBase)
# Make a copy of the existing features and delete them from the new file.
arcpy.CopyFeatures_management('PointInFile.lyr', infile[:-4] + "_selected.shp")
arcpy.DeleteFeatures_management('PointInFile.lyr')
arcpy.CopyFeatures_management('PointInFile.lyr', infile[:-4] + "_NewData.shp")


texts: How to convert point data to different sized polygons representing Marine protected areas?

texts:
I am trying to map Marine protected areas. I have collected data of where all these areas are. The data is expressed as point data (seen within the top image). Within the attribute tabled I have a column heading named Area. This identifies the area of the areas situated around each point of data. The areas vary is size and shape, therefore the buffer tool, wouldn't be very useful (additionally my CRS is within degrees not meters, so would be inaccurate to do).

I need a map which looks like the bottom image, So I was wondering if there is any way to use my area data to convert the points into polygons?


texts:
I was following this tutorial: http://www.acgeospatial.co.uk/k-means-sentinel-2-python/ on how to generate an unsupervised classification for sentinel-2 images.

img_ds = gdal.Open('data/products/sentinel_bands.tif', gdal.GA_ReadOnly)
band = img_ds.GetRasterBand(2)
img = band.ReadAsArray()
X = img.reshape((-1,1))
k_means = cluster.KMeans(n_clusters=8)
k_means.fit(X)
X_cluster = k_means.labels_
X_cluster = X_cluster.reshape(img.shape)

This generates the classification as a numpy array

Now, I would like to transform that array to a shapefile.
From my research I saw that you have to transform the array to a geodataframe with geopandas.
I was trying to do that here:

X_cluster_dataframe = pd.DataFrame(data=X_cluster[1:, 1:],
                                   index=X_cluster[1:, 0],
                                   columns=X_cluster[0, 1:])
X _cluster_geodataframe = gp.GeoDataFrame(X_cluster_dataframe,
                                          geometry=gp.points_from_xy())

But I do not know how to specify the geodataframe's geometry.

I know that there are tools that can transform rasters to shapefiles, but the idea is to implement an efficient way to classificate big sentinel-2 imagenery and it would be better if I can generate the shapefile for the classification without fisrt having to create the raster.


texts: gdalwarp causes ERROR 6: PNG driver doesn't support 6 bands in gdal2tiles

texts:
I've a problem with my gdalwarp function.

What I've done:

I built a .vrt file  from 6 .tif images:

gdalbuildvrt -overwrite /input/base.vrt /home/tiles/BASE/*.tif

My gdalinfo looks like that:

Driver: VRT/Virtual Raster
Files: temp/vrt/base.vrt
       /home/tiles/BASE/17_BASE_1.tif
       /home/tiles/BASE/17_BASE_2.tif
       /home/tiles/BASE/17_BASE_3.tif
       /home/tiles/BASE/17_BASE_4.tif
       /home/tiles/BASE/17_BASE_5.tif
       /home/tiles/BASE/17_BASE_6.tif
Size is 60007, 40007
Coordinate System is `'
Origin = (468000.000000000000000,5510000.000000000000000)
Pixel Size = (0.100000000000000,-0.100000000000000)
Corner Coordinates:
Upper Left  (  468000.000, 5510000.000)
Lower Left  (  468000.000, 5505999.300)
Upper Right (  474000.700, 5510000.000)
Lower Right (  474000.700, 5505999.300)
Center      (  471000.350, 5507999.650)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
Band 4 Block=128x128 Type=Byte, ColorInterp=Undefined

Everything works fine until now.
What I've done next is building one big .tif image from the .vrt file:

gdalwarp -multi -r cubic -s_srs EPSG:25832 -t_srs EPSG:3857 -cutline /shp/17.shp -crop_to_cutline -dstalpha -co compress=LZW /input/base.vrt /output/base.tif

But, now the gdalinfo from the output tif looks like:

Driver: GTiff/GeoTIFF
Files: temp/base.tif
Size is 47849, 25603
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (952949.608269182615913,6400383.903813457116485)
Pixel Size = (0.154577105681076,-0.154574816517696)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  952949.608, 6400383.904) (  8d33'37.77"E, 49d44' 3.40"N)
Lower Left  (  952949.608, 6396426.325) (  8d33'37.77"E, 49d42'40.66"N)
Upper Right (  960345.968, 6400383.904) (  8d37'36.96"E, 49d44' 3.40"N)
Lower Right (  960345.968, 6396426.325) (  8d37'36.96"E, 49d42'40.66"N)
Center      (  956647.788, 6398405.114) (  8d35'37.37"E, 49d43'22.03"N)
Band 1 Block=47849x1 Type=Byte, ColorInterp=Gray
Band 2 Block=47849x1 Type=Byte, ColorInterp=Undefined
Band 3 Block=47849x1 Type=Byte, ColorInterp=Undefined
Band 4 Block=47849x1 Type=Byte, ColorInterp=Undefined
Band 5 Block=47849x1 Type=Byte, ColorInterp=Alpha

And those Bands causes an error in the gdal2tiles function:

ERROR 6: PNG driver doesn't support 6 bands.  Must be 1 (grey),
2 (grey+alpha), 3 (rgb) or 4 (rgba) bands.

Does someone has any idea what happens? Why is the Band 2-4 Undefinied and Band 1 Gray? I think it should like the gdalinfo from the .vrt file (=> RGBA)?


texts: Connecting polygon accidentally deleted - filling in the gap

texts:
New QGIS user, Brighton 2.6.   I accidentally deleted an elongated (lots of "peninsulas" polygon that represents a stream-based habitat - two gaps on SW boundary and two on the northern boundary (leaving a small 'island' or floating polgyon.

It seems that one of the tools would allow quick closure of this area to draw back into as a record.   There are four gaps at the formerly contiguous boundary, but if I could sketch a line across these gaps and close the area back in, that would be ideal.


texts: Merging two Polygon objects of the same layers, snapping them together into one?

texts:
Currently I am converting an AutoCAD Draw project (DXF file) to QGIS where I have some several shapefiles as skeletons to combine with data. (Cut from one layer and paste them into a another. I have converted polylines (which are contours) to polygons. This is because "normally" these polylines are supposed to be closed. However, I encountered a problem. Two specific polygons which are supposed to be the same one, have not been converted correctly for unknown reasons.

How it looks in QGIS the two polygons are selected in yellow.

This what is supposed to look like. The black lines is suppose to be closed and go around the square object, this is what was not drawn that way. But you will get the point.

How do I merge them together?


texts: How to change the vertical datum of a geotiff DEM

texts: I am trying to do some change detection between SFM derived DEM products, but the DEM i would like to use as my base elevation (and one other) is in NAD83, the remainder are WGS84.  I would like to reproject the raster values (elevations) to WGS84 datum. Arcmap 10.4 cannot do vertical datum shifts, and VDatum does not seem to be changing the values within the raster when asked to do a datum shift, because it does not 'know' that the values are elevation.  Can I use GDAL or rasterio to tackle this problem?  I am new to geospatial analysis in python, but I am familiar with the language .


texts: Finding qgis2 folder to delete in ubuntu?

texts:
How can i find the qgis2 in qgis and where it is installed, i had one error on when qgis open this error will occurred "[ERROR] Can not make qgis.db private copy" some links told delete the folder qgis2, i can't find this folder.

How do I solve this problem?


texts: Filter polygons that share more than 2 vertices

texts:
Is  there a way to filter polygons that share more than 2 vertices with a feature?
Suppose I have a bunch of polygons, like in the image, but I just wanna see those that touches more than one vertices of the main feature.
I can only go till this point with filter rule in symbology, working with touches($geometry). If somebody could help me I would be very grateful.


texts: Seeking detailed instructions on how to install osmosis on MacOS 10.12

texts:
I want to extract some features from OSM, and gather that the best way to do this is by using Osmosis.
Can anyone provide me with instructions on how to install Osmosis for mac os 10.12?
so far I have dowloaded the osmosis zip file from here https://github.com/openstreetmap/osmosis/releases/tag/0.48.3
I have tried to follow this information:
https://wiki.openstreetmap.org/wiki/Osmosis/Installation#OS_X
however I do not understand the following instructions
" once a symbolic link to the osmosis executable has been created into a directory in the PATH.
Example:
$ PATH=$PATH:~/bin
$ ln -s ~/opt/osmosis/bin/osmosis ~/bin
The first line should be put into ~/.bash_profile configuration file, so that the addition of directory ~/bin to the PATH becomes permanent."
I'm not sure how to make the changes it suggests to the PATH? or even how to find the PATH...
apologies for the basic request, but I'm not experienced. I have tried looking for other links with alternative info but didn't have much luck.


texts: How to debug Quantum GIS PlugIn Visual Studio DLL?

texts:
I have developed a Quantum GIS plugin layer using Microsoft Visual studio and placed it in my plugins folder of my customized QGIS application.

I now intend to debug the DLL using Microsoft Studio.

I have Quantum GIS 1.8 sources and debug successfully compile and also try to debug my dll using this environment without success.

My question is: How can i debug a Quantum GIS DLL using Microsoft Studio?


texts:
I've been using QGIS 2.10 for a few weeks and have just started trying to export maps using Map Composer. This might be a very common question with a trivial solution but I just can't figure it out. For some reason, the map grid and coordinates aren't correct, and I cannot find anyway to fix them.

In the image below, I have created a plot of an example bathymetry raster with a grid that is using map units (both project and raster are in UTM 24S CRS), and an interval of "5" for both X and Y.

I have also inserted a distance scale with "1 maps units per bar unit" and a fixed width of "50 map units".

I have measured objects using this scale and it seems to be correct. Unfortunately, it completely disagrees with the map grid coordinates. I am very connfused with this, and mystified by how I can fix this.

Also notice that the grid coordinates start at zero on the lower-left corner. This is also completely wrong. Also don't know how to fix.


texts: Can I make PostGIS treat polygon edges like Google Earth does?

texts:
Consider the following geographic polygon in PostGIS:

POLYGON((60 -60,60 -50,179 -50,179 -60,60 -60))

The smallest latitude is -60, so initially I expected any point with a latitude less than -60 to fall outside of the polygon.

But it doesn't; a point at 136,-73 is considered to fall within the polygon:

SELECT ST_Covers(ST_GeogFromText('SRID=4326;POLYGON((60 -60,60 -50,179 -50,179 -60,60 -60))'), ST_GeogFromText('SRID=4326;POINT(136 -70)'));

returns true.

After some thinking, I realised that this is because PostGIS uses great circles to form the edges of polygons, and in this case, that makes the lower edge of the polygon 'bulge out' towards the south pole.

If I view the same polygon in Google Earth, the lower edge of the polygon runs along the -60 degree parallel, so apparently Google Earth uses a different manner to connect vertices than PostGIS.

What is the name of the manner that Google Earth uses, and can I somehow make PostGIS use that same manner?


texts: Excluding school-areas using Overpass-Turbo

texts:
As an example see: http://overpass-turbo.eu/s/PV6

How do I exclude the gyms on school-property?

I tried -relation["amenity"="school"] but I get

parse error: Unknown type "-"
static error: For the attribute "type" of the element "query" the only allowed values are "node", "way", "relation" or "area".
parse error: ';' expected - 'relation' found.


texts: Is there a python module to create Virtual Raster (.vrt) xml files? Rather than writing text strings to a file, I'd like to know if there is an automated way to do this.


texts: Editing attributes of a GeoJSON in Leaflet

texts:
I have a map in which I'm trying to create the capabilities of adding new features and editing attributes of new and existing features. I have managed to make the feature editable but I want people to be able to add new points to the same GeoJSON and edit attributes but I am completely stuck. I'm new to using leaflet and GeoJSON so I'm having a hard time figuring out the issue. This is the code I have right now for editing features:

 var drawnItems = businessLayer;

 map.addLayer(drawnItems);

 var drawControl = new L.Control.Draw({

     edit: {
         featureGroup: drawnItems
     }
 });
   map.addControl(drawControl);

I tried using solutions find at the following link:
https://stackoverflow.com/questions/32834401/edit-feature-attributes-in-leaflet

But I could not get it to work. I can post the rest of my code as well if that will help.

Anybody have any helpful links or insight that can help me resolve this issue?


texts: Error: Geometry constraint geom-type or SRID not allowed. How to ensure the correct geom-type?

texts:
Following the top answer on a previous question about storing geometries in a spatialite database, I am trying to update a spatialite table using:

engine.execute("SELECT AddGeometryColumn('your_table', 'geog', 27700, 'polygon', 'XY');")
# add a Spatialite geometry column called 'geog' to the table, using ESPG 27700, data type POLYGON and 2 dimensions
    # (x, y)
engine.execute("UPDATE your_table SET geog = GeomFromWKB(geometry, 27700);")

The code worked fine when I was using point data. However, when I use vector data I keep getting the error
your_table.geog violates Geometry constraint [geom-type or SRID not allowed]
I know the SRID is ok, so from the error I assume it is the geom-type that is incorrect. I have checked the geom-type in a couple ways: 1. opening the layer in QGIS and checking the porperties; and 2. opening the database in spatialite-gui and 'check geometries'. I have also tried using polygon, multipolygon and multipolygonz as the 4th argument in the AddGeometryColumn selection, but I get the same error. So, what is the best way to check which geometry it should be?


texts: Automating transfer of symbology information of raster to polygon feature using ArcMap

texts:
I have a raster dataset that is symbolized using Natural Breaks (Jenks) classification (9 classes). To reduce file size for web service display and download, I used the Reclassify tool to remap pixel values into 9 classes (class values were 1-9), and then used the Raster to Polygon tool to create a vector version of the dataset.

It was important to label polygon features based on the original raster symbology classification. So instead of a polygon labeled as 1, it would be labeled 10.3-24.2. To do this, the Remap Table was saved as an INFO table (functionality provided in the Reclassify tool window - see table below) and later joined to the polygon attributes table, in order to create a labeling field based on class range values.

This process works fine by manually saving the Remap Table as an INFO table (click save button in the Reclassify tool window), but i need to automate the process, since i need to create about 100 of these maps.

Is there a way to automate the transfer of symbology information (class range values) of a raster to a polygon feature?


texts: Which events are triggered when I draw, for example, a polygon? Like when I define the vertices of the polygon by clicking? And how can I capture these events? Or, better: Is there any documentation in the web that covers these issues?


texts: Converting WKT to WKB in Java using PostGIS?

texts:
How can I convert a WKT string to a WKB byte[] in Java using PostGIS functions?
I need to save that WKT string in a PostGIS database (geometry column.


texts:
I am trying to use osm2pgsql on Windows to import OSM data to a PostGIS database. I want to pass the EPSG:27700 via the --proj option as follows:

osm2pgsql --proj 27700 <other options not shown>

However, I get the following error:

osm2pgsql version 0.96.0 (64 bit id space)

Osm2pgsql failed due to ERROR: creation of CRS failed: No such file or directory

Any ideas as to why this is not working?

I used the Windows binary from osm2pgsql | History


texts:
I am trying to create a custom CRS in QGIS and the proj parameters given are as follows:

proj +ellps=sphere +a=127.4 +e=0 +proj=stere +lat_0=90 +lon_0=-32 +lat_ts=60 +x_0=8 +y_0=110

I used these parameters in Settings > Custom CRS in QGIS to create the Custom CRS. But when I created a grid using the new CRS it says 'CRS undefined'. What is the correct way to define this CRS in QGIS?


texts: dijkstra_sp_directed() not returning valid paths

texts:
I think I have found a major issue with pgrouting. I hope I am just doing something wrong due to my inexperience.

I am plotting the following query using OpenLayers.

select transform(the_geom, 900913) from dijkstra_sp_directed('ways', 52343, 39219, true, true)

On first glance this appears to work perfectly.

-->http://dl.dropbox.com/u/11502389/Screenshot-1.png

However if I successively label the nodes in each way, I learn that the path is not defined as I would expect.
I expected the below output.

Way 0: 1, 2, ...n
Way 1: 1, 2, ...n
.
.
.
Way n: 1, 2, ...n

What I actually get is a mess, some ways are ordered 1, 2, ...n and some are ordered in reverse n,... 2, 1

Here is a zoomed screenshot of the start of the path, this should be labeled 1, 2, 3, 4, from the screenshot you can see that its backwards.

-->http://dl.dropbox.com/u/11502389/Screenshot.png

Here is a zoomed screenshot of the end of the path, this is labeled correctly.

-->http://dl.dropbox.com/u/11502389/JourneyEnd.png

Why are some of the ways reversed? As far as I am concerned this is not a valid path. If a person was to follow this path
(1, 2, 3), (5, 4, 3, 2, 1) they would have to teleport from the end of the first way to the end of the second way and then walk backwards towards the end of the first way.

More information

The query I am using to label the nodes is below

select (ST_DumpPoints(transform(linemerge(the_geom), 900913))).geom, (ST_DumpPoints(transform(linemerge(the_geom), 900913))).path from dijkstra_sp_directed('ways', 52343, 39219, true, true)

The data is OSM data for Ireland imported using osm2pgrouting
and I followed this blog post to import the data.

Any help with this issue would be greatly appreciated.

Here is my attempt to modify the function dijkstra_sp_directed. I think my logic is correct. I must be doing something that SQL does not like because when I run it I get the below error

ERROR: record "previous_path_result" is not assigned yet

DECLARE
    r record;
    path_result record;
    v_id integer;
    e_id integer;
    geom geoms;
    query text;
    id integer;

    previous_path_result record;
    temp integer;

BEGIN
    id :=0;

    query := 'SELECT gid,the_geom,source,target FROM ' ||
          'shortest_path(''SELECT gid as id, source::integer, target::integer, ' ||
          'length::double precision as cost ';

    IF rc THEN query := query || ', reverse_cost ';
    END IF;

    query := query || 'FROM ' ||  quote_ident(geom_table) || ''', ' || quote_literal(source) ||
          ' , ' || quote_literal(target) || ' , '''||text(dir)||''', '''||text(rc)||'''), ' ||
          quote_ident(geom_table) || ' where edge_id = gid ';

    FOR path_result IN EXECUTE query
        LOOP
             geom.gid      := path_result.gid;
             geom.the_geom := path_result.the_geom;
             id            := id+1;
             geom.id       := id;

             IF NOT previous_path_result IS NULL THEN
                 IF NOT path_result.source = previous_path_result.target THEN
                     temp := previous_path_result.source;
                     previous_path_result.source := previous_path_result.target;
                     previous_path_result.target := temp;
                     geom.the_geom := ST_Reverse(geom.the_geom);
                 END IF;
             END IF;

             previous_path_result := path_result;
             RETURN NEXT geom;

        END LOOP;
    RETURN;
END;


texts: Importing Lat/Long coordinates as CSV file & only displaying single point in QGIS

texts: I am using QGIS and trying to import Lat/Long coordinates as delimited text layer. All points show up in my attribute table, however, only one point is displayed on the map? Any ideas?


texts:
I am trying to label polygons but it does not always work.

Either the polygon is quite large, then the label appears, or the polygon is small, then the label is not displayed.

I've already tried to lower label size, disable box about overlap labels, disable scale limits, yet, nothing works.

If I want the label being displayed, I have to use the node tool in order to make the polygon larger ...  this is not the good solution.

Has someone already met this problem?
(I have looked in other questions, but only polygons geometry validity is mentioned, and I have already checked it)


texts: Where to find datasets for creating CARTO VL maps?

texts: Is there any way to download example datasets like NYC taxi and origin-destination data? Like those that show flight paths from and to some cities. I think it'd help me greatly if I could somehow get those data and try to create CARTO VL animations on my own.


texts: ArcPy - Symbology Color Changes After Setting Class Break Values

texts:
I develop a python script tool to adapt the symbology of a set of shapefiles based on a given symbology template (layer file).

These are the properties that should be changed:

Number of classes and colors: using the ApplySymbologyFromLayer() function.
The symbology template has 4 classes set to natural breaks and the colors are set manually to specific RGB values (green, brown, red, blue), there is no color scheme used.

Class break values and labels: calculated based on statistical values of each shapefile (in the code below entered arbitrary values)

After using the ApplySymbologyFromLayer(), the shapefile is set to 4 classes, natural breaks, and the specified colors.
But when I then set the class break values and labels, the colors are reset to the default color scheme. I cannot use the ApplySymbologyFromLayer() afterwards, as then the classification scheme is reset.

EDIT: Obviously the problem are the manually set colors, as the colors remain consistent when a color scheme is used. But I need to apply the specific color values.

Why are the manually applied colors reset and is there a Workaround for this?

Here is my code:

import arcpy
inptShapefiles = arcpy.GetParameterAsText(0)
symbologyTemplate = arcpy.GetParameterAsText(1)

shapefiles = inptShapefiles.split(';')
for shp in shapefiles:
    name = os.path.splitext(os.path.basename(shp))[0]
    lyr = arcpy.MakeFeatureLayer_management(shp, name)
    arcpy.ApplySymbologyFromLayer_management(lyr, symbologyTemplate)
    lyr.getOutput(0).symbology.classBreakValues = [1, 100, 200, 300, 400]
    lyr.getOutput(0).symbology.classBreakLabels = ["Class A", "Class B", "Class C", "Class D"]
    arcpy.SaveToLayerFile_management(lyr, "E:/temp/" + lyr.getOutput(0).name, "RELATIVE")


texts:
There are several commercial and opensource tools for regression automation. While many of them could be employed to test GIS apps, they fall short, at times, dealing with the unique characteristics and use-cases of GIS apps and their imagery.

What are some tools for validating GIS apps?

I came across one such tool called Eggplant. Check out their video link in the bottom. They claim to not only compare data sets but also compare images. Seems interesting!

Are there other such commercial or opensource tools we can evaluate?

Thanks,
Ramesh


texts: Filtering polygons that intersect polyline dataset using ArcMap

texts:
I am trying to filter a polygon dataset that intersects with a polyline dataset.  I want all polygons that touch the polyline to be filtered out into a new shapefile.  I can't seem to make the intersect tool to work to do this as it won't allow polygon output with a polyline input.  I am working with ArcMap 10.6.
How do I make this dataset?


texts: How to work with .lyr files in ArcGIS Engine

texts:
I want to be able to load .lyr files created in ArcMap and ArcGIS Engine and save any updates in symbology in ArcGIS Engine back to the .lyr file.

What APIs should I use the only method I found so far is MapControl.addLayerFromFile() but the MapControl only let's me load a .lyr file?

I found an API during writing the question so I add an answer myself so it can be a wiki.


texts: Clipping raster with multiple polygons and naming the resulting rasters

texts:
I am trying to split a single-band raster (land cover) into smaller parts representing administrative units which are into another shapefile, containing their name as an attribute.

U I did this using SAGA "clip raster with polygon" (in QGIS 2.18) to split the raster and then naming the resulting rasters one by one (by hand) with the name of the administrative unit.

The problem is, I now need to do this process over ±2000 polygons, so doing it by hand is the last resort options.

I tried different options:
a) Using SAGA "clip raster with polygon", iterate over the polygon layer. Works good, but cannot specify an attribute to use a name.

b) Batch-process SAGA "clip raster with polygon" using individual administrative polygons, named after their administrative ID. Works fine, the name of the resulting rasters is as I want it, but the process is long and boring (batch-processing with 2000 rows, filled in part by hand, see image).

So, is there a way of speeding up the process? Is it possible to iterate over a layer, specifying an attribute to use as a name for the resulting layers? Is it possible to auto-fill the batch-processing window with layers open in the "layers panel"?


texts: Converting mercator meters without UTM zone?

texts:
I've been attempting to extract track information from Humminbird (Side Imaging Sonar) binary files and have hit a stumbling block with the coordinates it provides. The values are apparently in Mercator meters but no UTM zone is supplied in the output. I would like to know if anyone knows of a mercator meter conversion that doesn't require a UTM zone, and if so how do you do it? I've found one conversion which kind of does the job but not acurately enough (http://www.oziexplorer3.com/eng/eagle.html).

So basically I get the numbers  Lon: 16226043 Lat: -4557743. I know these should convert to 145.75528834210039 and -38.032398555522214 but can't work out how.

I've had a discussion about this issue on a more dedicated forum but came to a dead end (http://forums.sideimagingsoft.com/index.php?topic=5623.msg35417). I could be reading the binary section wrong but all other data in the line comes out fine.


texts: Performing MNF on multiple rasters in ENVI?

texts:
I am performing a forward MNF reduction in ENVI 5.6 on four different hyperspectral scenes (48 bands each). When I perform the MNF on each one, the Eigen values are different.
How should I go about handling this issue?
I considered mosaicking the images together, but apart from it not being best practice (according to ENVI help), you can only mosaic together images with up to 10 bands. Essentially, I am trying to produce a single coherent image with the least correlated bands.


texts:
Is there a way to use r.mapcalc to replace null with 0?

I tried if(isnull(map1), 0) but that produces a map that is entirely 0. I need only cells containing null to be altered!


texts: arcpy.SetParameterAsText fails in Python script tool

texts:
I am trying to create a derived output parameter for a Python script tool so that the output will show in a model (I want to use the output folder as an input for another tool). I want this output to be a folder that is the same as an input folder after geoprocessing has put the results (outShapeFiles) into that folder. Essentially, I have 4 initial parameters, one of which is the outFolder = sys.argv[4]. At the end of the code, I set `outFolder = arcpy.SetParameterAsText(5, "Output Features"). When I run it, everything works but the last part (the outShapefiles are created and put into the outFolder). I get this error message: "RuntimeError: Object: Error in setting parameter as text." Here is the code. Can you see how I can resolve this error so that my output parameter can be viable in a model? Note how I set outFolder at the end.

input_FC = sys.argv[1]
inField = sys.argv[2]
theFName = sys.argv[3]
outFolder = sys.argv[4] #Set as an input parameter in the script tool

script = sys.argv[0]
msg = "\nRunning: ... {}".format(script)
arcpy.AddMessage(msg)

arcpy.env.overwriteOutput = True
...

...
aMax = 1
for aVal in valueList:
  aMax = max(aMax, len(str(aVal)))
for aVal in valueList:
  if (str(aVal).isdigit()) and (not inType == "String"):
    fs = '"' + "%" + str(aMax) + "." + str(aMax) + 'i"'
    aSuffix = fs % aVal
    aVal = str(aVal)
  elif inType == "Double" and inScale == 0:
    aSuffix = str(aVal).replace(".0","")  ######
    aVal = str(aVal).replace(".0","")
  else:
    aSuffix = str(aVal)
    aVal = str(aVal)
  try:
    aSuffix = aSuffix.replace(" ","_")  #replace garbage in output files
    aSuffix = aSuffix.replace('"',"")
    aSuffix = aSuffix.replace("/","")
    aSuffix = aSuffix.replace("-","")
    outName = theFName + aSuffix + ".shp"
    outShapeFile = outFolder + "/" + outName
    outShapeFile = outShapeFile.replace("\\","/")
    #
    #Create a query and produce the file
    if (not aVal.isdigit()) or (inType == "String"):
      aVal = "'" + aVal + "'"
    whereClause = "%s = %s" % (inField, aVal)
    arcpy.MakeFeatureLayer_management(input_FC, "TempLayer", whereClause)
    arcpy.CopyFeatures_management("TempLayer",outShapeFile)
    arcpy.AddMessage("Output and query: " + outShapeFile + "  " +
    whereClause)
  except:
    whereClause = "%s = %s" % (inField, aVal)
    arcpy.AddMessage("Output and query: " + outShapeFile + "  " +
    whereClause + " did not work ")

arcpy.SetParameterAsText(5, outFolder) #This is where I try to
make outFolder an output parameter. I have set it to derived from the input "outFolder".

arcpy.AddMessage("\n  Processing complete" + "\n")
del arcpy


texts:
I try to create a interpolation with the follow script, the problem is that when the program tries to write QgsGridFileWriter, qgis crashes.
The problem is that the different examples are to lower version for qgis, and this script is to qgis3, someone knows if i have to write some parameter in output.writeFile(), like boolean of other versions?

#Generation Interpolate data#
 uri="file:///Z:/DataScience/TragsaMoscaOlivo2017/Tutoriales QGIS/Capa
 Predicciones/PREDICCION_oriental_I_S1.csvxField=LONG&yField=LAT&
 delimiter=%s"%(",")
 lyr=QgsVectorLayer(uri, 'Prediction','delimitedtext')
 values = (
 ('Low',0,0,'green'),
 ('Medium',1,1, 'yellow'),
 ('Large',2,2,'red'))
 QgsProject.instance().addMapLayer(lyr)
 # create a category for each item in values
 ranges = []
 for label,lower,upper,color in values:
   symbol = QgsSymbol.defaultSymbol(lyr.geometryType())
   symbol.setColor(QColor(color))
   rng = QgsRendererRange(lower, upper, symbol, label)
   ranges.append(rng)

# create the renderer and assign it to a layer
expression = 'pred_C' # field name
renderer = QgsGraduatedSymbolRenderer(expression, ranges)
lyr.setRenderer(renderer)
lyr.triggerRepaint()
#######################################################################
#Interpolation generation
layer_data=QgsInterpolator.LayerData()
layer_data.vectorLayer=lyr
layer_data.zCoordInterpolation=False
layer_data.interpolationAttribute=2 #Field
layer_data.mInputType=0 #points

tin_interpolator=QgsTinInterpolator([layer_data])
export_path ="Z:/DataScience/TragsaMoscaOlivo2017/Tutoriales
QGIS/pythoninterpolation/output4656.asc"
rect = layer2.extent()
res = 10
#ncol =int((rect.xMaximum() - rect.xMinimum())/res)
#nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)

#Grid generation and write
 output = QgsGridFileWriter(tin_interpolator,export_path,rect,300,300)
 print(output)
 output.writeFile()


texts: Assigning Elevation data from one point shapefile to another point shapefile in QGIS

texts:
I am working with QGIS. I have a point shapefile with all my elevation data and a point shapefile with the coordinates for my project. I would like to assign the elevations from the elevation data to my project file.

Is there a possibility to do that without creating a raster in-between?


texts: Execute a batch process of a spatial tool within an ArcGIS Toolbox

texts:
I have created a Data Interoperability Spatial ETL tool. I do have a Data Interoperability license.

I want to create a nightly batch process to run this tool.

As suggested, I have created a python script.

import arcpy,ConfigParser

# Read Config file
Config = ConfigParser.RawConfigParser()
Config.read('Config.txt')
TOOLBOXPATH=Config.get('ETL_Spatial','ToolboxPath')

class LicenseError(Exception):
    pass

try:
    if arcpy.CheckExtension("DataInteroperability") == "Available":
        arcpy.CheckOutExtension("DataInteroperability")
        print "Checked out \"DataInteroperability\" Extension"
    else:
        raise LicenseError
except LicenseError:
    print "Data Interoperability license is unavailable"
except:
    print arcpy.GetMessages(2)

# Import custom toolbox
arcpy.ImportToolbox(TOOLBOXPATH,'DataInterop')
arcpy.SpatialETLTool()

The value of TOOLBOXPATH is C:\Data Interoperability Tools\TransmissionDataRefresh\TransmissionDataRefresh.tbx.

But again I am getting below mentioned error while executing the script:-

Checked out "DataInteroperability" Extension
  Traceback (most recent call last):
  File "C:\Data Interoperability Tools\TransmissionDataRefresh.py", line 24, in <module>
    arcpy.ImportToolbox(TOOLBOXPATH,'DataInterop')
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\__init__.py", line 125, in ImportToolbox
    return import_toolbox(input_file, module_name)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\toolbox_code.py", line 434, in import_toolbox
    toolbox = gp.createObject("Toolbox", tbxfile)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\_base.py", line 381, in createObject
    self._gp.CreateObject(*gp_fixargs(args, True)))
IOError: The toolbox file C:\Data Interoperability Tools\TransmissionDataRefresh\TransmissionDataRefresh.tbx was not found.

Please help


texts:
I have a problem with a fixed layer.

I have to count accident by Link within 100m area.

When I use a Fixed Distance Buffer, there is a problem because some area is overlayed so they will count 1 accident to 2.

So I used the Node Tool.

I want to know is there any automatic way for this solution?

Key: All of the layers must not be overlayed.


texts: Using Trace geometric network in ArcPy script runs out of memory

texts:
I have written a short script to trace a geometric network representing rivers and creeks upstream from a number (10000~) of road intersections that are potential blockers to migrating fish. What we want to do is to evaluate where restorative measures would be most effective in terms of opening up longer stretches of water. My problem is that the script runs out of memory and crashes after 1500~ iterations.
The script basically selects the first junction in the network, traces upstream from there, adds the length of the resulting selected lines through another for loop, sums these up and updates the length field in the cursor, before it iterates to the next junction. Never mind the other fields in the row, in the end the script is supposed to calculate a bunch of other things for each junction too.
I haven´t had this issue before, even with loops of similar style involving far more iterations, which makes me wonder if there´s something particularly ineffective with how it´s written. The script looks like this:
import arcpy as ap

ap.env.workspace = r"X:/8853_Vandringshinder_trafikverket/ArbetsGIS/arbetsfiler/Geometric_final.gdb/"

punkter =  r"X:/8853_Vandringshinder_trafikverket/ArbetsGIS/arbetsfiler/Geometric_final.gdb/Geometric_final/Vag_vandringshinder_final"
linjer = r"X:/8853_Vandringshinder_trafikverket/ArbetsGIS/arbetsfiler/Geometric_final.gdb/Geometric_final/Geometric_final_Net"
fields = ('OBJECTID', 'E_Langd_NED', 'E_Langd_UPP', 'E_Tot_hinder', 'E_tot_artfynd')
arter = r"X:/8853_Vandringshinder_trafikverket/Inkommet_material/20201204_Artobservationer_Analysportalen/SpeciesObservations-2020-12-04-10-49-02/Observations.shp"

cursor = ap.da.UpdateCursor(punkter, fields)
count = 1

for row in cursor:
    query1 = "OBJECTID ="+ str(row[0])
    ap.SelectLayerByAttribute_management(punkter, "NEW_SELECTION", query1)

    linjer_trace = ap.TraceGeometricNetwork_management(linjer, "linjer_trace", punkter, "TRACE_UPSTREAM")

    Tot_langd_upp = []

    rows = arcpy.SearchCursor(linjer_trace)
    for rowdy in rows:
        langd = rowdy.getValue("Length")
        Tot_langd_upp.append(langd)

    row[2] = sum(Tot_langd_upp)
    cursor.updateRow(row)

    print("Row "+str(count)+" out of 10027")
    count+=1

Any thoughts on this?
I'm fairly new to python, and only do kind of basic this like this. I´ve also noticed I need to be in edit mode to run the script, which I don´t think I´ve experienced before. And - the script is run  in the python console in ArcMap 10.8.3


texts: Rotation field is missing in Layer Properties > Symbology > Advanced in QGIS 3.10. How to fix this?

texts: In the layer properties under symbology- advanced properties, the Rotation filed is missing.


texts:
I have road network (lines) and points (Address) . With the help of closestSegmentWithContext function, I join points to nearest lines in QGIS. Is there a way to to generate a nodes at every joining point as I have drawn manually in below figure for demonstration?

The code for joining line segements to nearest points is given below and shape files are available here:
roadlyr = QgsProject.instance().mapLayersByName('TR_ROAD')[0]
addrlyr = QgsProject.instance().mapLayersByName('ADDRESS')[0]

#List all features
adresses = [f for f in addrlyr.getFeatures()]
roads = [f for f in roadlyr.getFeatures()]

#Create empty vector layer
vl = QgsVectorLayer("LineString?crs={}&index=yes".format(roadlyr.crs().authid()), "myLayer", "memory")
provider = vl.dataProvider()

for adress in adresses:
    a = min([road.geometry().closestSegmentWithContext(adress.geometry().asPoint()) for road in roads], key=lambda x: x[0]) #For each road find closest address
    gLine = QgsGeometry.fromPolyline([QgsPoint(adress.geometry().asPoint()), QgsPoint(a[1])])
    f = QgsFeature()
    f.setGeometry(gLine)
    provider.addFeature(f)

QgsProject.instance().addMapLayer(vl)


texts:
I am using OpenLayers map by using the hosted JavaScript:

<script src="http://openlayers.org/api/2.13.1/OpenLayers.js"></script>

But my client has SSL installed and when I try to run my map page it shows:

(2)[blocked] The page at https://domain.com/rwd/ ran insecure content
   from http://openlayers.org/api/2.13.1/OpenLayers.js.

So I tried https one and it turns out openlayers doesn't have one

https://openlayers.org/api/2.13.1/OpenLayers.js

Then I download the Openstreet js and hosted in client server, but then all the styles and related images are lost. Although it does show the map, basically numerous warnings pop-up in the console and I'm afraid this might get rejected at play store or so. I'm developing a hybrid application which runs on server as well.

Warnings now :

The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://b.tile.openstreetmap.org/14/8743/5624.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://b.tile.openstreetmap.org/14/8742/5624.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8743/5623.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://b.tile.openstreetmap.org/14/8743/5625.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8744/5624.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://a.tile.openstreetmap.org/14/8742/5623.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8742/5625.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://a.tile.openstreetmap.org/14/8744/5623.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://a.tile.openstreetmap.org/14/8744/5625.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://a.tile.openstreetmap.org/14/8743/5622.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8742/5622.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8745/5624.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://b.tile.openstreetmap.org/14/8744/5622.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://b.tile.openstreetmap.org/14/8745/5623.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8745/5625.png.
The page at https://domain.com/rwd/#/customer-plot/234 displayed insecure content from http://c.tile.openstreetmap.org/14/8745/5622.png.

GET https://domain.com/rwd/js/lib/theme/default/style.css 404 (Not Found)
/*This the corresponding stylesheet that is loaded via Openstreet.js*/

I tried with cdn with SSL and still same issue:

https://cdnjs.cloudflare.com/ajax/libs/openlayers/2.11/OpenLayers.js

I downloaded the entire GitHub repo for Openstreet thinking this will solve the image and CSS issue.

https://github.com/openlayers/openlayers

While the above solved the CSS issues, the map however is shown by loading out images from another external site tile.openstreet.com. Seems like I have to dig into openstreet js as well..:(..

Important : I am using backbone.js.

Note : I have gone through OpenLayers 2.12 and http basic authentication woes and that has not helped me. I do not have any control over the server configuration. I just have access to a folder where I need to set up the website an everything else works fine, but this SSL is troublesome.

Changed real website address to avoid Google linking it.


texts: Why does QGIS Profile plugin lack profile graph?

texts:
I've searched but can't seem to find a similar issue with this plugin - the plugin starts ok and seems to run but I don't get any view of the profile or some other buttons in the UI as I've seen elsewhere - please see image. I do get the the table of heights which I can export into Excel for graphing but I want to see the profile in QGIS.
I'm running Dufour on a Windows 8 machine and all install using OSGeo4W.


texts: I am having a problem in creating new signatures in Signature Editor in ERDAS 2014. Whenever I click on the icon (i.e., after selecting the aoi), it shows an error which says "unable to open associated image: Image association may be edited with Edit menu (Image Association...)". After I had set the associated image from the edit menu of the Signature Editor and tried to create new signature, it displayed "Viewer image does not match associated image for signature set. Image association may be edited with Edit menu (Image Association...)". PLEASE HELP!!!


texts:
This is seemingly pretty simple.  I have a csv file I want to add to the map. It is comma delimited and there is no lat/lon information. I will join it to a shapefile later. I used other examples from stackoverflow and this documentation in the QGIS cookbook to do it. I can add it with the following code:

uri = 'file:///Users/ep9k/Desktop/SandraMonson/CountByZip.csv?delimiter=,'
csv_file = QgsVectorLayer(uri, 'Patient Data', 'delimitedtext')
QgsProject.instance().addMapLayer(csv_file)

This adds the file to the map. But in the attribute table it reads all data as column headers. The attribute table looks like this:

ZipCode, PatientCount, StateOrProvince should be the column headers, with data below. How can I read the first line as column headers and the rest as data in those columns?

Another thought I have, after looking at the QGIS cookbook documentation, I went into PyQGIS API documentation and into Module>Core Library> and looked for QgsVectorLayer and QgsProject, and it seems counterintuitive as a csv file is not a vector layer.


texts: Do Navigation Apps Shrink the Hippocampus?

texts:
A study shows London cab drivers' hippocampi grow compared to a control group of similar non-drivers.  I presume these cab drivers don't rely on navigation apps.

  ...acquisition of an internal spatial representation of London was
  associated with a selective increase in gray matter (GM) volume in
  their posterior hippocampi and concomitant changes to their memory
  profile.

Are there any studies examining the effect of Navigation app usage on hippocampus size?

I searched, but could not find any similar studies.


texts: Error in OpenLayers: c.de is not a function

texts:
I've used vector tile in OpenLayers 4.1 to show vectors on top of OSM base map. The problem is, it shows a white screen and throws the "c.de is not a function" in the console. I debugged the code and found out it's related to vectorTileSource part but I can't get any detail. What's the problem with my code? The script part of the code attached here:

<script>

var baseLayer = new ol.layer.Tile({
    source: new ol.source.OSM()
});

var fill = new ol.style.Fill({
    color: 'rgba(0,0,0,2)'
});
var stroke = new ol.style.Stroke({
    color: 'rgba(0,0,0,4)'
});
var circle = new ol.style.Circle({
    radius: 6,
    fill: fill,
    stroke: stroke
});
var vectorStyle = new ol.style.Style({
    fill: fill,
    stroke: stroke,
    image: circle
});
var vectorTileSource = new ol.source.VectorTile({
    format: new ol.format.TopoJSON({
        defaultProjection: 'EPSG:4326'
    }),
    projection: 'EPSG:3857',
    tileGrid: new ol.tilegrid.createXYZ({
        maxZoom: 19
    }),
    url: 'http://{a-c}.tile.openstreetmap.us/vectiles-water-areas/{z}/{x}/{y}.topojson'
});
var vectorTileLayer = new ol.layer.Vector({
    source: vectorTileSource,
    style: vectorStyle
});

var lookAt = ol.proj.transform([51.34, 35.65], 'EPSG:4326', 'EPSG:3857');
var view = new ol.View({
    center: lookAt,
    zoom: 12
});

var map = new ol.Map({
    target: 'map',
    layers: [baseLayer, vectorTileLayer],
    view: view
});


texts:
The purple is the original polygon to be smoothed and I would like to smooth the polygons as the red line using Python.

I tried simplify but looks like I need a buffer and simplify the buffer. But I want to simplify the polygon directly. I am open to any other methods using Python.

import geopandas as gpd

d = gpd.read_file('poly.shp')
x = d.buffer(0.0001)

tolerance = 0.00005
simplified = x.simplify(tolerance, preserve_topology=True)
simplified.to_file('simplified.shp')


texts: generic "There is a problem with your connection" error in CartoDB

texts:
I changed the cartoCSS for a visualization, but in a pretty straightforward way (fyi: changed some numerical cutoffs to other numbers, did not touch anything else). I am getting a "There is a problem with your connection" error on the visualization ever since, for many hours now, in multiple browsers.

Please let me know if there is anything to do about this apart from some deep plumbing only tech support can do. CartoDB staff is very helpful, but the turnaround because of the time zone difference with Madrid is getting frustrating.

Here is the new piece of cartoCSS:

#location_of_winners25mb [ hufvalue <= 15530191054167961601] {
marker-width: 25.0;
marker-fill: #850200;
}
#location_of_winners25mb [ hufvalue <= 226926540] {
marker-width: 23.3;
marker-fill: #B10026;
}
#location_of_winners25mb [ hufvalue <= 87005948] {
marker-width: 21.7;
marker-fill: #F11810;
}
#location_of_winners25mb [ hufvalue <= 46218709] {
marker-width: 20.0;
marker-fill: #FC4E2A;
}
#location_of_winners25mb [ hufvalue <= 27412874] {
marker-width: 18.3;
marker-fill: #FD8D3C;
}
#location_of_winners25mb [ hufvalue <= 18665512] {
marker-width: 16.7;
marker-fill: #FF9900;
}
#location_of_winners25mb [ hufvalue <= 12041280] {
marker-width: 15.0;
marker-fill: #FFA300;
}
#location_of_winners25mb [ hufvalue <= 7283027] {
marker-width: 13.3;
marker-fill: #FFCC00;
}
#location_of_winners25mb [ hufvalue <=  3711385] {
marker-width: 11.7;
marker-fill: #fed976;
}
#location_of_winners25mb [ hufvalue <= 1316117] {
marker-width: 10.0;
marker-fill: #FFFFB2;
}

And here is the old one:

#location_of_winners25mb [ hufvalue <= 88302118874.01] {
   marker-fill: #d73027;
}
#location_of_winners25mb [ hufvalue <= 83885841] {
   marker-fill: #f79272;
}
#location_of_winners25mb [ hufvalue <= 32596503] {
   marker-fill: #fed6b0;
}
#location_of_winners25mb [ hufvalue <= 17639686] {
   marker-fill: #fff2cc;
}
#location_of_winners25mb [ hufvalue <= 9809490] {
   marker-fill: #d2ecb4;
}
#location_of_winners25mb [ hufvalue <= 4741029] {
   marker-fill: #8cce8a;
}
#location_of_winners25mb [ hufvalue <= 1592812] {
   marker-fill: #1a9850;
}


texts: Determining if route self intersects using Java?

texts:
I have a sequence of coordinates (lat,lng) defining a route.

How can I determine whether this route is self intersecting?

I am trying to implement a Java method for this, but I have been unsuccessful so far.


texts:
Is it possible to use PostgreSQL with QGIS and ArcGIS? i.e. one database for different clients.

Is it anything that I have worry or any problems to use the DBMS from two different software?

I know that I need PostGIS and ArcSDE and I can't manipulate or store raster data with PostGIS 1.5 in PostgreSQL from QGIS.

What would you recommend or what to avoid?


texts: How to geocode addresses that are outside Address Ranges of source layer?

texts: Im using TIGER lines to geocode a table of addresses. Say for instance the address range for a street is 1100-1200. What would be the best way to match an address such as 1094?


texts: QGIS calculate the 3D surface area of a region?

texts:
I would like to calculate the true surface area of a region from DEM (2D elevation data).

Is it possible to do this with QGIS?

And if yes, how?

I managed to obtain a 3D visualization of the elevation data using nviz (a GRASS-module) but I do not see any possibility to calculate the 3D surface area.


texts: Open source desktop GIS for high school curriculum?

texts:
I am a physics teacher and contemplating to integrate GIS to our high school's curriculum. A de facto software looks to be ArcGIS, which is a proprietary software.
I wonder what are  open-source alternatives that are suitable for high school students.

Are there any open-source desktop GIS and corresponding books that can be integrated into a high school curriculum?


texts: Merging multiple 16-bit image bands to create a True Color Tiff

texts: I'm trying to merge three Tiff band images in python, red, green and blue into a True color image. PIL does not support the file depth as they are 16-bit files with 11-bit values. I have seen it mentioned to use a sequence of Image objects. However after attempting to research this I'm unsure how to approach it.


texts: gdal_translate errors with VRT: GetBlockRef failed and IReadBlock failed

texts:
I am projecting global rasters with the Homolosine. gdalwarp reprojects it nicely, but it does not compress the output, that must be done afterwards with gdal_translate. To reduce computation time, I tried applying gdalwarp into a VRT and then the gdal_translate into GeoTIFF.

The warping carries out fine:

$ gdalwarp -overwrite -s_srs EPSG:4326 -t_srs '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' -multi -co "BIGTIFF=YES" -te -20037500 -6729000 20037250 8600750 --config CHECK_WITH_INVERT_PROJ TRUE -tr 250 250 -of VRT input/TEST.tif output/TEST.tif.vrt
Creating output file that is 160299P x 61319L.
Processing input file input/TEST.tif.
Using internal nodata values (e.g. 65535) for image input/TEST.tif.
Copying nodata values from source input/TEST.tif to destination output/TEST.tif.vrt.

Checking with gdalinfo, everything looks in place:

$ gdalinfo output/TEST.tif.vrt
Driver: VRT/Virtual Raster
Files: output/TEST.tif.vrt
       input/TEST.tif
Size is 160299, 61319
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Interrupted_Goode_Homolosine"],
    UNIT["Meter",1]]
Origin = (-20037500.000000000000000,8600750.000000000000000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  SCALE=0.0001
Corner Coordinates:
Upper Left  (-20037500.000, 8600750.000)
Lower Left  (-20037500.000,-6729000.000)
Upper Right (20037250.000, 8600750.000)
Lower Right (20037250.000,-6729000.000)
Center      (    -125.000,  935875.000) (  0d19'37.24"W,  8d24'25.59"N)
Band 1 Block=512x128 Type=UInt16, ColorInterp=Gray
  NoData Value=65535

The final step with gdal_translate:

$ gdal_translate -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "BIGTIFF=YES" -of GTiff output/TEST.tif.vrt output/TEST.tif
Input file size is 160299, 61319
0ERROR 1: Too many points (441 out of 441) failed to transform, unable to compute output bounds.
ERROR 1: output/TEST.tif.vrt, band 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0

Is this behaviour expected? Could it be a bug?


texts: Generate Tessellations looping through all shapefiles from folder in ArcGIS Pro ModelBuilder?

texts:
I am trying to create Hexagonal grids for each shapefile in a folder. I opted for ArcGIS Pro ModelBuilder tool and was successful in creating a loop that reads and loads shapefile in the memory. But afterwards, I am trying to grab the extents of each shapefile and want to generate Hexagonal grids for each polygon. Struggled a lot and still have not figured out where the issue is.

Attached copy of model that I built for reference.


texts: Different prime meridian in equal earth projection with QGIS?

texts:
Is it possible to use the equal earth projection with a prime meridian other than Greenwich (0°) in QGIS?

Adding Equal Earth Projection in QGIS? shows how to use the equal earth projection with QGIS, but the prime meridian is "fixed" at 0° (Greenwich)


texts: How to prop up/convert a 2D line to 3D using DEM in PostGIS?

texts:
I have planar (2D) polylines, which I'd like to convert to 3D ones using a DEM of the same region. I am wondering if this can be done in PostGIS?

I guess I can program a SQL function that reads all the vertices of a polyline/LineString, and then promotes each vertex to 3D one by one. But I just wanted to ask first to avoid reinventing wheels.

A related question (Joining 2D polylines and DEM into 3D polylines using QGIS/GRASS/SAGA?) asked about how to do this in QGIS but I need to work in PostGIS.


texts:
To avoid carrying books and optimize things while traveling, I like to build a Google Map with sights and cafés/restaurants. I usually print it on paper or keep a PDF in my phone so that I don't pay for a data connection.

But single Google Map is unable to take a bunch of addresses as input and doesn't offer a way to number locations, it's very time-consuming. I had to find and download numbered thumbtacks. Here's an example of Prague I started today.

Are there add-on's to either Google Maps or OpenStreetMap that provide a better solution? Ideally, I could just paste a list of addresses, it would add all the locations on a map, each with a number, and let me re-order the locations easily.

Thank you.

Edit: Google Maps has the option of importing date from a KML, KMZ, or GeoRSS file. Could this be the solution?


texts: Trimming service area polygons in ArcGIS Pro?

texts:
Does anyone know how to trim polygons for a service area network analysis in ArcGIS Pro?

It seems to have a default of 100 meters and I would like to not trim them.

I can't use ArcMap, which easily allows for this. Unfortunately it keeps crashing with the complexity of the network I'm trying to run.


texts: What data source is behind the ArcGIS GeoEnrichment service?

texts: I would like to know if there is a reference to the data of the GeoEnrichment service, available with a ArcGIS Online account. I only found this post of 2013 where it seems there is no reference. Does anyone know anything about it? I would like to know which data I am using. In particular, the data I would like to enrich are from Milan (Italy).


texts: Conversion to Leaflet from Google Maps Throws Image Overlay Out of Alignment?

texts:
I'm in the middle of transitioning away from Google Maps and to Leaflet due to the recent change in Google Maps usage policy.  I've converted all but one of my map pages, but unfortunately, that one is the most popular.

I must be missing something obvious, because the image overlay displays correctly on Google Maps, but taking the exact same bounds/images and putting them into Leaflet causes a shift.  I've read in other answers that such subtle shifts are caused by coordinate conversion problems, but everything I'm reading says that Google Maps and Leaflet use the same coordinate system.  What am I missing?

Here's an example of the Google Maps page:  https://www.rabbitears.info/contour.php?appid=8a441da96f4148f5b2c9cf6212b6c0ea&site=1&map=Y

And here's an example of the Leaflet page:  https://www.rabbitears.info/contour_leaflet.php?appid=8a441da96f4148f5b2c9cf6212b6c0ea&site=1&map=Y

You can see the overlay has shifted slightly northward, enough so that the image (television station coverage area) no longer lines up with the terrain.


texts:
I need to automaticly overlay/snap the red line (gps-track) exact on the lines of my roads-shape.
Is that possible?


texts: Setting proper rendering order because offset causes odd corners for lines

texts:
I used offset on these lines that are on top of each other. How can I make it so the red line is not interrupted by the orange?
I can put the orange under the red, but there is another corner where the opposite then would occur.

Each line segment is a new attribute. These lines are on top of each other but used offset on one of them on this section. How can I prioritise them without it affecting the whole line?


texts: Access ArcGIS Server directories in geoprocessing service

texts:
I would like to save a file in the output folder on the ArcGIS Server and then provide a link for the user so he can download it.

I have two problems:

I can't figure out how to access the output folder in arcpy - because the geoprocess has to run locally first and only then it is uploaded to the server - how do i tell ArcGIS to use a folder on the server on my local machine, it doesn't make any sense to me.

I understand that it is possible to download a file via rest once it is the output folder - but what exactly is the URL I'm supposed to give the client? The ArcGIS Server seems to change the name of the directory I create.


texts: QGIS .adf file saved as GTiff changes elevations

texts:
I have got two .adf DEMs which I need in GTiff format so they can be merged. When I 'Save As' the layer as a GTiff, the resulting output has elevations which are inconsistent with the original .adf file. For example: the max elevation in the first .adf file is 2943.76 m and the same file converted to GTiff format has a max elevation of 2881.86 m.
Where is this discrepancy coming from and how can I ensure this doesn't happen?


texts:
I am attempting to work through "Getting to know Web GIS" by Pinu Fu, and I think it may be out of date with the current arcgis online. I am attempting to add an image from a URL (http://esrimapbook.esri.com/GTKwebgis/chapter2/blue_down.gif) and the channge symbol window (change style>counts and amounts (color), Options>Symbols) Yet unfortunatly it does not have an add an image option (see below):

am I looking in the wrong place?

No blogs I have found seem to help, how do you add an image to the symbols?


texts:
I am currently developing an inspection map using AGOL/ArcMap and the Collector app. I have a series of potential deficiency sites i.e. upgradient, downgradient, lateral, etc... all possible at the same area. In my map I have domains setup for the user to pick a site location (as described above)… is there a way for me to link the domain with a symbol/color? i.e. the user picks "upgradient" and the symbol automatically turns to a red circle, or if the user picks "downgradient" it turns to a green circle.

As of now the only way I can figure to do this is after exporting the data into ArcMap and adjusting the symbology on my desktop.


texts: Configuring mod_tile and renderd for memcache

texts:
I have my mod_tile.config set up like this

# This is the Apache server configuration file for providing OSM tile support
# through mod_tile

LoadModule tile_module modules/mod_tile.so

<VirtualHost *:80>
    ServerName tile.openstreetmap.org
    ServerAlias a.tile.openstreetmap.org b.tile.openstreetmap.org c.tile.openstreetmap.org d.tile.openstreetmap.org
    DocumentRoot /var/www/html

# Specify the default base storage path for where tiles live. A number of different storage backends
# are available, that can be used for storing tiles.  Currently these are a file based storage, a memcached
# based storage and a RADOS based storage.
# The file based storage uses a simple file path as its storage path ( /path/to/tiledir )
# The RADOS based storage takes a location to the rados config file and a pool name ( rados://poolname/path/to/ceph.conf )
# The memcached based storage currently has no configuration options and always connects to memcached on localhost ( memcached:// )
#
# The storage path can be overwritten on a style by style basis from the style TileConfigFile
    ModTileTileDir memcached://
...

and my renderd.config set up like this

[renderd]
num_threads=4
tile_dir=memcached://localhost:11211
stats_file=/var/run/renderd/renderd.stats

[mapnik]
plugins_dir=/usr/lib/mapnik/3.0/input
font_dir=/usr/share/fonts/truetype
font_dir_recurse=1

[ajt]
URI=/hot/
TILEDIR=memcached://localhost:11211
XML=/home/renderaccount/src/openstreetmap-carto/mapnik.xml
HOST=localhost
TILESIZE=256
MAXZOOM=20

yet I am receiving this error:

renderd[10095]: Loading parameterization function for debug: init_storage_backend: initialising memcached storage backend at: memcached://
ERROR: init_storage_memcached: Support for memcached has not been compiled into this program

How do I compile support for memcached into mod_tile and renderd?


texts: Mapping a NetCDF in QGIS in the correct location on Earth

texts:
I was wondering if anyone has encountered NetCDFs being mapped in the "wrong" location by QGIS? I'm using QGIS 3.16 and I've tried 2.18 too. I think the problem has to do with the longitude units of the NetCDF being in degrees east and QGIS not understanding how to deal with them properly. When I add other layers (vector or raster) they map correctly to my area of interest (British Columbia). My map display is set to WGS 84 (EPSG: 4326) and I import my file.nc as a mesh layer (I've tried raster too but it doesn’t make a difference) and set the CRS of the mesh layer to WGS 84 (4326) too after it has been imported. The NetCDF file imports and visualizes without any problems, it is just in the wrong place (i.e. it doesn’t appear anywhere near my study area. When I zoom into the mesh layer, The Coordinate box at the bottom of the display gives the longitude (degrees east, approx. 225), latitude (approx 50), and then when I zoom to a vector layer in my study area the coordinate box shows the longitude (approx.-128) and latitude (approx. 50).
I don't think the problem lies with the NetCDF file as longitude degrees east are common, I think it is how QGIS is mapping it that may be the problem, like it doesn't know how to convert the longitude degrees east.
The NetCDF has a single variable (Chlorophyll), with 4 dimensions (time, altitude, longitude and latitude). I downloaded it from https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMBchla1day.html using this bounding box: lat degrees north: 48, 55, lon degrees east: 227, 234
Also, I have tried mapping it in ArcGIS Pro, using the NetCDF to raster tool, and it works perfectly. Unfortunately, I need to be able to use it in QGIS as my ArcGIS Pro trial is about to expire.
I thought perhaps subtracting the longitude from 360 would solve the problem, so I used xarray to do so. QGIS crashed on trying to import this modified NetCDF so I assume that was not a good solution.
Any advice?


texts: How to remove redundant points using Open Source GIS (preferably)?

texts:
We have points table including vehicle routes. In huge traffic conditions, devices collecting data on route may output redundant points where at any coordinate same points may repeat. We can easily search and remove those points, but we have another issue like this.

Since points have threshold like 5 meters (for example) among them, when in traffic output may occur so many close (adjacent) points that we cannot seem to analyse. Also we can't adjust device, we have to move on with what we have.

We need to implement an algorithm or script to remove these redundant points. Or at least we can join them into one as we calculate some fields using functions like avg, sum.

How can we achieve this?

We have coordinates in MapInfo table but as points we can move them into anything, your answers do not have to be in one platform.

Open source solutions would be ideal.


texts: ARCMAP Symbology: Category by Expression?

texts:
I'm trying to categorize symbology by a particular portion of a column (the left 3 characters).  This layer is based on a SQL view that I don't have access to modify -- and I don't want to create a copy because it needs to be updated daily.  (So I suspect creating a calculated field isn't possible)

I want it to look something like this:
ID     Name      Category
1      ABC123    ABC
2      ABC828    ABC
3      DEF381    DEF

Then color by the category column (which doesn't exist currently).  I know I could manually change the symbol and description on each item, but there are hundreds and more show up randomly.  The first 3 letters should be consistent though.

Any way to accomplish this without database access, while retaining the ability for the layer to update dynamically as the data changes in the source?


texts: Is ArcGIS Online suitable for storing rasters which need further processing?

texts:
I have a geo-referenced RGB aerial image, which I've added to a new MXD. Now I wish to share it with other people, so I publish the MXD to ArcGIS Online.

Unless I'm mistaken, the only option for publishing the raster is to store it as a Tiled layer.

Is it possible to perform further analysis on this tiled layer on ArcGIS Online (for example, to calculate the NDVI) or is the tiled layer on ArcGIS Online only suitable for visualisation?


texts: Can OVERLAPBDYDISJOINT or ON ever occur between two polygons?

texts:
Oracle uses a series of specific spatial relationships to indicate the how different shapes interact with each other.

Two of the relationships it describes are somewhat odd:

  OVERLAPBDYDISJOINT: The interior of one object intersects the boundary and interior of the other object, but the two boundaries do not intersect. This relationship occurs, for example, when a line originates outside a polygon and ends inside that polygon.

and

  ON: The interior and boundary of one object is on the boundary of the other object. This relationship occurs, for example, when a line is on the boundary of a polygon.

Am I correct that these can never happen with two shapes that are the same inherent dimensionality? That is, a polygon can never have this relationship with another polygon, and a line can never have this relationship with another line. It can only happen if I have a point and a polygon, a line and a polygon, or a point and a line. Or have I misunderstood?


texts: How put Leaflet Filelayer data into editable layer so I can edit after Upload

texts:
I am using “togeojson” and “Leaflet.FileLayer” for File upload in my leaflet draw using below Code:

    L.Control.FileLayerLoad.LABEL = '<img class="icon" src="folder.svg" alt="file icon"/>';
    var style = {color:'red', opacity: 1.0, fillOpacity: 1.0, weight: 2, clickable: false};
  var LoadFile =  L.Control.fileLayerLoad({
     position: 'topright',
    layer: L.geoJson,
    // See http://leafletjs.com/reference.html#geojson-options
    layerOptions: {style: {color:'red'}},
    // Add to map after loading (default: true) ?
    addToMap: true,
    // File size limit in kb (default: 1024) ?
    fileSizeLimit: 1024,
    // Restrict accepted file formats (default: .geojson, .kml, and .gpx) ?
    formats: [
        '.geojson',
        '.kml',
        '.gpx'
    ]
}).addTo(map);

How do I access upload file data so I can put it into editable layer and edit uploaded draw and save as geojson?

I have previously load geojson to leaflet draw and edite json + save to database but not success to get upload from directory.


texts: Is there a Python expression for Does Not Contain that I can use in a Update Cursor? I am looking for the '@' symbol in a column within a table (not a FC) and want to delete all entries without one, is there a way accomplish this?


texts: Problem 'unable to open datasource" with ogr2ogr to reproject shapefile

texts:
I have no idea why I get the following error:

FAILURE:
Unable to open datasource `Sydney2.shp' with the following drivers.

I have already checked successfully the shape with ogrinfo.

I am using the following command:

ogr2ogr -f "ESRI Shapefile" Sydney.shp Sydney2.shp -s_srs EPSG:32656 -t_srs EPSG:4326

Any ideas?


texts: Removing a Lock From a Field in ArcGIS Online

texts:
I'm trying to remove the lock on a map layer in ArcGIS Online that belongs to a map being used in field collection. It's missing some data. I've tried un-indexing it and that did nothing.
Does anyone have any suggestions?
Picture below added for clarification.

I've tried configuring popup in map viewer to allow for the field to be edited, but it's grayed out and I'm unable to do so.


texts: Error exporting contour lines from Google Earth Engine

texts:
I would like to export these contour lines as .SHP but GEE says that variable collection must be a FeatureCollection.

Here is the complete code


texts: Access Elements in SVG with Openlayers 3 using forEachFeatureAtPixel()

texts:
I am trying to access elements (id, style) of a SVG image in Openlayers 3. I tried to link different examples, but the problem is that displaying a SVG is only possible with ImageStatic() and the SVG properties get lost. I am using forEachFeatureAtPixel() to react to different mouse events, which worked fine for a KML layer. How can I access the SVG content with forEachFeatureAtPixel()? A minimum working example is given below:

 <!DOCTYPE html>

    <title>Static Image</title>
    <link rel="stylesheet" href="http://openlayers.org/en/v3.16.0/css/ol.css" type="text/css">


texts: Automating the updating of the attribute table of a shapefile that is linked to a CSV file

texts:
I know this has been discussed before, and I used the code provided by @Joseph in How to join data from excel to an attribute table in QGIS without creating duplicates for my Macros.

I am trying to achieve what @Joseph said, but I am attaching a screenshot of the error that I got when I did what you said. Below is my complete macros code.

'Noviembre_P01' is my shapefile
'project_2_gis_mapping.csv' is the csv file that I want to be able to connect automatically to this shapefile so that whenever my csv is updated, the attribute table in the shapefile also updates itself.
'LOC_CODE' is the common attribute to my .shp and .csv files.
'first_bloo' and 'diagnosis' are the 2 fields I am interested in mainly, out of all the fields in the csv file.

I also did not create a join like he mentioned.

Conversely, will using 'Watch File' in the 'Add Delimited Text Layer' dialogue box server the exact same purpose?

def openProject(): from PyQt4.QtCore import *
import qgis
from qgis.core import QgsMapLayerRegistry, QgsVectorJoinInfo

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    # Change to your shapefile name
    if layer.name() == "Noviembre_P01":
        qgis.utils.iface.setActiveLayer(layer)
        shp = qgis.utils.iface.activeLayer()

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    # Change to your csv name
    if layer.name() == "project_2_gis_mapping.csv":
        qgis.utils.iface.setActiveLayer(layer)
        csv = qgis.utils.iface.activeLayer()

# Set up join parameters
shpField='LOC_CODE'
csvField='LOC_CODE'
joinObject = QgsVectorJoinInfo()
joinObject.joinLayerId = csv.id()
joinObject.joinFieldName = csvField
joinObject.targetFieldName = shpField
shp.addJoin(joinObject)

# Define fields to update and joined fields to copy values from
ip1 = shp.fieldNameIndex('first_bloo')
ip1_join = shp.fieldNameIndex('project_2_gis_mapping.csv_first_bloo')
ip2 = shp.fieldNameIndex('diagnosis')
ip2_join = shp.fieldNameIndex('project_2_gis_mapping.csv_diagnosis')

shp.startEditing()
for feat in shp.getFeatures():
    shp.changeAttributeValue(feat.id(), first_bloo, feat.attributes()[first_bloo_join])
    shp.changeAttributeValue(feat.id(), diagnosis, feat.attributes()[diagnosis_join])
shp.commitChanges()

# Remove join
shp.removeJoin(csv.LOC_CODE())
    pass


texts: Using a ShapeFile spatial index with ogr2ogr and ogr.DataSource.ExecuteSQL()

texts:
I created a spatial index foo.qix for a foo.shp with ogrinfo -sql "CREATE SPATIAL INDEX ON foo"  foo.shp However, when I do a query with ogr2ogr and look at it in strace, it doesn't use that index:

strace -e open,access  /opt/gdal-custom/bin/ogr2ogr -f CSV /vsistdout/  foo.shp -dialect sqlite -sql "SELECT * FROM foo WHERE ST_Intersects(GEOMETRY,BuildCircleMbr(0, 0, 150000, 4326))" 2>&1 | grep foo
open(" foo.shp", O_RDONLY) = 3
open(" foo.shx", O_RDONLY) = 4
open(" foo.dbf", O_RDONLY) = 4
open(" foo.cpg", O_RDONLY) = -1 ENOENT (No such file or directory)
open(" foo.CPG", O_RDONLY) = -1 ENOENT (No such file or directory)
open(" foo.prj", O_RDONLY) = 5

Once I get this working, if there is a speed improvement, I'd also like to apply spatial index to a similar query from python via ExecuteSQL on the shapefile.


texts: Apply Grid Shift Transform file to GeoServer

texts:
Trying to apply a grid shift file to GeoServer.  Have looked at the documentation but all it seems to say is to put the files in the user_projections, which I've done.
I'm trying to use OSTN15_NTv2_OSGBtoETRS.gsb to convert EPSG:27700 to EPSG:3857.  I've downloaded the NTv2 grid files and added to the  user_projections folder, however it doesn't seem to have any affect when testing in the Reprojection console.
How do I 'make' GeoServer use the more accurate grid shift?


texts:
I've seen raster attribute tables mentioned in several places here at GIS.SE and at ESRI websites.

However, I don't understand what is a raster attribute table (RAT). The resources I found so far mostly explains how to build/delete a RAT and errors associated in the process, rather than what it is and what it's used for.

What I've figured out so far, is that RAT is only defined if the raster in question is categorical (in other words integer valued).

What is a raster attribute table? Can it have arbitrary number of columns and as complex as vector attribute tables?

How is RAT different from the attribute table of a vector dataset, and what is RAT typically used for?


texts: Availability of Frequent LiDAR Data in UK?

texts:
I am looking for LIDAR data of grasslands located all around UK with a spatial resolution of 1m and a frequency of at least once a week.

I've already had a look at the Environmental Agency but their data is far too infrequent.

The aim is to use the acquired data to monitor the biomass of a pasture field over different seasons.


texts: Export features to pdf using JSPDF and Openlayers 3?

texts:
I built a map with different type of data, WMS layer, Vector Layers ( features extracted them from GeoJson file).

Used Technologies :
- Openlayers 3
- JSPDF

How can i export just the selected feature from the map ? i parsed the whole of net asking about that, but without solution.


texts: How to create a x-y scatterplot from attributes in QGIS 2.0?

texts:
There are a few old packages that allow the creation of x-y scatter plots from attributes in vector data

(Scattergram for 1.7) http://plugins.qgis.org/plugins/scattergrammidentify/

(SPQR for very early editions) http://www.maths.lancs.ac.uk/~rowlings/Software/Spqr/

But does anyone know of a package for 2.0 that would do this?  or a way to rig up the old ones?

Thanks


texts: Does FME have spatial "Difference" like QGIS / "Erase" in ArcGIS and Mapinfo?

texts:
I am looking for an FME equivalent to the difference/erase tool in QGIS/ArcGIS/Mapinfo. Basically I want to see the spatial change between 2 datasets.
https://gis.stackexchange.com/questions/76839/distinction-between-difference-and-clip-tools-in-qgis#:~:text=Difference%20creates%20a%20new%20feature,not%20overlap%20the%20%22clipper%22.
An old process in FME is https://community.safe.com/s/question/0D54Q000080hcVvSAI/erasing-subtracting-and-merging-polygon-layers
Regards,


texts: Geoserver SLD creation using a local image

texts:
Using this documentation I have made the below SLD file, attempting to use an svg that I obtained from here I have moved the svg into /var/lib/tomcat7/webapps/geoserver/data/styles but the image will not render when previewing in the geoserver style editor (the return is a grey square).  I've gone through this SE post and tried to reverse engineer it and have also used a direct path the result are the same.  What am I missing?

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <!-- a Named Layer is the basic building block of an SLD document -->
  <NamedLayer>
    <Name>default_point</Name>
    <UserStyle>
    <!-- Styles can have names, titles and abstracts -->
      <Title>Default Point</Title>
      <Abstract>A sample style that draws a point</Abstract>
      <!-- FeatureTypeStyles describe how to render different features -->
      <!-- A FeatureTypeStyle for rendering points -->
      <FeatureTypeStyle>
        <Rule>
          <Name>rule1</Name>
          <Title>platform</Title>
          <Abstract>Platform icon</Abstract>
            <PointSymbolizer>
                <Graphic>
                    <ExternalGraphic>
                        <OnlineResource xlink:type="simple" xlink:href="26 SUSPENDED-OILSHOW.svg"/>
                        <Format>image/svg+xml</Format>
                    </ExternalGraphic>
                    <Size>20</Size>
                </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

I was able to get it to work with a direct path xlink:href="file:////var/lib/tomcat7/webapps/geoserver/data/styles/26%20SUSPENDED-OILSHOW.svg"  Why will it not work with a relative path?


texts: pyqgis 3 QgsVectorLayer.startEditing() doesn't trigger vertex display

texts:
In pyqgis 2 standalone application, the following code will trigger the display of vertex as red cross mark, but in pyqgis 3.0.3 standalone application, it doesn't:

layer.startEditing()
mapCanvas.refresh()

why is that? Do we have to somehow manually trigger the vertex markers to be displayed?


texts: QGIS 3 processing script is recognized as a model

texts:
I refer to the QGIS 3 (QGIS 3.4.4-1, 64Bit standalone installer) processing script shown below. If I copy the script, do Create New Script... from the processing toolbox, paste it to the empty script and save it as, let's say xlsx_vector_writer.py the script is somehow added to scripts, but as a model like this:

Note the Algorithm:ID 'script:model' which to me is a contradiction in itself, and should be 'giswg:xlsxvectorwriter', as stated in the source code below.

When I right-click it, it behaves like a script AND a model as well:

one may verify, that this is the script just saved by selecting 'Edit script...'.

Q: Is anyone able to reproduce this? What might it be that causes this weird behavior?

This is the script:

from PyQt5.QtCore import QCoreApplication
from qgis.core import *
import processing, datetime

class XLSXVectorWriterProcessingAlgorithm(QgsProcessingAlgorithm):
    INPUT = 'INPUT'
    MIDFIX = 'MIDFIX'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return XLSXVectorWriterProcessingAlgorithm()

    def name(self):
        return 'xlsxvectorwriter'

    def displayName(self):
        return self.tr('XLSX Vector Writer')

    def group(self):
        return self.tr('giswg')

    def groupId(self):
        return 'giswg'

    def shortHelpString(self):
        return self.tr("Export Vector as XLSX with filename as projectfolder_midfix_date.xlsx")

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.INPUT,
                self.tr('Input layer'),
                [QgsProcessing.TypeVectorAnyGeometry]
            )
        )

        self.addParameter(
            QgsProcessingParameterString(
                self.MIDFIX,
                self.tr('The middle of this fine output filename (fix)')
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        input_vec = self.parameterAsVectorLayer(
            parameters,
            self.INPUT,
            context
        )

        midfix = self.parameterAsString(
            parameters,
            self.MIDFIX,
            context
        )

        feedback.setProgressText('get project folder')
        projectfolder = QgsExpressionContextUtils.projectScope(QgsProject.instance()).variable('project_folder')
        prefix = projectfolder.split('\\')[-1].lower()
        feedback.setProgressText('get date')
        today = datetime.datetime.today().strftime('%Y%m%d')

        feedback.setProgressText('make filename')
        filename = '{0}_{1}_{2}.xlsx'.format(prefix, midfix, today)

        feedback.setProgressText('write xlsx')
        crs = QgsCoordinateReferenceSystem() #instead of None
        w = QgsVectorFileWriter.writeAsVectorFormat(input_vec, '{0}/{1}'.format(projectfolder, filename), 'utf-8', crs, 'XLSX')

        return {}


texts: OpenTripPlanner (Error: 406) No transit times available

texts:
I try using OpenTripPlanner to create a webmap with multimodal routing function.

I followed the instruction on: http://opentripplanner.readthedocs.io/en/latest/Basic-Usage/

I downloaded otp-1.1.0-shaded.jar from: https://repo1.maven.org/maven2/org/opentripplanner/otp/
I grabbed GTFS data from: https://code.google.com/archive/p/googletransitdatafeed/wikis/PublicFeeds.wiki (download for Berlin:  http://www.vbb.de/de/datei/gtfs-daten.zip)
Also got latest osm.pbf file from: http://download.geofabrik.de/ (download for Berlin: http://download.geofabrik.de/europe/germany/berlin-latest.osm.pbf)
Started the local server with: java -Xmx2G -jar otp-1.1.0-shaded.jar --build C:/user/name/otp --inMemory,
whereas berlin-gtfs.zip and berlin-osm.pbf are located in C:/user/name/otp

The server is starting and I can access the map via http://localhost:8080/
But when I want to use the routing function for transit, it always gives me the following error:

  No transit times available. The date may be past or too far in the
  future or there may not be transit service for your trip at the time
  you chose. (Error 406)

My thoughts are that the GTFS data are somehow broken or not useable for OpenTripPlanner.

So I used the GTFS Feed Validator (https://github.com/google/transitfeed/wiki/FeedValidator)

The Results of the GTFS FeedValidator:

Agencies: Verkehrsverbund Brandenburg-Berlin, ...
Routes: 1363
Stops: 12845
Trips: 204915
Shapes: 0
Effective: June 06, 2013 to December 14, 2013

Found these problems:
19528 warnings

87 Duplicate IDs
1 Expiration Date
1810 Invalid Values
5313 Overlapping Trips In Same Blocks
5 Stops Too Closes
12310 Too Fast Travels
2 Unrecognized Columns


texts: Leaflet-omnivore load multiple KMLs & center

texts:
I'm new with leaflet-omnivore and looking for the optimal way to load a number of KMLs (each of them containing a polygon) and center the map in their "bounds"

I have tried adding the loaded polygons to a feature group and then try to get the bounds of that in order to center the map. The polygons are showing on the map successfully but the map does not zoom/center

    var polygonGroup = new L.FeatureGroup();

    // Load kml and add to feature group
    $.each(kmls, function (i, kmlUrl) {
        var polygonLayer = omnivore.kml(kmlUrl);
        polygonGroup.addLayer(polygonLayer);
    });

    // Add group to map
    map.addLayer(polygonGroup);

    // Attempt to center map
    map.fitBounds(polygonGroup.getBounds());`

Can anyone point me to the right direction?


texts: I am searching for a way to lock the size of the circle markers in Leaflet for R, addCircleMarkers().  I need the markers to retain the size they are relative to the polygon they are within upon zooming.  At the moment, when you zoom, the markers stay the same size relative to your screen.


texts: Problems printing OpenLayers Plugin Maps To Scale

texts:
In the print composer I can print all maps to scale and satisfactory.
However, when I add the OpenLayers, usually Bing Aerial, to the project I cannot print it out to scale. I click the scale box, set it to the relevant scale.

The image below is a section of the actual size with OS mapping:

The next image is with the OpenLayers.

Both images cropped to remove client info. It is almost double the size but in-between scales. Is there any way to correct this?

Using Windows 7, QGIS 1.8.0.


texts: Is it possible, and if so how, to effectively classify vegetation types using false near infra red imagery?

texts:
Background : I use near infra-red imagery to delineate wetlands and water bodies with good success. It's also obvious to distinguish between new healthy vegetation and older growth and land areas with no vegetation but i need to devide vegetation into further groups.  To look at other examples, I downloaded the latest NRCAN (Natural Resources Canada) CANVEC (Canada Vector) vegetation data (polygons with attributes) and the data is incredibly detailed given the coverage (Canada) and scale (1:50000).  It is my understanding the classification was done with using Landsat 7 & 8 multispectral imagery in the near infra-red spectrum. I believe this data is also continually updated as new imagery becomes available.  There are at least 14 different vegetation types (screen capture with types below). How is this possible? It seems that to do my project manually would be possible but also an enormous undertaking at the detail / scale I need (approx 1:2000).

Questions:

How can I efficiently classify vegetation using custom false NIR high resolution imagery into similar categories?

Is there a fairly accurate automated process or is this done visually / manually?

What is a typical technique one may use to classify vegetation in this way?

CanVec Vegetation categories as classified from imagery:

Typical example of our custom imagery. I can delineate this area manually and perhaps work out some categories but this is only a fraction of our data and it would take VERY long time to finish the entire area.


texts: Arcmap 10.3 addin is not working on 10.2.1

texts:
I have made an add-in using Visual Studio 2013 and ArcObjects 10.3 .Net SDK. It is running fine on 10.3. I want to run this on 10.2.1.  I have change the following settings:

<Targets>
<Target name="Desktop" version="10.2" />
</Targets>

It is installed but when I click on it ArcMap crashes without any error. I have tried on two different systems also. I have also tried by creating a blank form on button click to check but it didn't work and crashes.


texts:
I have a table that has an ID, point geometry and time. It's basically GPS coords and some other data. I want to be able to select all rows that interact (come close in space and time) to my target ID.

For example I get my target track, John's points, with a simple query and then I can create a line from those points, buffer it and fetch all other rows that fall within that buffered track.

In the image below I have all of John's points for a month and all points from other ID's for the same month that fall within 10 km of his track. This works fine if I limit to a day or two.

The problem I'm experiencing is that as the time limit increases I'm getting more and more bad matches; the points of John's track that happen on the 15th are surrounded by tracks happening on the 1st or 30th.

So how should I be doing this? Create CTE of John's points and then some kind of moving window over that to only fetch points that are within 1 km and 1 hour of John's point X?

I'm not sure what tools to use...

SQL I use to get image (although I colour John's points differently):

WITH target_track AS
    (SELECT DISTINCT ON (name)
  name,
  ST_MakeLine(gps_point) AS track,
  ST_Buffer(ST_MakeLine(POSITION), 0.1) AS buff_track
     FROM gps_history
    WHERE
    name = 'john'
    AND
     event_time BETWEEN '2019-01-01' and '2019-01-30'
    GROUP BY name)

SELECT gps_history.* FROM
gps_history, target_track
WHERE
ST_Within(gps_point, target_track.buff_track)
AND
event_time BETWEEN '2019-01-01' and '2019-01-30'

I use group by and distinct to get a single row for the target_track CTE and then use that to find all the points that are in that time/space region.


texts: Add Features from KML to a GeoPackage using GeoTools

texts:
Using GeoTools v13.2 I need to add features to a GeoPackage from KML. So far I can get the features in a SimpleFeatureCollection which can be added to the GeoPackage:

public void addFeaturesToGeoPackage(InputStream featuresKML, GeoPackage geopackage)
{
    Parser parser = new Parser(new KMLConfiguration());

    SimpleFeature f = (SimpleFeature) parser.parse(featuresKML);
    Collection placemarks = (Collection) f.getAttribute("Feature");
    List<SimpleFeature> featureList = new ArrayList<SimpleFeature>(placemarks);

    //Get SimpleFeatureCollection
    SimpleFeatureCollection featureCollection = DataUtilities.collection(features);

    //Get the feature schema
    SimpleFeatureType schema = featureCollection.getSchema();

    //TODO: Create FeatureEntry from schema
    FeatureEntry featureEntry = createFeatureEntry(schema);

    //Create the geometry table
    geopackage.create(featureEntry, schema);

    //Add all the features
    geopackage.add(featureEntry, featureCollection);

}

I'm not sure how to create the FeatureEntry as what I'm currently using gives an arbitrary name to the "geometryColumnName" attribute which causes a java.lang.IllegalArgumentException:

Caused by: java.lang.IllegalArgumentException: Geometry column featureVectors does not exist in schema

This is what I was able to come up with based on the documentation and examples I could find through some Google searches but I would like to know the right way to accomplish this/get the FeatureEntry.


texts: Converting shapefile from projected coordinates using pyproj and fiona

texts:
I am new to GIS. I am trying to map out wards for the city of Milwuakee using shapefiles found on their county website county website. I am following Converting projected coordinates to lat/lon using Python? with some success. My code is:
from pyproj import Proj, transform
# wisconsing EPSG:32054
# epsg:4326 is for the entire world, wgs 84...not obvious
inProj = Proj(init='epsg:32054')
outProj = Proj(init='epsg:4326')
x1,y1 = 2560131.496875003, 406816.434375003
x2,y2 = transform(inProj,outProj,x1,y1)
print(x2,y2)

with output,
-65.70220967836329 43.08590211722421

Problem is this is wrong. The lon/lat for Milwaukee is -87.863984 and 42.920816.
Secondly, how can I do this programmatically for the entire shapefile. I would like to plot this into basemap. When I try to follow Shapefile projection for Matplotlib Basemap I get an error
Code:
with fiona.open("ward2012/ward.shp") as shp:
    ori = Proj(init='epsg:32054' ),
    dest= Proj(init='EPSG:4326',preserve_units=True)
    with fiona.open('ward2012/MKE_wards_lat_lon.shp', 'w', 'ESRI Shapefile', shp.schema.copy(), crs=from_epsg(4326))as output:
        for point in shp:
            x,y =  point['geometry']['coordinates']
            point['geometry']['coordinates'] = transform(ori, dest,x,y)
            output.write(point)

error:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-139-a5079ab39f99> in <module>()
      4     with fiona.open('ward2012/MKE_wards_lat_lon.shp', 'w', 'ESRI Shapefile', shp.schema.copy(), crs=from_epsg(4326))as output:
      5         for point in shp:
----> 6             x,y =  point['geometry']['coordinates']
      7             point['geometry']['coordinates'] = transform(ori, dest,x,y)
      8             output.write(point)

ValueError: not enough values to unpack (expected 2, got 1)


texts: Solving route in multimodal network dataset does not use all forms of transport available?

texts:
After having completed Exercise 2 in the official Esri Network Analyst Tutorial, I decided to find a route in the multimodal dataset, to see how it functions. After adding several network locations along an intentionally long route, I pressed Solve. I was quite surprised to find that the route chosen took walking along the streets as the only transportation mode available (despite railway stations and lines being included in the dataset).

I placed my origin and destination points quite far apart close to metro entrances, so I assumed that the directions would include walking to that first metro entrance, travelling on the train to the last metro entrance, exiting the station and walking to the destination point. The route ArcMap gave me indicated I would have to walk 3 hrs 46 mins along the streets to get to my destination!

I looked at the Analysis Settings, but I am not familiar enough with Network Analyst to know if something needs to be changed in order to get the results I want.

How do I ensure that the route chosen uses all available forms of transport to find the shortest route?


texts: Counting NoData values in raster using ArcGIS Desktop?

texts:
I have created a raster in ArcGIS Desktop 10 with a large number of holes (NoData) in it.

I want to fill the holes using something like the the following:

rOut = arcpy.sa.Con(arcpy.sa.IsNull(rGround),
arcpy.sa.FocalStatistics(rGround, arcpy.sa.NbrCircle(FillRadius,"CELL"), "MEAN"),
rGround)

It is desirable to fill the holes by averaging across a circle with a relatively small radius, but a small radius leaves larger holes unfilled, so I need to iterate this operation.

To do that I would like to be able to get the number of NoData values in my raster, once it reaches zero I can stop filling.