texts: Can I split up v.in.ogr into multiple v.clean commands?

texts:
I'm trying to load a large set of polygons into GRASS with "v.in.ogr" and am running out of memory at the "Break Boundaries" stage.  I'm wondering if I can break into steps all of the different things that v.in.ogr does to a data set when it is run without the "-c" flag?  For example, I know that it must run the equivalent of "v.clean tool=bpol" and "v.clean tool=rmdupl", but I'm not exactly sure what to run after those 2 to complete the polygon cleaning that normally happens.  I figure that I might be able to save some memory if I don't have to run all the operations in one go.

Does this makes sense?

Thanks.

Ok, so it was like digging for gold to find this, but in a presentation given by some people who seem to know GRASS pretty well, they said that the following steps were needed to clean polygon geometry after import:

Cleaning procedures for area import

break polygons
remove duplicates
break boundaries
remove duplicates
clean boundaries at nodes
change dangles to lines
remove bridges

Link to presentation here:  http://geostat-course.org/system/files/GRASS_geostat_landau_2011_intro.pdf

Wait, no... there's more to this apparently.  I tried those steps listed above (or attempted to at least). I think I'm worse off now than when I started.

v.in.ogr -c
Importing map 52 features...
Number of boundaries: 52
Number of centroids: 52
Number of areas: 42
Number of isles: 42
Number of incorrect boundaries: 10
Number of centroids outside area: 10

1. break polygons
  v.clean type=boundary tool=bpol
Number of boundaries: 270
Number of centroids: 52
Number of areas: 2
Number of isles: 4
Number of incorrect boundaries: 268
Number of centroids outside area: 50

2. remove duplicates
  v.clean type=boundary tool=rmdupl
Number of boundaries: 190
Number of centroids: 52
Number of areas: 40
Number of isles: 4
Number of incorrect boundaries: 126
Number of centroids outside area: 12

3. break boundaries
  v.clean type=boundary tool=break
Number of boundaries: 5147
Number of centroids: 52
Number of areas: 151
Number of isles: 0
Number of incorrect boundaries: 5132
Number of centroids outside area: 40
Number of areas without centroid: 139

4. remove duplicates
  v.clean type=boundary tool=rmdupl
Number of boundaries: 4803
Number of centroids: 52
Number of areas: 153
Number of isles: 0
Number of incorrect boundaries: 4781
Number of centroids outside area: 38
Number of areas without centroid: 139

5. clean boundaries at nodes
  (No idea what-the-flip this means.  Clean with what?  Do it by hand?)

6. change dangles to lines
   v.clean type=boundary tool=chdangle thresh=-1
Number of boundaries: 4803
Number of centroids: 52
Number of areas: 153
Number of isles: 0
Number of incorrect boundaries: 4781
Number of centroids outside area: 38
Number of areas without centroid: 139

7. remove bridges
  v.clean type=boundary tool=rmbridge in=eq1_6 out=eq1_7
Number of boundaries: 4420
Number of centroids: 52
Number of areas: 153
Number of isles: 22
Number of incorrect boundaries: 4352
Number of centroids outside area: 38
Number of areas without centroid: 139

Guess I'm off to read the source for v.in.ogr now.


texts: Extract By Mask - Raster Horizontal Shift

texts:
Is there a way to snap a raster/DEM to a mosaic dataset?

I am using the Extract By Mask tool in ArcPro on a mosaic dataset. The output raster's cell values do not match the original mosaic dataset it came from. The cells seem to have shifted after the Extract By Mask Tool was run. The only way I know how to fix this problem is by changing "Snap Raster" in Environments, but have not yet found a way to snap a raster to a mosaic dataset.


texts: Appending Data to FeatureLayer in ArcGIS API for JavaScript?

texts:
I am working on a web app where I am making calls to a database api and it returns a JSON Object. I've managed to do so successfully as well as any manipulations that need be done.

Is there a way to add new fields to a FeatureLayer and the subsequently populate with the relevant data?

I do not necessarily want to edit the feature layer, just attach data on the fly so that I can then use it for analysis such as heat maps and what not. Have not really found anything concrete in the api documentation.

Ive tried the following code that results in an error. Fairly certain it is because i am trying to set the value before I add the Field.

featureLayer = new FeatureLayer(....);
fieldName = x;
value = y;

featureLayer.graphics.attributes.fieldName = value;


texts: Querying Geospatial data in Athena database, getting rows in range of 5 miles from specified coordinates

texts:
I have database table Building with these columns: name, lat, lng

How can I get all Buildings in range of 5 miles from specified coordinates, for example these:

-84.38653999999998

33.72024

?

https://docs.aws.amazon.com/athena/latest/ug/geospatial-functions-list.html


texts: QGIS - Add layer that is Editable and with Add Polygon Feature selected

texts:
I want to add a layer in QGIS that is editable and has Add Polygon Feature selected (by default). I can easily accomplish the editable part by this command:

layer = QgsVectorLayer("Polygon?crs=epsg:4326", "temporary_points", "memory")
layer.startEditing()
#Something more here?
QgsProject.instance().addMapLayer(layer)

A bonus question is why it takes extremely long time to get the saving popup when I have that as a part of my code? (and none when I run it from the Python Console)


texts: How to change a schema name within a MXD with ArcObject .NET?

texts:
I have a MXD file where all layers point to an Oracle Spatial database. Recently our DBA changed the schema name. So now my MXD file has broken links, each time I open it, all my layers have a exclamation mark beside them.

In order to bulk change the new schema name in the MXD file, I've written an .NET ArcObjects application. But since all my layer data links are broken in my MXD, I do not have the capability to get a reference on IFeatureClass from a IFeatureLayer.

I've found on the web that is normal if a IFeatureLayer is broken within the MXD that the IFeatureClass always returns null.

Actually, all MXD's layers point on the schema "A" and now I want them to point on schema "B". How to do that programmatically?

Can you help me to get the good object that should give me access to change the QueryDef of my layers? Here is the code that I've written until now:

         //Exit if no map document is selected
        string sFilePath = openFileDialog1.FileName;
        if (sFilePath == "") return;

        IMapDocument ipMapDoc;
        ipMapDoc = new MapDocumentClass();

        ipMapDoc.Open(sFilePath);

        IMap map = ipMapDoc.get_Map(0);

        IEnumLayer layers = map.get_Layers();
        ILayer layer;
        while ((layer = layers.Next()) != null)
        {
            if (layer is IFeatureLayer)
            {
                IFeatureLayer featureLayer = layer as IFeatureLayer;
                IFeatureClass featureClass = featureLayer as IFeatureClass;

                if (featureClass == null)
                {
                    Console.WriteLine("FeatureClass where FeatureLayer name is " + layer.Name + " is null");
                }

            }
        }


texts: Rotating point symbol to underlying line for QGIS 3?

texts:
The expression below worked perfectly under QGIS 2.18. One simply entered it as an expression for the valve symbol and all valves rotated based on the underlying pipe feature. However it no longer works on QGIS 3.2 due to changes in the qgsGeometry functions.

How do I do this in QGIS 3?

degrees(azimuth(start_point( geom_from_wkt( geomnearest( 'WaterPipes02_Vermont', '$geometry'))), end_point( geom_from_wkt( geomnearest( 'WaterPipes02_Vermont', '$geometry')))))+90


texts: One-to-Many table relate viewed through ArcGIS Server

texts:
I want to create a Table Relate with a one-to-many relationship between a layer and a non-spatial table based on a common field.

My question is, if I publish this data to ArcGIS Server, will an identify operation on that layer return all of the related values in the non-spatial table?

I want to try this, but I don't have sample data that meets these requirements (if anyone knows where I can find some, it would be very helpful).

Here is an example from ESRI of an Identify operation in ArcMap on two tables with a one-to-many relationship.  I need to know if the results in ArcGIS Server would be consistant.

I am using ArcGIS 10 and SQL Server 2008 R2.


texts: I have been using QGIS with an Apple Magic Mouse.  Despite disabling all the point and click features and gestures in the systems preferences, the maps constantly zoom in and out extremely and uncontrollably with the merest touch of the cursor to the screen.  What could be causing this?


texts: I want to calculate a cost surface (anisotropic) using the SAGA tool. I am using a 'slope' raster for the cost grid entry, an 'aspect' one for the direction of max cost. I have to put the destination points in, but it is a raster layer. So I tried to rasterize them using the rasterize (raster to vector) tool, but when I run it nothing happens. I tried to modify the size of the pixels, but again, the layer I end up with is non-existent. I am wondering if I what I am doing is correct to end up with this anisotropic surface, and also, can I perhaps use the interpolation tool to rasterize the points?


texts: PostGIS: MultipolygonZ to floating point raster

texts:
I am trying to convert a table of MultipolygonZ geometries to a floating point raster where the pixel values of the output raster will be derived from the Z values of the MultipolygonZ geometries. The ultimate goal is to perform map algebra between the resultant floating point raster and a DEM (already loaded in PostGIS). I have already tried to perform the analysis directly between the MultipolygonZ geometries and the DEM without success - see PostGIS: intersection of a MultipolygonZ and a floating point raster.

Here is an example of one of the geometries that I am working with - it is essentially an inclined plane that represents water levels of a flood plain interpolated between two cross sections:

SELECT ST_GeomFromText('MULTIPOLYGON Z (((277378.7412 5048118.3481 22.46,278037.7398 5047686.44410001 22.46,277929.677100001 5047698.6316 22.46,277380.212900001 5048034.458 22.46,277105.837500883 5048202.1531993 22.46,277167.5633 5048248.2075 22.46,277173.661101116 5048252.75649908 22.46,277378.7412 5048118.3481 22.46)))')

Is this possible using the ST_AsRaster(), ST_ExtractToRaster() or any other functions/methods within PostGIS? I am able to create a raster using ST_AsRaster(), but I can't get the output raster to contain the interpolated Z values.


texts: Python user defined function written for ArcGIS Pro returns None when it should return an address?

texts:
I am rewriting a 2.7 Python script for compatibility with ArcGIS Pro and organizing it with user defined functions so that I can debug portions of the script as I go. I am new to the process, but I thought it would be a good idea to pull the function out and test it alone. It's supposed to search a database and find an address location based on a SQL query. Right now the location is hardcoded for testing. When I run the script and print the selection variable, I get "None". Any suggestions or advice for debugging this?

import pdb
import arcpy
import os
import string
from io import StringIO
from io import BytesIO
import getpass

def select_address():
    workspace = "G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\BSC_Maps\\AutoLayerFiles.gdb"
    arcpy.env.workspace = workspace
    arcpy.env.overwritOutput = True
    location = '1900 McCall Rd'#arcpy.GetParameterAsText(0)
    expression = "FULL_STREET_NAME LIKE '" + location + "%'" #this is an SQL expression and can be modified as necessary
    aprx = arcpy.mp.ArcGISProject("G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\BSC Maps\\BSC_Template_Auto_Generate_PRO_v1.aprx")
    for m in aprx.listMaps():
        #print(m)
        for lyr in m.listLayers():
            #print(lyr)
            if lyr.name == "LOCATION.addresspoint":
                CoaAddresses = lyr
                #print(CoaAddresses)
                selection = arcpy.SelectLayerByAttribute_management(CoaAddresses,"NEW_SELECTION", expression)
                print(selection)
                for row in arcpy.SearchCursor(selection): #retrieve address of selected feature
                    selectedAddress = row.FULL_STREET_NAME
                    return(selectedAddress)

selection = select_address()
print(selection)

EDIT: I have tried applying the suggestions made in the comments and tried the following, but I still get None.

import pdb
import arcpy
import os
import string
from io import StringIO
from io import BytesIO
import getpass
workspace = "G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\BSC_Maps\\AutoLayerFiles.gdb"
arcpy.env.workspace = workspace
arcpy.env.overwritOutput = True
location = '7808 HART LN'#arcpy.GetParameterAsText(0)
aprx = arcpy.mp.ArcGISProject("G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\BSC Maps\\BSC_Template_Auto_Generate_PRO_v1.aprx")

def select_address(address):

    expression = "FULL_STREET_NAME LIKE '" + address + "%'" #this is an SQL expression and can be modified as necessary
    for m in aprx.listMaps():
        #print(m)
        for lyr in m.listLayers():
            if lyr.name == "LOCATION.addresspoint":
                CoaAddresses = lyr
            #print(lyr)
                with arcpy.da.SearchCursor(CoaAddresses, ['FULL_STREET_NAME'], where_clause=expression) as cursor:
                for row in cursor:
                    print(row[1])

selection = select_address(location)
print(selection)


texts: Is there an automatic way to find buildings in a raster file?

texts:
Am trying to instal BREC plugin in QGIS 2.4 and 1.8 for the purpose of automatically detect polygons (buildings) from a raster (TIF) file but i can't find that plugin, its not available in plugin library. How can I get that plugin in QGIS and its working steps?

Is there any other way to automatically find the polygons (buildings) from a raster file (TIF) in QGIS?


texts: XY event layer doesn't work with ArcGIS Server 10.1

texts:
I have a geoprocessing tool written in Python that contains operations on an XY event layer. It worked perfectly when running from ArcCatalog. When published as geoprocessing service on ArcGIS Server 10.1, however, it never worked. The JSON response when executed as REST service was

"error": {
  "code": 400,
  "message": "Unable to complete operation.",
  "details": [
   "Error executing tool.: ERROR 999999: Error executing function.\nThe XY event source name has not been setup correctly.\nFailed to execute (CopyFeatures).\n\nFailed to execute (NewCuringModelExpPy).\nFailed to execute (New_Curing_Model_ExpPy)."
  ]
 }

I have been struggling with this for days. Has anyone else ever come across this issue before or how was that solved?

Thanks!

Alex


texts:
I am developing tool on ArcGIS ArcMap 10.0 with .NET and Python technology.

In TOPOTORASTER tool , I give a line feature class. the line feature class path something like this:

inputOne = "C:\Users\DELL\AppData\Local\Temp\tmpA378\tmpA379.gdb\PerpendicularLine"

inPointElevations = inputOne
print os.path.dirname(inputOne)
pathAfterJoin = os.path.join(inputOne)

arcpy.TopoToRaster_3d(inPointElevations,"D:\\Temp\\topo1234")

but in python IDE (pycharm) I have seen that feature class path is not valid.
feature class path is broken and it becomes like:...C:\Users\DELL\AppData\Local\Temp    mpA378  mpA379.gdb ....

space is taken tempA378 folder ...

How to efficiently give input feature class path ??

I am trying like :

arcpy.TopoToRaster_3d(r'inPointElevations',"D:\\Temp\\topo1234")

but not work..


texts: Leaflet Layergroup - hide one layergroup on first page render

texts:
Using Leaflet.js, I have two LayerGroups:

var LayerSeventeen = new L.LayerGroup();
var LayerEighteen = new L.LayerGroup();

And their respective markers bound to those layersGroups:

var marker_reptile_house =  L.marker([33.92641, -83.387], {icon: buildingPin}).addTo(map);
LayerSeventeen.addLayer(marker_reptile_house);

var marker_gator_slider_deck =  L.marker([33.92669, -83.38705], {icon: animalPin}).addTo(map);
LayerEighteen.addLayer(marker_gator_slider_deck);

And the javascript to run it in a script named "using-leaflet.js":

// more stuff above...
map.setView([33.927106, -83.387013], 17)

map.on('zoomend', onZoomend);
function onZoomend(){
    if(map.getZoom()>=17){
    map.addLayer(LayerSeventeen);
    };
    if(map.getZoom()>=18){
    map.addLayer(LayerEighteen);
    };
    if(map.getZoom()<17) {
        map.removeLayer(LayerSeventeen);
    };
    if(map.getZoom()<18) {
        map.removeLayer(LayerEighteen);
    };
};

// more stuff below...

This all works as described in this example (not mine):
http://geohacker.github.io/leaflet-snippets/

Except all of the markers show up when the page first renders at zoom level 17. No zooming is involved.

How do I get it to render only LayerSeventeen on the first page render?


texts:
I am new to spatstat package and I am trying to make a point patten map with a shape file and some case points. I was able to plot the shape file plot but I just don't know how to make the points on the map plot through the spastat object even I tried points() function or plot() function. What am I missing?
pacman::p_load(sf, sp, spatstat)

# shape file
manchester_ward <- st_read("https://raw.githubusercontent.com/RUMgroup/Spatial-data-in-R/master/rumgroup/data/wards.geojson")
Fallowfield <- manchester_ward %>% dplyr::filter(wd16nm == "Fallowfield")
Fallowfield_pp <- as.owin(Fallowfield)
plot(Fallowfield_pp)

# points
crime <- structure(list(crime_type = c("Burglary", "Burglary", "Burglary",
"Burglary", "Burglary", "Burglary"), geometry = structure(list(
    structure(c(-2.242511, 53.44558), class = c("XY", "POINT",
    "sfg")), structure(c(-2.237017, 53.445735), class = c("XY",
    "POINT", "sfg")), structure(c(-2.2334, 53.442408), class = c("XY",
    "POINT", "sfg")), structure(c(-2.242048, 53.448871), class = c("XY",
    "POINT", "sfg")), structure(c(-2.242048, 53.448871), class = c("XY",
    "POINT", "sfg")), structure(c(-2.243242, 53.444339), class = c("XY",
    "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = -2.243242,
ymin = 53.442408, xmax = -2.2334, ymax = 53.448871), class = "bbox"), crs = structure(list(
    input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(crime_type = 1L), .Label = c("constant",
"aggregate", "identity"), class = "factor"), row.names = c(NA,
6L), class = c("sf", "data.table", "data.frame"))

crime_pt <- matrix(unlist(crime$geometry), ncol = 2, byrow = T)
k1 <- ppp(x = crime_pt[,1], y = crime_pt[,2],
          window = Fallowfield_pp, check = T)
plot(k1)


texts:
I would like to use OSGeo4W shell with Python3 but when typing python3 I get the following error:

Fatal Python error: Py_Initialize: unable to load the file system codec
  File "C:\OSGEO4~1\apps\Python27\lib\encodings\__init__.py", line 123
    raise CodecRegistryError,\
                            ^
SyntaxError: invalid syntax

How can I use Python3?


texts: Which is the model semivariogram for binary data?

texts: I have a binary data (1,0) to calculate spatial distribution of soil bacteria. I have to use indicator kriging (presence/absence) and fitting model of semivariogram. Then, which is the model semivariogram for binary data?


texts:
I have some tiff files and I'm trying to generate a vrt index with gdal tool, gdalbuildvrt.

I run this command:

gdalbuildvrt orto_index.vrt -srcnodata "0 0 0" -vrtnodata "0 0 0" *.tif

and throws the next warning:

Warning 6: gdalbuildvrt does not support rotated geo transforms.
Skipping 2335-C.tif

So there are some tiff that are skipped.

If we compare the information about a skipped file with another that it has not skipped, we could see the next difference:

In the skipped file, we can see the next line:

GeoTransform =
  345770.7884398972, 0.6000000238418568, 3.473043927400512e-15
  6139634.558238117, 1.757016748640284e-14, -0.6000000238417672

In the non skipped file, the other one:

Origin = (223499.236031105130678,6213152.835421450436115)
Pixel Size = (0.600000023841858,-0.600000023841858)

How can I make to get a vrt file including all files, I mean, without skip some files?


texts: Selecting connected lines between start and end point in QGIS?

texts:
Is there a method in QGIS to select all connected lines that are between a start and end feature (point or selected lines)?

I imagine a workflow going something like this:

Choose the starting line feature
Choose the ending line feature
Run a process that will select all features (between the start and end) that are connected.

The end result would look something highlighted yellow lines in the picture below, assuming I had pre-selected lines at the beginning and end of the now highlighted arcs.

This is an alternative method: How do I select the longest connected lines from PostGIS ST_ApproximateMedialAxis?


texts: What is the difference between Sub and a Function in Mapbasic?

texts:
In MapBasic it seems like Functions can have side effects so I don't really see the difference between using the following:

Declare Sub MyProcedure (myParameter As String) and then use myParameter after Call MyProcedure(myParameter)

or

Declare Function MyProcedure () As string and do something like myParameter = MyProcedure()

Am I missing something?


texts: Why I cannot export OSM map as a .pdf/.jpeg in QGIS 2.4?

texts:
I have installed on my Mac OSX 10.9.5 the version QGIS 2.4.

My problem is that I cannot export any OpenStreetMap with the Composer and I do not understand why. The thing is that I imported the OSM map with the OpenLayers plugin... everything normal, OK. Then I started to edit some vectors... and then I tried to export the map with the Composer... but it just does not work! I can load everything there, but when I press the export buttons, then the window becomes blank. And it creates a new .pdf for example but with 0 bytes...

Does anyone have an idea about what I did wrong, or what is wrong? I may have missed some key details that I did not know that should be done...but still...Is there any chance to have the background of the OSM map exported somehow like .pdf? What I can do?

Or better, another question, is it even possible to export the OSM imported with OpenLayers Plugin with the Composer?... or I should download the data on my computer first?


texts: QGIS Tile Layers Plugin crashing when adding custom map

texts:
I am trying to add a base map to QGIS using the TileLayerPlugin. However, whenever I execute the plugin QGIS crashes. Specifically, I am trying to add this basemap to QGIS.

I created a .tsv file in TextWrangler that reads as follows:

Thunderforest   Thunderforest   http://{s}.tile.thunderforest.com/pioneer/{z}/{x}/{y}.png

I went to Web >> TileLayerPlugin >> Add Tile Layer and then added the directory with the .tsv file. I select the .tsv file (see screenshot) and click Add. When I do that QGIS crashes. I also tried inserting this code into the tsv file with no success

Thunderforest   Thunderforest   https://[abc].tile.thunderforest.com/cycle/{z}/{x}/{y}.png


texts:
I'm using GRASS on OS X. After installing GRASS I can no longer open my terminal to a prompt.

On opening the terminal, the following line runs by itselt

'/Applications/GRASS-6.4.app/Contents/MacOS/grass.sh'; exit

And GRASS starts but itself. I can't find how to stop terminal from doing this. I just want my shell prompt back.


texts:
I'm trying to install MapServer version 7.0.6. I also want php MapScript so I'm doing it by compiling from source. I've installed two versions of PHP (5.6 and 7) on my system (Linux Mint 18.1) but I've downgraded to 5.6 since as far as I'm concerned version 7 is still in development.  So everything worked well until I got this error:

[ 82%] Building C object
mapscript/php/CMakeFiles/php_mapscript.dir/mapscript_error.c.o
In file included from
/usr/local/src/mapserver/mapscript/php/php_mapscript.h:37:0,
              from
/usr/local/src/mapserver/mapscript/php/mapscript_error.c:32:
/usr/local/src/mapserver/mapscript/php/php_mapscript_util.h:267:1:
error: unknown type name ‘zend_object_value’
zend_object_value mapscript_object_new(zend_object *zobj,
 ^
In file included from
/usr/local/src/mapserver/mapscript/php/mapscript_error.c:32:0:
/usr/local/src/mapserver/mapscript/php/php_mapscript.h:461:8: error:
unknown type name ‘zend_object_value’
extern zend_object_value mapscript_object_new(zend_object *zobj,
zend_class_entry *ce,
    ^
/usr/local/src/mapserver/mapscript/php/php_mapscript.h:463:8: error:
unknown type name ‘zend_object_value’
 extern zend_object_value mapscript_object_new_ex(zend_object *zobj,
zend_class_entry *ce,
    ^
/usr/local/src/mapserver/mapscript/php/mapscript_error.c: In function
‘mapscript_throw_exception’:
/usr/local/src/mapserver/mapscript/php/mapscript_error.c:53:10:
warning: return from incompatible pointer type [-Wincompatible-
pointer-types]
return zend_throw_exception(mapscript_ce_mapscriptexception, message,
0 TSRMLS_CC);
       ^
/usr/local/src/mapserver/mapscript/php/mapscript_error.c: In function
‘zm_startup_mapscript_error’:
/usr/local/src/mapserver/mapscript/php/mapscript_error.c:114:37:
error: too many arguments to function
‘zend_register_internal_class_ex’
 mapscript_ce_mapscriptexception =
 zend_register_internal_class_ex(&ce,
 zend_exception_get_default(TSRMLS_C), "Exception" TSRMLS_CC);
                                 ^
 In file included from /usr/include/php/20160303/main/php.h:40:0,
             from
 /usr/local/src/mapserver/mapscript/php/php_mapscript.h:35,
             from
 /usr/local/src/mapserver/mapscript/php/mapscript_error.c:32:
 /usr/include/php/20160303/Zend/zend_API.h:292:28: note: declared here
  ZEND_API zend_class_entry
  *zend_register_internal_class_ex(zend_class_entry *class_entry,
  zend_class_entry *parent_ce);
                          ^
 mapscript/php/CMakeFiles/php_mapscript.dir/build.make:62: recipe for
 target
'mapscript/php/CMakeFiles/php_mapscript.dir/mapscript_error.c.o'
 failed
 make[2]: ***
 [mapscript/php/CMakeFiles/php_mapscript.dir/mapscript_error.c.o]
 Error 1
 CMakeFiles/Makefile2:547: recipe for target
 'mapscript/php/CMakeFiles/php_mapscript.dir/all' failed
 make[1]: *** [mapscript/php/CMakeFiles/php_mapscript.dir/all] Error 2
 Makefile:127: recipe for target 'all' failed
 make: *** [all] Error 2

I've also installed php5.6-dev on my system.

Is there any fix for this?


texts: What GIS application can open and use .osm files?

texts:
A quick question that I don't see an answer to elsewhere: if I download a .osm dataset from Planet.osm (for example) which GIS applications could consume/display this data without translation?

For example, this page says ArcGIS 10 "can use OpenStreetMap as a basemap" - but I'm not sure if that means it will open .osm files or read the data from a service?

So, are there any GIS apps (proprietary or open source) which will read .osm files?

Or perhaps taking a step back: is OSM data meant to be consumed like this, or is OSM meant to be used mostly for dynamic display of mapping on web pages?


texts: Exporting OpenLayers.Layer.WMS to GeoJSON / KML

texts:
I'm looking the way to exporting an especific region of a WMS layer to some XML format (GeoJSON, GML, KML) using OpenLayers and Geoserver. I've found how to export features objects with OpenLayers.Format.GeoJSON.write(), but this method only allow converting feature objects or geometry objects.

An example of my idea: I'm viewing a map of a certain country with a couple of WMS layers with some data. Then I zoom to a particular state of that country and I want to export what I'm viewing to GeoJSON or KML for use that in a desktop GIS application (like QGIS).

Does anyone know how to do this?

Thanks.


texts: How to calculate the difference in elevation between two datasets of points?

texts:
I have two datasets of different points in the same coordinate system, which describe the same line (road), so I want to see the difference in elevation between them.

I've been told that I should transform them to lines and then find the intersection points and calculate the difference in elevation there. But how can I do it in QGIS? Is there any plugins that I can use to somehow present the points as a route and then compare the two routes?

The first dataset derived from NRTK (vrs) measurements during a kinematic positioning on a road (1518 points). The second dataset derived from measurements by total station 30 years ago (100 points) so there is no common point between the datasets. I need to describe the elevation difference of the road, so I have to find a way to compare the datasets.
Both datasets represented as E,N (projected), H (orthometric).


texts:
I am trying to run the Calculate Adjacent Fields ToolBox on one of my Shapefiles, but the process is taking hours to run, and I have yet to complete the process. My table is 18,000 rows, and so far the most I've gotten through was 2,000 after 5 hours.  I tried creating and index on the field I want to run the Calculate Adjacent Fields on, (MN_80).   But that does not seem to improve performance.

Has anyone run into this issue before, and what were some solutions you found to speed up the process?


texts: Find shortest path and save it with geometry

texts:
With the guide of  Network analysis library

I found the shortest route on the network I'm analyzing. The script I used is as follows:

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(594869.256,4279586.388z)
pStop = QgsPoint(594011.309,4279760.098)

tiedPoints = director.makeGraph(builder, [pStart, pStop])
graph = builder.graph()

tStart = tiedPoints[0]
tStop = tiedPoints[1]

idStart = graph.findVertex(tStart)
tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)

idStart = tree.findVertex(tStart)
idStop = tree.findVertex(tStop)

if idStop == -1:
  print("Path not found")
else:
  p = []
  while (idStart != idStop):
    l = tree.vertex(idStop).inArc()
    if len(l) == 0:
      break
    e = tree.arc(l[0])
    p.insert(0, tree.vertex(e.inVertex()).point())
    idStop = e.outVertex()

  p.insert(0, tStart)
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor(Qt.green)

  for pnt in p:
    rb.addPoint(pnt)

Unfortunately, the result I get is a path highlighted in a different color. Is there a way to save the path as a shapefile containing the length of the path as information? Or, even better, save this information in a text file?

I'm following this suggestion:

lastPoint = None
d = 0
for point in p:
    if lastPoint is not None:
         dx = point.x - lastPoint.x
         dy = point.y - lastPoint.y
         d += sqrt(dx*dx+dy*dy)
    lastPoint = point
print "distance is "+str(d/1000.0)+"km"

but I get an error:
File "path/to/scipt.py", line 52, in
    dx = point.x - lastPoint.x
TypeError: unsupported operand type(s) for -: 'builtin_function_or_method' and 'builtin_function_or_method'


texts: QGIS Random Points in Polygon PRESERVING Attribute Table?

texts: I have a complicated shapefile of overlapping polygons. I need to generate random points inside each polygon, which is easily accomplished using Vector -> Research Tools -> Random Points Inside Polygons, but this obliterates the attribute table of the generated point shapefile which is not okay. Under more simple conditions, I could use NNJoin to rejoin the attributes, but because of the multiple overlapping polygons this is not acceptable as I cannot be sure the nearest polygon is actually the correct source of data. Can anyone suggest a method to generate a single random point in each polygon that has a 1:1 (ie., perfect) transfer of the attributes from the polygon to the random point, for each polygon?


texts:
I'm trying to recruit local police departments into a new GIS certificate that I am teaching in Florence South Carolina.  None of the Police departments have any experience with GIS, but they are all interested in the technology.

I'm making an image to advertise to them what they could learn.

Question: what are the hot topics in GIS for Police?  What should I include that would really catch their attention.


texts:
What are the arithmetic rules that govern Nodata values? Regardless of how these null values are defined? I am using  ArcMap Raster Calculator.

I am thinking specifically of operations like 1*Nodata (or 1*null), 1/Nodata, 1+Nodata,1-Nodata, Nodata-1 etc...


texts: Using PYQGIS add WMS layer from geoserver which is proxy protected

texts:
I have a wms layer hosted on a Geoserver used to display on a web map application using OpenLayers. I will get the map image using the proxy URL which sends a request to the Geoserver and returns the wms layer image as a response. But, I couldn't use the same URL in pyqgis to load the layer as the Raster layer in QGIS. I could load my localhost GeoServer data successfully using the QgsRasterLater function, but since my remote server's GeoServer could only be accessed via proxy URL, I am getting an invalid layer response.
Function to load localhost layer (Works successfully):
urlWithParams = "crs=EPSG:4326&format=image/png&layers=roads&styles=roads&url=http://localhost:7080/geoserver/demo_workspace/ows"
        rlayer = QgsRasterLayer(urlWithParams, 'roads', 'wms')
        print(rlayer.results())
        if not rlayer.isValid():
            print("Layer failed to load!")
        else:
            QgsProject.instance().addMapLayer(rlayer)

Function to load remote GeoServer layer which is proxy protected(Not working):
Case 1:
urlWithParams = "SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image/png&TRANSPARENT=true&LAYERS=demo_worksapce:layer1&CRS=EPSG:3857&TILED=true&WIDTH=256&HEIGHT=256&STYLES=&BBOX=8844681.416934315,1408887.3053523675,8922952.933898335,1487158.822316388&url=https://demo.com:2080/api/v1/proxy/image"
        rlayer = QgsRasterLayer(urlWithParams, 'roads', 'wms')
        print(rlayer.isValid())
        if not rlayer.isValid():
            print("Layer failed to load!")
        else:
           QgsProject.instance().addMapLayer(rlayer)

Case 2:
urlWithParams = "FORMAT=image/png&LAYERS=demo_worksapce:layer1&CRS=EPSG:3857&STYLES=&url=https://demo.com:2080/api/v1/proxy/image"
        rlayer = QgsRasterLayer(urlWithParams, 'roads', 'wms')
        print(rlayer.isValid())
        if not rlayer.isValid():
            print("Layer failed to load!")
        else:
           QgsProject.instance().addMapLayer(rlayer)

Both cases are failing. Kindly help me.


texts: I am using the Georeferencer plugin to overlay an older map over a OpenLayers google streetview map. I have set all my GCP points to be accurate however once I click the begin georeferencing button which should load my older map as a new vector layer only my red gcp points load and my map is no where to be found. Do you guys have any ideas as to how I can fix this? This is my first time posting here so if I left out important info feel free to ask for more information.


texts: Projection Issues - OpenLayers / MapServer / NAD83 State Plane FT

texts:
I am having projection issues with OpenLayers version 2.13 with MapServer and the Proj4 library.  Pulling layers in NAD83 State Plane Feet projections on top of an OSM or MapQuest Base Map in Google EPSG:900913 / EPSG:3857 is working perfectly in most locations.  However there are some states where an identical setup does not work, specifically Alaska and Vermont.

Alaska and Vermont are two cases where there is not an EPSG code for NAD83 State Plane Feet (only meters), but there are Esri codes in Feet. In Alaska I have tried ESRI:102631 (NAD 1983 StatePlane Alaska 1 FIPS 5001 Feet) and this works well in MapServer alone but the layers do not show in OpenLayers.

We have not upgraded the MapServer installation and Proj4 in a couple of years, so maybe a new version will help, but since things are working well in MapServer I think the issue is in OpenLayers. I have also tried creating a custom EPSG definition with no luck.

I am not sure what to try next for troubleshooting. Are there any issues in OpenLayers when using Esri codes? Should we be using NAD27 in these states? That would be OK except that we have a lot of imagery in NAD83 Feet that should probably not be converted to different a projection.

Here is the definition we are using in OpenLayers (NAD 1983 StatePlane Alaska 1 FIPS 5001 Feet):

Proj4js.defs["ESRI:102631"] = "+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=-36.86989764583333 +k=0.9999 +x_0=4999999.999999999 +y_0=-4999999.999999999 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs";


texts:
I am having some problems on downloading ASTER data.
When I download a Landsat-8 scene, it cames with the full bands on it (14 band files), an metadata file and a quality file.

I need to download ASTER data, with all ASTER bands on it.
I have tried to download ASTER_L1T_xxxx.hdf, but it only gives me 8 bands as per image below:

I have used Glovis (https://glovis.usgs.gov/) to download the image since it seems to be easier to use.


texts: How likely is it that Smartphone Camera GPS coordinates are wrong?

texts:
According to the coordinates printed directly on this photo of my late son Kelvin here:

https://www.dropbox.com/sc/w6tvwuvbosuh9rk/AADzT-3vDXcFQ9LnKK1BYMZUa

...and the Exif data from the .jpg:

GPSLatitudeRef: N
GPSLatitude: 24, 58, 34.6
GPSLongitudeRef: E
GPSLongitude: 121, 19, 37.768

...he was in Taiwan on Jan. 5, 2012 (two and a half months before he died).

However, I don't think he was ever in Taiwan, or anywhere near there. Is there any explanation for why a camera would generate the wrong coordinates? Maybe the whole thing is "messed up" because the supposed date of the image, as indicated by all the Exif data extractable:

GPSLatitudeRef: N
GPSLatitude: 24, 58, 34.6
GPSLongitudeRef: E
GPSLongitude: 121, 19, 37.768
Software: HTC Laputa Browser
DateTime: 2012:05:10 18:28:30

...indicates the picture was taken 5/10/2012, but he died 3/17/2012.

I can imagine the date of the camera being off as more likely than the calibration of the coordinates, though. The Bing map for "No. 23, Xinghua Rd, Taoyuan City Taoyuan County, Taiwan 330" (which is what I see when I enter "24° 58' 34.60" N, 121° 19' 37.77" E" into my browser).

The "software" entry from the Exif data says "HTC Laputa Browser", and the map does show a couple of "HTC"s on the map, whatever that is. I'm intrigued and confused. Can anybody shed any light on this?


texts: Exporting PostgreSQL table to FileGDB: wkbUnknown layer geometry type

texts:
I'm attempting to export a PostGIS (1.5.3) database table as a File geodatabase using the new FileGDB OGR / GDAL 1.9.0 driver:

$ ogr2ogr -overwrite -sql "SELECT * from db_pts LIMIT 1" \
    -f "FileGDB" /tmp/test.gdb \
    PG:"host='HOST' dbname='DB'"

I get this error:

ERROR 1: Error (-1): FGDB layers cannot be created with a wkbUnknown layer geometry type.
ERROR 1: Terminating translation prematurely after failed translation from sql statement.

The geometry column in the table I'm exported is set up properly and is in the geometry_columns table.  I think the problem has to do with the way the -sql option is returning the layer geometry type.  For example, accessing the same database table with ogrinfo:

$ ogrinfo -ro -sql "SELECT * from db_allobs LIMIT 1" PG:"host=HOST dbname='DB'"
Layer name: sql_statement
Geometry: Unknown (any)
Feature Count: 1
...
OGRFeature(sql_statement):0
...
  POINT (-151.472835 69.635245)

It's as though the query doesn't look in geometry_columns when building the layer, but later, when exporting the rows it knows that the geometry being exported is a POINT (in this case).

When I use a similar command to produce a shapefile, all is well, including the layer geometry, which is set to point.

Anyone have any ideas?


texts: Creating list of intersected polygons for each line using QGIS?

texts:
I have a line layer representing transport routes.

I would like to add a virtual layer listing every city crossed by each line. I would only need the id of each city to be listed (column Code_INSEE).

The result should be

Line 1 | 12001, 12003, 12122|
Line 2 | 23001, 23002, 23004, 23120,23012|
etc.

This question is linked to the solution proposed in this topic: Automatic calculation field based on line modification within QGIS


texts: Merging and Clustering Rasters using QGIS?

texts:
I have 3 rasters for three regionalized variables (r1, r2, r3). The rasters have the same projection, origin and cell spacing.

Can I merge this into ONE "multivariate" raster? (ie x,y,r1,r2,r3)

I would then like to then use a clustering/classification algorithm. My desired product is a 1 map with classified polygons.

I come from a SAS background and am new to QGIS.


texts:
I am trying to do a conflation. This is my method, and my problem:

METHOD

Step 1: I have two tables with road geometries. I am using ST_Dwithin to subselect geometries from table 'A' that lie within 500 meters of geometries from table 'B'. Let's call this selected part of table 'A', 'sel_a'. (This just to whittle the massive table 'A' into just useful geometries). I think I got this query right (waiting for results at the moment).

Step 2: Create points from geometries in table 'sel_a'; I am querying start, end and mid points for each geometry.

Step 3: Compute the closest line from table 'B' for each point obtained above. Then I can, for each geometry in table 'A' attribute a geometry from table 'B', such that the geometry from table 'B' was selected as being closest to at least two points created in step two. (Hit me in the comments if this is not clear.)

PROBLEM

The only way I can think of doing step 2 is using this query (which is pretty straightforward), but I'm not sure how to combine the results into two columns:

SELECT DISTINCT ON (sel_a.id)
     sel_a.id AS A_ID,
     ST_Line_Interpolate_Point(sel_a.geom, 0.0) AS sp,
     ST_Line_Interpolate_Point(sel_a.geom, 0.5) as mp,
     ST_Line_Interpolate_Point(sel_a.geom, 1.0) as ep,
   FROM sel_a

Which would result in a table that looks like this:

A_ID  sp          mp        ep
1     <start pt>  <mid pt>  <end pt>
2     <start pt>  <mid pt>  <end pt>
3     <start pt>  <mid pt>  <end pt>
...   ...         ...       ...

Now how do I then execute step 3? I will need to somehow transform the table that the above query returns into something that looks like this:

A_ID  Point_frm_A closest_B_ID
1     <start pt>  <to be computed>
1     <mid pt>    <to be computed>
1     <end pt>    <to be computed>
2     <start pt>  <to be computed>
2     <mid pt>    <to be computed>
2     <end pt>    <to be computed>
...   ...         ...

How do I do that? Or (preferably), can I modify the query above to directly result in the table structure I am looking for? (The latter will avoid an extra transposing step).


texts: ArcPy: 999999 ERROR when using AddJoin_management

texts:
I want to execute three joins in total. The first two joins work, but when I'm trying to join both resulting datasets I get the beloved 999999 ERROR.

# sr tables for joins
sr1Table = sr1Path + ".dbf"
sr2Table = sr2Path + ".dbf"
# create nw feature layer for join
nw1Layer = featureFolder + "nw1_layer"
nw2Layer = featureFolder + "nw2_layer"
arcpy.MakeFeatureLayer_management(nw1Path, nw1Layer)
arcpy.MakeFeatureLayer_management(nw2Path, nw2Layer)

# execute join on ID, sr1_layer will be joined to the nw1_layer
arcpy.SetProgressorLabel("Joining NW1 with SR1 ...")
arcpy.AddJoin_management(nw1Layer, "ID", sr1Table, "ID")
# execute join on ID, sr2_layer will be joined to the nw2_layer
arcpy.SetProgressorLabel("Joining NW2 with SR2 ...")
arcpy.AddJoin_management(nw2Layer, "ID", sr2Table, "ID")
# execute join on ID, nw2_layer will be joined to the nw1_layer
arcpy.SetProgressorLabel("Joining NW1 with NW2 ...")
arcpy.AddJoin_management(nw1Layer, "ID", nw2Layer, "ID")

Is this because I again have to create Feature Classes of nw1Layer + nw2Layer with CopyFeatures_management, then again use MakeFeatureLayer_management on the new Feature Classes and use these for the 3rd join? I try to avoid this, because I work with huge Shape Files and creating those Feature Classes would cost many hours of time.


texts: How to use CIMWFSServiceConnection in ArcGIS Pro

texts:
Does anyone know how to add a WFS service to a map using the CIMWFSServiceConnection class?
The WMS and WMTS equivalents work perfect, but when using CIMWFSServiceConnection, the layer that's added to the map has a red exclamation mark in front of it, and the Source entry on the layer's properties dialog shows Data source information not available
This is the code I used for adding a WMS service, the only differences are the URL, and the class used in line 2:
var serverConnection = new CIMInternetServerConnection { URL = "https://example.com/wms?service=WMS" };
connection = new CIMWMSServiceConnection { ServerConnection = serverConnection };
LayerFactory.Instance.CreateLayer(connection, MapView.Active.Map);

Adding a WFS through the user interface (i.e. Favorites > Add Item > New WFS Server) also works without problems.


texts: What are the possible listeners and event types for an openlayers map (ol.Map)?

texts: I cannot find this in the docs anywhere but it should be an easy answer for anyone who knows OpenLayers on the inside. In OpenLayers 4.2, what are the possible key values for listeners: during instantiation of an ol.Map object and what are the possible event values that can be passed to e.g. ol.Map.getView().on(event, function, opt_this) to trigger it?


texts: GDAL (Or Any Open Source GIS SDKS) | Read ESRI Geometric Network

texts:
I am looking at some legacy ArcObjects code that reads & parses a geometric network.  I no longer have access to ESRI license(s) at any level - so legacy code needs to be reworked since ArcObjects are off the table...

Are there any open source libraries that can read geometric networks?  Don't really care what language (python, c#, java,..etc).  My specific needs are:

Get Edge Ids of items (i.e. edges) particpating in geometric network
Correlate edge to feature in source feature class


texts: Maplex Street Placement exporting labels to Illustrator as separate letters?

texts:
I make city maps with ArcMap 10.2.2 that I export to Illustrator for finalization.
To label streets, I use the Maplex Street Placement as below:

As I need to be able to rework labels with Illustrator easily, I use the 'straight' placement to avoid curved labels.

The problem is even using the straight placement, in Illustrator side, I get some separated letters labels:

Most labels arrive in separated straight words, and some of them are texts on path. Both cases are convenient to edit, but separated letters case involves rewriting the whole street name.

Does anyone know how to prevent Maplex/Illustrator export from creating separated letters labels?


texts: Calculating proportional change in surface volume

texts:
I am using the Surface Volume tool (3D Analyst/Functional Surface) in ArcGIS 10.1 to calculate the volume between a reference plane and two water table surfaces (year 1 and year 2) that I interpolated from monitoring data (i.e., water table elevation points).

I specified that the tool should calculate the volume ABOVE the plane (up to the surface). My x,y, and z units are all in meters, so Z-factor = 1. I used the same reference plane height for both surfaces (it was set to a height that would not intersect either surface).

I am interested in calculating the proportional (percent) change in surface volume between years, not in the raw volume numbers per se. After running the tool, I did the following to calculate percent change: ((Year 2 vol - Year 1 vol)/ Year 1 vol)*100

I ran the Surface Volume tool multiple times using the same two surfaces (year 1 and year 2), but changed the height of the reference plane (the plane never intersected either surface). I expected the percent change to remain the same, regardless of the reference plane height. However, the percent change in volume between the two years changes each time I alter the reference plane height.

Here is an example of my data:

Percent change in volume between year 1 and 2:

-0.31% (ref plane at 725 m)
-2.10% (ref plane at 760 m)
-5.46% (ref plane at 763.7 m)

As I raise the plane, it gets closer to my surfaces, and the raw volume numbers decline.

Why does the percent change differ as I change the reference plane height?

Perhaps it is related to the overall downhill trend in my water table data.


texts: How to create a solar intensity map for a garden plot?

texts:
I'm interested in creating a raster surface that depicts the average solar intensity per unit area for a small garden plot.  I imagine the inputs into the model would be the following:

Plot area--for simplicity sake, let's assume the garden plot is
circular.
Average tree height of surrounding forest edge
Solar intensity at various sun angles
The average sun angle in relation to the garden plot over time.

I ultimately want to show the best location for plant growth within garden plots that are surrounded by trees.  Could anyone point me toward a solution for this analysis?  I'm sure someone has done this already, although I cannot find any references.

Edit:

For example purposes, the attached image shows a 1 acre garden plot surrounded by an even-aged stand of conifers.  To further simplify this analysis, let's assume we're interested in only the solar zenith and intensity during the gardening season from June to October.  Felling trees is not an option!


texts:
I'm writing an Android navigation application using Java in my own time as a learning exercise. It reads Shapefiles (with some pre-processing by me) downloaded from http://download.geofabrik.de/europe.html.

I can't find elevation/topographical data in Shapefile format to download anywhere.

Can anyone suggest a free source?

Alternatively, I could invest time in learning another GIS file format and writing code to read it. I'd like to choose a format that can be read quickly by devices with low processing power, and is also likely to be widely used in future.

I want to display elevation data in some form in the navigation application I am writing. I'd like to try various approaches, so data in the form of spot heights or contour lines would be useful.


texts: Import 3D textured object into ArcGIS Pro

texts:
I want to import a 3D mesh which is exported from Photoscan. I tried with Collada export (.dae) and 3ds models (.3ds).

I created a global scene in ArcGIS and imported the file with the function "import 3D objects", and defined a file database in order to put the multipatch features in.

Now, the imported object has a solid color (see printscreen below) but I want to add the texture to the 3D object (which was exported in JPEG from Photoscan). The function "replace the multipatch" does not work, as it used to with ArcScene.


texts:
I'm working on a project where we have a large number of very detailed polygons being displayed on a image using Openlayers.  The polygons are added to the map through a GML feed.

We are having end user performance issues on some of our more complicated polygon feeds.  The issues are two fold:

1) The browser will freeze up when loading the polygons from GML initially.  After the loading it seems to perform pretty well is there some way to stream the data in slower as to not freeze the browser when the polygons are initially being added to the map?

2) The data payload to the user is too large.  We would like to simplify the actual amount of polygon data that is sent to the user.  We have our source polygon data in a mysql database, is there an algorithm we could use that would simplify our source polygon data to reduce the number of vertexes?

Some of this may be very basic, sorry I'm pretty new to this stuff.


texts:
I'm trying to send a JSON file to my Node.js app via Protocol.HTTP. The file is correctly sent, but I can't access to an array. Here is the relevant code :

Client side

(...)
    var vectorProtocol = new OpenLayers.Protocol.HTTP({
                    url: '/getcoords',
                    format: vectorFormat,
                    readWithPOST: true,
                    params: {
                        "code_company": "8501",
                        "data_company": [
                            {
                                "origin": "2013P00109",
                                "type": "LJ",
                                "naf": "5610A",
                            },
                            {
                                "origin": "2013P00110",
                                "type": "FJ",
                                "naf": "5481"
                            }
                        ]
                    }
                });
(...)

Server side, I try to build an array with only the "origin" field of my array "data_company":

function getCoords(params, callback) {
      var arrOrigin = params.data_company.map(function(d) {
         return d.origin;
      });
      (...)
}

And I get this error :

  TypeError: Object [object Object] has no method "map"

It seems that my "data_company" array is not parsed. Is there a reason for that ? And a workaround ?

Thanks


texts: GeoPandas and Cartopy how to dynamically zoom to Multiline

texts:
I am new to GeoPandas and Cartopy.
I am trying to generate map plots for different regions dynamically.
When generating the image from geopandas ('df.plot()') the size is correct, but when I generate a plot using Cartopy the Multiline appears very small and I cannot find a way to fit it to the figure. What I am doing wrong?
import psycopg2 as db
import matplotlib.pyplot as plt
import geopandas
from cartopy import crs as ccrs

con = db.connect(connection_string)
sql = 'select geom from table'
df = geopandas.GeoDataFrame.from_postgis(sql, con,crs= 21781)
df.plot(color='black')

crs_epsg = ccrs.epsg('21781')
df_epsg = df.to_crs(epsg='21781')
fig, axs = plt.subplots(1, 1, subplot_kw={'projection': crs_epsg},figsize = (5,2.5))
axs.add_geometries(df_epsg['geom'], crs=crs_epsg,facecolor="none", edgecolor='black')
extend = [650000,690000,22000,260000]
axs.set_extent(extend,crs=crs_epsg)
axs.outline_patch.set_visible(False)
axs.background_patch.set_visible(False)
axs.background_patch.set_alpha(0)
fig.patch.set_alpha(0)

plt.show()

Here the image of GeoPandas directly:

And this is the image made with GeoPandas and Cartopy:

Here one of the Multistrings: Sample MultiString


texts:
I am trying to analyse time-series images from MODIS (MOD13Q1) in TIMESAT to clear the noise and get an output.
With EVI images I got also QA. I have to use it because there is often a cloud cover in my working area. I saw a few answers how to weight a certain range of values but I am not sure about it, especially because there are only 3 rows in TIMESAT. Can you give me an example of how to weight a MODIS data?
After reading many manuals and articles I still feel uncertain about my work and I have a couple of questions.


texts: QGIS custom form : mixing Qt Widgets and QtQuick QML Widgets

texts:
On QGIS 3.16 I am trying to create a custom form with some Python logic for feature navigation and a QML widget to display nicely images.
I embedded my QtQuick QML file in the custom form python file with this:
def add_qml_widget_to_form():
    """Add QML widget in form."""
    global myDialog  # This is set in the main form function
    # manually drawn vertical layout in the corresponding .ui file
    vbox = myDialog.findChild(QVBoxLayout,"photo_layout")
    view = QQuickView()
    view.setResizeMode(QQuickView.SizeRootObjectToView)
    qml_url = "C:\\Path\\to\\files\\CustomForm.qml"
    path = QDir.fromNativeSeparators(qml_url)
    view.setSource(QUrl.fromLocalFile(path))
    if view.status() == QQuickView.Error:
        for error in view.errors():
            QgsMessageLog.logMessage(error.description())
    else:
        container = QWidget.createWindowContainer(view, parent=myDialog)
        container.show()
        vbox.addWidget(container)

The QML file uses a Repeater which fetches data from a table Relation:
 Repeater {
    model: expression.evaluate("relation_aggregate( 'test3','array_agg', coalesce(\"url\",'')||'ÿ'||coalesce(\"photo_type\",'')||'ÿ'||coalesce(\"orientation\",''))")

Unfortunately, the expression doesn't evaluate correctly (no output), and I guess it is missing the form scope with layer and feature information to work.
I would like my code to work the same as when I use the QML Widget in a drag-and-drop form which is feature-aware and evaluate correctly the expression.
Any ideas how I could deal with that ?
Another option would be to switch entirely:

to QML, but I have seen little documentation on how to use PyQt functions to connect to signals for my navigation needs.
to QtWidgets, but everyone seems so extatic on QML it doesn't seem the right direction.


texts: Applying symbology to group of layer using ArcPy in ArcGis 10.1?

texts:
I am creating a python script and one of the things i need to do is apply symbology from layer to a group of shapefile in my mxd . But when i specified the folder of symbology layer, symbology deosn't working, without any error message.

import arcpy.mapping
mxd=arcpy.mapping.MapDocument("CURRENT")
def importer ():
    # Set outWorkspace
    outWorkspace = r'c:\Users\HAIDRA\Documents\Desktop\essai'
    mxd = arcpy.mapping.MapDocument("current")
    style = askdirectory ()
    styl_path.set(style)
    # set the workspace
    arcpy.env.workspace = style
    list = arcpy.ListFiles("*")
    print list
    # Execute MakeFeatureLayer for each input shapefile
    for li in list:
        #make a layer from a feature
        styllist = arcpy.MakeFeatureLayer_management (li, li +"_lyr")
        for styl in styllist:
            print styl
    # Script arguments
    styllist = styl_layers
    styl_layers = arcpy.GetParameterAsText()
    layer = Layers_to_Symbolize
    Layers_to_Symbolize = arcpy.GetParameterAsText()
    # Process: Apply Symbology From Layer
    for UpdateLayer in Layers_to_Symbolize:
            arcpy.AddMessage("Updating: " + UpdateLayer)
            arcpy.ApplySymbologyFromLayer_management(UpdateLayer,styl_layer)
    Output_Symbology = styl_layer
    arcpy.RefreshTOC()
    del mxd


texts: Postgres: How do I get the closest distance for each element within a point coordinate table from a polyline table? And if this is possible what version of Postgres is this doable.


texts: How to use `layer.getFeatures()` in a standalone script?

texts:
I am following the QGIS Cookbook and this post where it concerns reading attributes from layers. I would like to exploit this method and implement it into a standalone script.

Essentially, I want to read the first feature of a shapefile from a field called Rank, take the value of this feature (all values in the Rank field are the exact same) and include it into a QGIS algorithm.

I wonder if there's an alternative to using layer.getFeatures() as, I believe, this is used with the canvas when QGIS is running. I have the following basic script but doesn't like fname.getFeatures():

path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

for fname in glob.glob(path_dir + "*.shp"):
    for feature in fname.getFeatures():
        fid = 1
        iterator = fname.getFeatures(QgsFeatureRequest().setFilterFid(fid))
        idx = fname.fieldNameIndex('Rank')
        Rank_value = idx[1]

        #Do some processing, finish with field calculator algorithm

        general.runalg("qgis:fieldcalculator", fname, 'Rank', 1, 10, 0, False, Rank_value, path_res + fname)

Is this the correct method or would it be much easier to read the Rank_value from a shapefile and write to another programatically?

EDIT:

Following @gcarrillo's answer, the script does not produce any results. I have modified the script to match the one suggested but not sure why no output files are produced.

To try and clarify exactly what I'm looking for:

A value from the Rank field of a shapefile is read.
A grid is created.
The same shapefile is clipped onto the grid (it would then have more polygons, thus more features).
The Rank field for this newly clipped shapefile is updated with the value read at the start.
Next shapefile...

path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

for fname in glob.glob(path_dir + "*.shp"):
    vLayer = QgsVectorLayer(fname,"any name","ogr")
    idx = fname.fieldNameIndex('Rank')
    Rank_value = vLayer.uniqueValues(idx)[0]

    outputs_1=general.runalg("qgis:creategrid", 1000, 1000, 24108, 18351.157175, 258293.802316, 665638.226408, 1, 'EPSG:7405', None)
    outputs_2=general.runalg("qgis:clip", outputs_1['SAVENAME'], fname, None)
    outputs_3=general.runalg("qgis:fieldcalculator", outputs_2, 'Rank', 1, 10, 0, False, Rank_value, path_res + fname)


texts:
As presented in a talk at FOSS4G Mapbox Studio allows to create Mapbox vector tiles and export them as a .mbtiles file.

The mapbox-gl.js library can be used to dynamically style and render Mapbox vector tiles on client (browser) side.

The missing part: How can I self-host Mapbox vector tiles (.mbtiles) so that I can consume them with mapbox-gl.js?

I know that Mapbox Studio can upload the vector tiles to the Mapbox server and let it host the tiles. But that's no option for me, I want to host the vector tiles on my own server.

The TileStream approach below turned out to be a dead end. See my answer for a working solution with Tilelive.

I tried TileStream which can serve image tiles out of .mbtiles files:

My webpage uses mapbox-gl v0.4.0:

<link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.4.0/mapbox-gl.css' rel='stylesheet' />
<script src='https://api.tiles.mapbox.com/mapbox-gl-js/v0.4.0/mapbox-gl.js'></script>

and it creates a mapboxgl.Map in a JavaScript script:

  var map = new mapboxgl.Map({
    container: 'map',
    center: [46.8104, 8.2452],
    zoom: 9,
    style: 'c.json'
  });

The c.json style file configures the vector tile source:

{
  "version": 6,
  "sprite": "https://www.mapbox.com/mapbox-gl-styles/sprites/bright",
  "glyphs": "mapbox://fontstack/{fontstack}/{range}.pbf",
  "constants": {
    "@land": "#808080",
    "@earth": "#805040",
    "@water": "#a0c8f0",
    "@road": "#000000"
  },
  "sources": {
    "osm_roads": {
      "type": "vector",
      "url": "tile.json"
    }
  },
  "layers": [{
    "id": "background",
    "type": "background",
    "paint": {
      "background-color": "@land"
    }
  }, {
    "id": "roads",
    "type": "line",
    "source": "osm_roads",
    "source-layer": "roads",
    "paint": {
      "line-color": "@road"
    }
  }]
}

... with the following TileJSON specification in tile.json:

{
  "tilejson": "2.1.0",
  "tiles": [
    "http://localhost:8888/v2/osm_roads/{z}/{x}/{y}.png"
  ],
  "minzoom": 0,
  "maxzoom": 12
}

... which points to my TileStream server running at localhost:8888.
TileStream has been started with:

node index.js start --tiles="..\tiles"

... where the ..\tiles folder contains my osm_roads.mbtiles file.

With this setup, I can open my webpage but only see the background layer. In the browser network trace I can see that tiles are indeed loaded when I zoom in, but the browser JavaScript error console contains several errors of the form

Error: Invalid UTF-8 codepoint: 160      in mapbox-gl.js:7

Since vector tiles are not .png images but rather ProtoBuf files, the tiles URL http://localhost:8888/v2/osm_roads/{z}/{x}/{y}.pbf would actually make more sense, but that doesn't work.

Any ideas?


texts: Inconsistent Spatialite Update Behavior; SRID not allowed?

texts:
I have two tables: crime and (police) stations. The latter has only 4 records with (name, lat, lng) fields. The former has over 500K records with two California Zone 2 (ft) coordinate columns named X_Coord and Y_Coord, respectively. I use R to populate the SQLite database and am building a script to turn it spatial and prepare it for analysis.

My problem is that I can populate the crime geometry, but not the station geometry. I'm using the stand-alone spatialite executable. I run the init_spatialite.sql script. Then I run the four SQL commands below, with the last one failing: "SQL error: 'stations.geom' violates Geometry constraint [SRID not allowed]". I ran a boolean and count check to make sure the crime geometry matched the coordinates (they do). So for whatever reason, the same approach doesn't work on the stations table. I don't understand why I'm getting the constraint violation (the SRID exists in the metatable and is correct).

This may be neither here nor there, but the MakePoint function doesn't exist, even though it's in the documentation (and saw it as another possible approach to take).

SELECT AddGeometryColumn('crime',    'geom', 2226, 'POINT', 2);
SELECT AddGeometryColumn('stations', 'geom', 4326, 'POINT', 2);
UPDATE crime    SET geom = GeomFromText('POINT(' || X_Coord || ' ' || Y_Coord || ')', 2226);
UPDATE stations SET geom = GeomFromText('POINT(' || lng     || ' ' || lat     || ')', 4326);

EDIT: I followed the same routine in the spatialite-gui and it worked just fine. This is odd behavior, to say the least. I'll try this from home where I have (sudo) admin right, for today I've been practicing at work where I am a limited user.


texts:
We have a new install of ArcGIS Portal running with the WebAdaptor using IIS.  If we login directly, it works just fine.  However, we have a mandate to use a CAS login process, and when we enable that, we end up with a redirect loop upon login.

What definitely happens is that, when you go to the main page, it does redirect to the CAS server, and then redirects back to the home page.  However, it seems to not find the cookie with the CAS-generated token there, and redirects back to the signin page.  The signin page finds it in the session, I guess (at least, it definitely doesn't go back to the CAS server), and redirects back to the home page, and the loop continues until the browser calls a stop to it.

The CAS client is the DotNetCasClient, and its configuration in the Web.config for the WebAdaptor app has:

<configuration>
  <configSections>
    <section name="casClientConfig" type="DotNetCasClient.Configuration.CasClientConfiguration, DotNotCasClient" />
  </configSections>
  <casClientConfig casServerLoginUrl="https://loginHost/cas/login"
               casServerUrlPrefix="https://loginHost/cas"
               serverName="https://portalServer"
               redirectAfterValidation="true"
               gateway="false"
               renew="false"
               singleSignOut="false"
               ticketValidatorName="Cas10" />
...
  <system.web>
...
    <authentication mode="Forms">
       <forms cookieless="UseCookies"
          domain="our-domain" <!-- CAS and portal are both in same domain -->
          loginUrl="https://loginServer/cas/login" <!-- same url as above -->
          name="CasLoginTicket"
          slidingExpiration="true"
          timeout="30" />
    </authentication>

...

We've tried tweaking most of these parameters.  Changing redirectAfterValidation makes it never return from the CAS server.  Changing singleSignOut to true and adding
    serviceTicketManager="CacheServiceTicketManager"
changed nothing.

Changing the ticket cookie name to "esri_auth" got rid of the loop, but resulted in no login.

And I checked on browser cookie settings; we're accepting 1st- and 3rd-party cookies, as well as session cookies.

Any ideas?

I don't have all the details, since I didn't implement it, but we were able to talk with the CAS server via SAML.  That took care of the problem.


texts:
We are building an app starting from the ArcGIS Tax Parcel Viewer.  I have added buttons to allow users to switch to Street View or Pictometry when they click on the map.  All that is working fine, but now I want to change to a custom cursor after pressing a button so they know when these options have been turned on.

Using CSS like this:

cursor: url(../images/cursors/Pictometry.cur), default;

I have been able to turn on the custom cursors over pretty much every DIV in the page except the map.  Even when I enable it for the entire <html> or <body> tags, as soon as the cursor goes over the map, it reverts to the default look.  So far, I have been unable to locate any CSS or HTML elements that I can change to enable this (either in the code or programmatically).  Does anyone know if this can even be changed, and if so, where?


texts: problem overlaying osm with os-derived data in QGIS - both British National Grid

texts:
I am currently working alongside city council datasets derived from Ordnance Survey; this is projected in their ArcSDE system as "British National Grid" or EPSG:27700.  My problem is that I am trying to overlay data downloaded from OpenStreetMap - which I do using osm2pgsql to export into a new PostgreSQL db creating PostGIS layers with the projection system ESPG:27700.

However, when I try to overlay these, the two data sources are offset from each other by 50m.  The city council state that the OSM data is wrong, siting that it does not overlay well onto their aerial photography; whereas if I add the datasets to something like MapBox Streets, it would appear that the Ordnance Survey data is the one that is offset.

Which data source is wrong?  How does such an offset arise?  What methods can be used to align & project them correctly?

Three questions that have been burning in my mind for the last few weeks; any advice or insight on the problem would be most appreciated.

Yours in geography,

Paul


texts: Creating minimum bounding geometry using the convex hull method in R

texts:
I would like to create a small enveloppe around a set of polygons that encloses around them. I would like to avoid using the bounding box function (example below) because it would be too big. Rather, I would like to use a convex-hull like algorithm to create a polygon just around the desired polygons. At the very least, an oblique rectangle/diamond would be an improvement to cut out all the corners.

In QGIS, the function is Minimum bounding geometry with the geometry type Convex hull.

Sample code for image below

data(meuse)
coordinates(meuse) = ~x+y
meuse<-as(meuse, "sf")
meuse_poly<-st_buffer(meuse, dist=meuse$elev*15)
meuse_poly$ID<-c(1:nrow(meuse_poly))

plot(st_geometry(meuse_poly))
plot(st_geometry(meuse_poly[which(meuse_poly$landuse=="Bw"),]), add=T, col="red")
plot(st_geometry(st_as_sfc(st_bbox(meuse_poly[which(meuse_poly$landuse=="Bw"),]))), add=T, border="blue")


texts:
As far as I know neither PROJ4 nor ESRI tools can apply the Peirce quincuncial projection.

Does anybody know what libraries/softwares can manage it?


texts: Selecting more than 1,000 features from AGOL hosted feature service?

texts:
When data from a cloud-hosted feature service is loaded into ArcMap, it queries the first 1,000 features "to enhance performance"... I, however, need to see all the features.  How can I add all features to my map? I want to do a selection on the service directly according to this blog post:
https://blogs.esri.com/esri/arcgis/2013/10/10/quick-tips-consuming-feature-services-with-geoprocessing/

  Max Records

  If you’ve already performed some analysis on a hosted feature service
  from arcgis.com, you may have encountered the following warning:

  WARNING 001406: Query operation has exceeded the feature service
  transfer limit

  Only the first 1000 records will be used by the geoprocessing tools.
  Hosted feature services have their maximum query set to 1000 records.
  1000 is also the default maximum with a local ArcGIS for Server.
  However you have the option to increase this value in your local
  deployment. To work with a feature service that has more than 1000
  records you will have to iterate through the records 1000 (or less) at
  a time. Remember that geoprocessing tools honor selection sets. You
  can manually select 1000 records at a time, or you can use a
  geoprocessing tool like Select Layer by Attribute. Tools like Merge
  and Append can be used to combine all your features into a single
  feature class depending on your needs.

However, when I go to use the tool, Select Layer By Attribute, I am unable to use the hosted feature service directly as my Input layer. Is there a specific way I have to do this? Has anyone created a custom tool or script to do this?


texts: How to show multiple features selection button/dialog box on the map using .NET?

texts:
If you use the Edit tool to click a location on the map with overlapping features, a small icon  Selection chip pops up with a dialog box previewing the selectable features where you clicked, so you can distinguish the correct feature from the overlapping ones.

How to get this "selection chip" from code (using .NET)? What interface I should use?


texts: Splitting shapefile into features and keeping attributes using ArcPy?

texts:
I want to split a shapefile so each feature becomes its own shapefile. The below works:

with arcpy.da.SearchCursor(catch10, ["FID", "SHAPE@"]) as cursor:
     for row in cursor:
      out_name = "Sup_" + str(row1[0]) #
      arcpy.FeatureClassToFeatureClass_conversion(row1[1],outSup, out_name, **field_mapping=**)

But the attributes don't come with it. It gives back a shapefile with FID (0), Shape, ID (0). I have tried different entries under field mapping -- long strings copied from results and shorter {"Name" etc) -- but nothing has changed it. I'm missing something


texts: Coordinate Reference System in QGIS: Project vs Layer

texts:
I do not understand why there can be in QGIS different CRS for project and for layers. I have worked previously with subsurface petroleum software and there CRS "stuff" was lot more strict in that sense that data was almost always tied to the CRS of the project.
Will the data be prone to inconsistency if there is different CRS for project and layers?


texts: How to download all OSM data within a boundingbox with Overpass?

texts:
Is it possible to download openstreetmap by direct link through http://overpass-turbo.eu without removing any data?

Basically I am trying to download the openstreetmap data with defined boundaries in the same way that it is currently possible through openstreetmap.org:

http://api.openstreetmap.org/api/0.6/map?bbox=11,50,11.1,50.1

It's just the output should be .geojson instead of of .osm.

I am not sure I understand the Overpass API syntax.
I need something similar to this:

[out:json];
node(11,50,11.1,50.1)
out body;

But how should I make the download link from it (assumed the upper syntax is correct, which it isn't)?

https://overpass-api.de/api/interpreter?data= [out:json](48,11,48.19,11.73)... ?


texts: How to reclassify grey scale tiff in qgis

texts: I downloaded a raster image of the sea floor age. The values in band 1 ranges from -32768 to 32767. I need to reclassify it to no. of million of years to make sense of the data. How can I do so? I try to use Raster Calculator but I do not know how to input the syntax. Any kind soul able to help? Thanks.


texts:
I am following this tutorial http://workshops.boundlessgeo.com/tutorial-wordmap/ on the web about using setting up sql view on geoserver. I have set up my sql view correctly and have the javasctipt part working. When I enter a location in the textbox and press enter, the map will only display features of that particular location. However, I would like to go further so that my map can also zoom into the features of that location.

This is my code:

// Start with map of startWord
wms11.mergeNewParams({viewparams: "word:"+startWord});

//Add all layers to the map
map.addLayers([osm,wms11]);

// Text field component. On 'enter' update the WMS URL
var textField = new Ext.form.TextField({
    value: startWord,
    listeners: {
    specialkey: function(fld, e) {
        // Only update the map when the user hits 'enter'
        if (e.getKey() == e.ENTER) {
            wms11.mergeNewParams({viewparams: "word:"+fld.getValue()});
            var bounds = wms11.getDataExtent();
            alert(bounds);
            map.zoomToExtent(bounds);});
            }
        }
    }
});

The map was zoomed into an incorrect area and the alert(bounds) shows the bounds is undefined. What is going wrong here?

Thanks in advance!


texts:
Im trying to load a tile service in my QGIS 2.8 map canvas.
Talking about the following service: https://tiles4.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Dijkenkaart_van_Nederland/MapServer

Ive tried several methods but nothing seems to work.
As you can see by clicking the link, there are json, api and tile urls.

There must be a way to make this service useful?


texts: How to obtain elevation differences between point and raster cell in R

texts:
I have an empty raster and point shapefile (points = black dots in the image). For the analysis I use a Digital Elevation Model (DEM) with the same extent and cell size as the empty raster.

I would like to fill the empty raster cells with the elevation difference between every raster cell of the empty raster and the raster cell where the closest point is situated.

The green lines show distance between every raster cell the respective closest point.

Result: a raster with every cell filled difference in elevation between the cell elevation and the elevation from the cell where the point is located.

The first thing to know is: Which is the closest point?
The distance between a raster cell and the closest point can be calculated with the distanceFromPoint tool:

r <- raster('./Raster/DEM/DEM_Refocus_5kmBuffer.tif')
s <- readOGR('./Shapefiles und CSV/River/Rive_Punkte_3035.shp')

dfp <- distanceFromPoints(r, s)

But what I actually need is the elevation differnce between every raster cell of the raster and the the points. Do you have an idea to obtain this?


texts: Deny steps with count more than 3 in osm2po

texts:
wtr.tag.highway.steps =     1, 1, 1,   foot

In osm2po I can deny stairs with count 3 steps like this:

wtr.deny.step_count.3 = foot

But, I would like to prohibit stairs with a number of stages over 3 (so 4 and more).

wtr.deny.step_count.>3 = foot ---> that does not work in osm2po

wtr.allow.step_count.[1|2|3] = foot ---> that does not work in osm2po


texts: SLD layer style colour ranging from x -> y

texts:
I am new to web mapping. I would like to create a SLD style that colours the points if they fall between certain numbers for example 0-500. I know I could manually create each colour for each value and assign it. However this is not practical when you reach high numbers of values that need colouring. One option is to group them into small groups of 10, but that still gives me 50 colours to do.

I have a few hundred points, served from GeoServer via WMS, that rance from 0-500.

How would I go about doing this with SLD? The image below shows an example of what I would expect in a desktop client.


texts: Using expression operators in if/else conditions without getInfo in Google Earth Engine

texts:
İ have a two-segment list with id and value and i want to make a new list from imageList using second element value.
i want to get the first values in first column until 100 in second colum and save the first three values in first column to another list.
var imageList = ee.List([
  [101, 102],
  [102, 103],
  [103, 104],
  [106, 100],
  [104, 105],
  [105, 106],
  ])

i define starting point with firstImageId :
var firstImageId = 101

i get list length with getInfo:
var size = imageList.length().getInfo()

And running the function:
var imageList = getImageList(firstImageId,imageList)
print('New image List :', imageList)

And there are main function and helped function finding for indexof key:
// functions
function getImageList(firstImage, list) {
var cropImageList =ee.List([]);
var currentImageId= firstImageId

if(firstImageId) {
for (var i = 0; i<size;i++){
    print(i, currentImageId)
    var indexofkey = indexOf(list, currentImageId)
    var currentImageValueID= ee.List(list.get(indexofkey)).get(1)
    currentImageValueID= currentImageValueID.getInfo();
    if (currentImageValueID ==100 ) {
      print ('Break at : ',currentImageValueID)
      break;
    }
    cropImageList = cropImageList.add(currentImageValueID)
    currentImageId= currentImageValueID
    }
   }

  return cropImageList
}
// Thanks to Daniel Wiell for indexOf function
function indexOf(list2D, value) {
var array = ee.Array(list2D)
var keys = array.slice(1, 0, 1).project([0]).toList()
var values = array.slice(1, 1, 2).project([0]).toList()
return keys.indexOf(value)
}

The running code link is below:
https://code.earthengine.google.com/565d3ce033edcb427597a68dbea80889

The code is working but running slow because i use getInfo: I must use an expression operator to break for loop.
My question is, how can i define if/else conditions with local variables without using getInfo for below section:
        currentImageValueID= currentImageValueID.getInfo();
        if (currentImageValueID ==100 ) {
          print ('Break at : ',currentImageValueID)
          break;
        }


texts: Export a vector layer in QGIS 3 processing script only as GPKG

texts:
I'would export from a processing script (QGIS 3) a vector layer as GPKG. I'm using self.addParameter(QgsProcessingParameterFileDestination(...)) to define the destination file. This method allows user to select the proper extension file (SHP, GPKG, ...) but I'd export the output vector of my script just in GPKG because I'm using QgsVectorFileWriter to create the output, which admits just one type of extension per time.
Or I'would save the QgsVectorLayer object that I created in the script with any extension selected by the user.


texts: Join by location on Qgis when target layer has multiple layer

texts:
I posted here a while ago about spatial joins. I have a huge dataset (20,000?) of bird and mammal ranges.
I have another dataset of all the countries.

In qgis, I was attempting to link the two shapefiles so that I could make a dataset containing all the countries the animals live in.

I thought I had it sorted with 'join by location' on qgis, but I realised if the animals range was in more than one country, it only saved the first country.

Can anyone tell me how I could do this differently so that I can get all the countries, not just the first?

Edited:both datasets are shapefiles.


texts: Watching folder for arrival of files to batch process?

texts:
Is there a way to maintain a watch on a folder to enable automatic file type conversion? I have a program that dumps DXF files - I batch convert all files to both shapefile and mapinfo format.

It would be great to automate it further so that it automatically undertakes the conversions.


texts:
Question: Is it possible to automate the functionality of right-click layer in TOC > data > export raster using arcpy in ArcGIS Pro? Specifically including the clipping extent options? I have reviewed several similar SE posts (most relevant here), but all of them suggest using "copy raster", which does not have the aforementioned clipping options.

Background: I'm trying to download data from an image service (located here). I only want to download subset of the data, similar to a clip or extract by mask (see image below)

If I right-click and do data > export raster with the full AOI, i get the following error.

If I tile up the AOI into much smaller pieces it works perfectly. However, I then face the issues of having to manually right click and extract hundreds of times. I would instead like to automate this with arcpy, but 'export raster'. doesn't seem to be an actual geoprocessing tool.

edit: I've also tried using extract by mask, which works but is far slower. ~5 minutes per 'tile', as opposed to just a few second with when using 'Export Raster'


texts:
I'm trying to follow this example on Mapbox blog, about using 3D shapes in PostGIS.

The database-query is this:

DROP TABLE tilemill_buildings;
CREATE TABLE tilemill_buildings
AS SELECT
ogc_fid,
CAST(size_sqm AS INT) as size_sqm,
CAST(ROUND(CAST(ST_ZMax(wkb_geometry)/2 AS NUMERIC),0) AS INT) AS z_xlarge,
CAST(ROUND(CAST(ST_ZMax(wkb_geometry)/2.5 AS NUMERIC),0) AS INT) AS z_large,
CAST(ROUND(CAST(ST_ZMax(wkb_geometry)/3 AS NUMERIC),0) AS INT) AS z_med,
CAST(ROUND(CAST(ST_ZMax(wkb_geometry)/4 AS NUMERIC),0) AS INT) AS z_small,
wkb_geometry AS wkb_geometry
FROM building_3d
ORDER BY
ST_YMax(Box3D(wkb_geometry)) DESC;

...but I get this error:

ERROR:  function st_zmax(bytea) does not exist
LINE 5: CAST(ROUND(CAST(ST_ZMax(wkb_geometry)/2 AS NUMERIC),0) AS IN...
                        ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

I have tried to re-install PostGIS twice.  I have run the CREATE EXTENSION postgis; command.

Why am I getting this error, and how to fix it?


texts: How to map using more 'accurate' DATUM than WGS84?

texts:
WGS 84 is not very accurate as it is a compromise for the whole earth.

Assuming theoretically I'm surveying using very accurate IMU and very accurate GPS (in WGS 84 of course). I think I'll get a contradicting results:

While the GPS (as accurate as it is) will always throw me back to that 'compromise for the whole earth', the IMU on the other hand will very accurately draw a relative local 'true' map.

Am I correct here, that in such cases I need to convert the GPS's WGS84 coordinates into some more accurate DATUM dedicated for the region I'm in (ITRF?), in order to align IMU's relative map to some global frame? Because the global frame of that more proper DATUM to my region is closer to the local map the IMU's outputed rather than the rough estimation of the regular WGS 84?

EDIT:

The black line is the result of the IMU.

The blue line is WGS84

The green line is GDA2020 (assuming it's Australia).

coord 1 blue != coord 1 green

coord 1 blue is what I get from the GPS

coord 1 green since it's GDA2020 is more acculturate represents reality (the black line).

Am I correct here?


texts: Does the location of a DEM affect hillshade outputs?

texts:
I have a generic digital surface model (DSM) and I want to simulate illumination in different parts of the country on the same day, all using the same DSM. (i.e., how does illumination differ in Seattle and San Antonio on March 4?) I am using the ArcMap v10.4 hillshade function, but there are other implementations of hillshade functions in various software packages.

The ArcMap hillshade function requires the solar elevation and azimuth as inputs. These can be calculated for a given location and date/time using NOAA's solar position calculator.

I want to know if the spatial location of the DSM affects the output of the hillshade function, or if the elevation and azimuth are the only variables that affect the output.

To test this, I created three separate rasters as inputs: my original projected DSM, the DSM with the projection removed, and the DSM projected but shifted an arbitrary number of X and Y units. In R:

require(raster)

#Input one
DSM <- raster("DSM.tif")

#Input two
DSM_nocrs <- raster("DSM.tif")
crs(DSM_nocrs) <- NULL

#Input three
DSM_shifted <- shift(DSM, x = 100000, y = 100000)

I inputted these three rasters into the ArcMap hillshade function using the same solar elevation and azimuth. I saved the outputs to my working directory. I then compared the raw values of each output in R.

DSM_hillshade <- raster("DSM_hillshade.tif")
DSM_nocrs_hillshade <- raster("DSM_nocrs_hillshade.tif")
DMS_shift_hillshade <- raster("DMS_shift_hillshade.tif")

all.equal(values(DSM_hillshade), values(DSM_nocrs_hillshade))
[1] TRUE

all.equal(values(DSM_hillshade), values(DSM_shift_hillshade))
[1] TRUE

As you can see, R indicates that the values of the three hillshade rasters are the same even though one wasn't projected and one had a different spatial location. This suggests that the spatial location of the DEM/DSM has no impact on the output of the hillshade function. Perhaps then, for my analysis, I do not need to change the spatial location of my generic DSM and only need to change the solar elevation and azimuth.

My question: Does this seem like a fair interpretation? Is there any reason to believe that the spatial location of a projected DEM/DSM has an impact of the output of a hillshade function?


texts: Requesting styled Google Map tiles via URL

texts:
I'm not working with the Google Maps API but I'd still like to be able to request styled map tiles. Is it possible to request a styled Google Map tile via a URL similar to the slippy map URL format for OSM, eg:
http://[abc].tile.cloudmade.com/your_CloudMade_API_key/2/256/zoom/x/y.png

The '2' in the URL above represents the map style ID. I'm not sure if there's any such identifier in the Google Maps API?

Thanks for any advice.
Ed Saunders


texts:
I am doing the infamous qgistutorials:

Two raster layers (converted from vector to raster) having ("pseudo") band values 0 and 1, were added up with + in the raster calculator.
This resulted in a "merged raster" with values 0, 1, 2; with 2 being the cross section of two 1s of the input rasters.

To remove the 2, the merged raster was put into the raster calculator. again with the expression merged raster > 0.
As a result the latest raster has only 0, 1, (and the cross-section 2 was not represented in the raster image anymore).

Here is, where I have problems with comprehension. What does the calculator do in the last step?
I know that > is a bigger sign. But would not it take all the values from the raster > 0, i.e. take 1 and 2 and remove 0 information?


texts:
In order to help our client manipulate a map to show disparities in public school course offerings in Illinois at the school district level, we’re looking for a current year vector shapefile of the school district map, exactly like what we see here:

https://www.illinoisreportcard.com/District.aspx?districtId=15016299025

The thought would be to illustrate course deserts through creating an interactive and / or static “heatmap” at the district level (technically showing only high school and unit districts, not ALL).

Can anyone help me navigate to a map of this kind? I am struggling to figure out how to find this because I’m just not that savvy with this sort of thing.


texts: Empty PostGIS geometry displaying correctly in QGIS

texts: On a table, with pgadmin (III) I see just emptiness in the geometry column of a table, while the layer is displaying properly in QGIS (2.14). Any idea on why and how to display the geometry in pgadmin ?


texts: Is it possible to only load Leaflet vector tiles at specific zoom levels?

texts:
I am "getting to know" Leaflet + TileStache, and wondering if there is a way to not have GeoJSON tile layers go back to the server on every single zoom level, but rather fetch every, say, 3 zoom levels ... and instead do client-side clipping and zoom-level specific styling (like start rendering labels, etc). of features it already has loaded.

This would save a lot of traffic....

I see a similar question came up a while back:

https://stackoverflow.com/questions/28904160/how-to-use-chunky-coarse-grained-zoom-levels-with-leaflet-maps

but never seems to have been definitively resolved.


texts: Geotools 12.3 - Create DataStore from WFS Service ArcGIS Server

texts:
I'm trying to create a DataStore in Java with Geotools 12.3.
With WFS services in geoserver works fine, but with WFS services in ArcGIS doesn't work at all.

This is the exception:

WARNING [net.refractions.xsi.element] (default task-78) Starting schema with ServiceExceptionReport element.
ERROR [stderr] (default task-78) java.lang.NullPointerException
ERROR [stderr] (default task-78)    at org.geotools.data.wfs.v1_0_0.WFS_1_0_0_DataStore.getSchema(WFS_1_0_0_DataStore.java:368)
ERROR [stderr] (default task-78)    at org.geotools.data.wfs.v1_0_0.WFSFeatureSource.getSchema(WFSFeatureSource.java:204)

This is my Java code:
"Main Class"

public String getDataSelection(String url) throws IOException {

    Map mapResponse = new HashMap();
    String response = null;

    double x1 = strToDbl(this.getBboxT()[0]);
    double y1 = strToDbl(this.getBboxT()[1]);
    double x2 = strToDbl(this.getBboxT()[2]);
    double y2 = strToDbl(this.getBboxT()[3]);

    try {

        String url = "http://sigserv02.anla.gov.co:6080/arcgis/services/Alertas_tempranas_deforestacion/MapServer/WFSServer?service=WFS&VERSION=1.0.0&REQUEST=GetCapabilities";
        int layerName = 0;
        DataStore ds = configWFSDataStore(url);

        FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
        FeatureSource fs = ds.getFeatureSource(layerName);
        FeatureType schema = fs.getSchema();
        CoordinateReferenceSystem worldCRS = DefaultGeographicCRS.WGS84;
        String geometryPropertyName = schema.getGeometryDescriptor().getLocalName();
        CoordinateReferenceSystem targetCRS = schema.getGeometryDescriptor().getCoordinateReferenceSystem();
        ReferencedEnvelope click = new ReferencedEnvelope(x1,x2, y1, y2, worldCRS);

        ReferencedEnvelope bbox = click.transform(targetCRS, true);
        Polygon clickPolygon = JTS.toGeometry(bbox);
        Filter filter = ff.intersects(ff.property(geometryPropertyName), ff.literal(clickPolygon));

        Query query = new Query(layerName, filter);
        query.setMaxFeatures(1500);

        FeatureReader<SimpleFeatureType, SimpleFeature> source = ds.getFeatureReader(query, Transaction.AUTO_COMMIT);

        CoordinateReferenceSystem CRSobjetivo = CRS.decode("EPSG:4326");

        if(source.hasNext()){
            .......
        }
        else{
            response = "empty";
        }
    } catch(Exception e){

        e.printStackTrace();
        ..........
    }
    return response;
}

configWFSDataStore Method:

 public static DataStore configWFSDataStore(String url) throws IOException {
    DataStore dataStore = null;

    BBConectionWFS conectionWFS = new BBConectionWFS();
    conectionWFS.setGetCapabilities(url+"service=WFS&VERSION=1.0.0&REQUEST=GetCapabilities");
    dataStore = conectionWFS.dataStoreWFS();

    return dataStore;
}

BBConectionWFS Class:

public DataStore dataStoreWFS() {
    DataStore dataStore = null;
    try {

        Map<String, String> connectionParameters;
        String getCapabilities = "http://sigserv02.anla.gov.co:6080/arcgis/services/Alertas_tempranas_deforestacion/MapServer/WFSServer?service=WFS&VERSION=1.0.0&REQUEST=GetCapabilities";

        String variableCapabilities = "WFSDataStoreFactory:GET_CAPABILITIES_URL";
        connectionParameters = new HashMap<String, String>();
        connectionParameters.put(variableCapabilities, getCapabilities);

        dataStore = (new WFSDataStoreFactory()).createDataStore(connectionParameters);

    } catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }

    return dataStore;
}


texts: Query WMS GetFeatureInfo with known latitude and longitude

texts:
My issue is that normally, when you want to query a WMS layer, you need to pass width & height of a map, and x,y (in pixels) which describe a point on the map where user clicked.

I want to do this in other way - I want to query this WMS layer with longitude and latitude known before. I assume, that by using some ArcObjects libraries I will be able to convert latitude,longitude and projection to width,height,x,y so that I can use GetFeatureInfo. Can someone please explain how to achieve this?

I'm using ArcObjects SDK 10.1 for .NET.


texts: How to make labels in a fixed position for a road layer in QGIS 3.2

texts: I am working on a road layer containing hundreds of roads and want to print it out in a A2 size with all labels. I want the labels to be fixed in the original position and not depend on scaling of the map. Is there a solution?


texts: Representing channel characteristic data

texts:
I have stream channel cross profiles of a small stream, as well as water depth data from dive loggers we installed the beginning of summer. There is only one dive logger in the stream. I am trying to graphically represent the change in water depth over time within the channel cross section profiles we measured, but I am having difficulty in going about it. I thought using Hec-Ras would be a good start, but being unfamiliar with the program I've become stuck.

How do I go about representing this data?


texts: Using an Expression to Return the Name of Field being Symbolized, within the Layer being Symbolized? (QGIS)

texts:
I'm working in QGIS, formatting maps in their layout view. I'd like to swiftly and accurately make >50 nearly identical maps using a single layer feature with a large attribute file. Each map would be the same map, except for the field (column of attribute table) whose data is being symbolized in the map layout.
I'd like for the name of the field to be automatically updated and displayed on the map in a label box in layout view (or some analogous text display).
It appears inserting an expression into the label item is the answer to my dilemma but there don't appear to be any expressions which return the name of the field being symbolized within the layer being symbolized.  Is this correct?
I have access to ArcGIS desktop and online, should I switch to their programs for this niche functionality? (I'd prefer not if possible haha)


texts: Get x coordinates in a new field in ArcPy?

texts:
In desktop you can do that by: calculate geometry and then select X_Coordinate_of_centroid. How to do this in ArcPy?

Tried:

arcpy.AddField_management("f1.shp","X_Coord","DOUBLE")
arcpy.CalculateField_management("f1.shp","X_Coord") #how to make this work?

It needs also the expression parameter.