texts:
There are a lot of questions for this method, but does not fit in my case.
I have designed a network of roads, manually, and this is optimal. Max longitudinal grade of roads is 12%. (These are forest roads).
I tried to calculate least cost path with r.walk, SAGA Least cost path, and I got a result like this.

This is good result for walking, over the ridge, but not for trucks.
Is there another way to obtain the approximate result with my solution?

On this link I have set an example of the whole project with layers [http://www.mediafire.com/download/me8vgx77vu7vtzw/LeastCost.rar], so if someone wants to try. If someone has a good example to explain it would be good.


texts: Add NEW WMTS from geoportal.gov.pl in OpenLayers v3.9.0?

texts:
i try to add New wmts raster service  from geoportal.gov.pl
https://mapy.geoportal.gov.pl/wss/service/PZGIK/ORTO/WMTS/StandardResolution?version=1.3.1&SERVICE=WMTS&REQUEST=GetCapabilities.
or
https://mapy.geoportal.gov.pl/wss/service/PZGIK/ORTO/WMTS/StandardResolution?version=1.1.1&SERVICE=WMTS&REQUEST=GetCapabilities
Now i use old version (code).
How can i use new service. how to set new  resolutions from new?
//old version
//copy from  https://gis.stackexchange.com/questions/293636/add-wmts-layers-to-openlayers

var myLayer= new ol.layer.Tile({
                opacity: 1,
                name: "" + array_params.layer_name + "_" + array_params.layer_index + "",
                source: new ol.source.WMTS({
                  url: 'https://mapy.geoportal.gov.pl/wss/service/WMTS/guest/wmts/ORTO',
                  layer: 'ORTOFOTOMAPA',
                  matrixSet: 'EPSG:2180',
                  format: 'image/jpeg',
                  projection: new ol.proj.Projection({
                    code: "EPSG:2180",
                    units: "m",
                    extent: [144907.16581514146, 129171.69217334315, 890379.1658151413, 812697.4119729949]
                  }),
                  tileGrid: new ol.tilegrid.WMTS({
                    origin: [100000, 850000],
                    matrixIds: ["EPSG:2180:0",
                      "EPSG:2180:1",
                      "EPSG:2180:2",
                      "EPSG:2180:3",
                      "EPSG:2180:4",
                      "EPSG:2180:5",
                      "EPSG:2180:6",
                      "EPSG:2180:7",
                      "EPSG:2180:8",
                      "EPSG:2180:9",
                      "EPSG:2180:10",
                      "EPSG:2180:11",
                      "EPSG:2180:12",
                      "EPSG:2180:13"],
                    resolutions: [2116.6709,
                      1058.33545,
                      529.167725,
                      264.5838625,
                      132.29193125,
                      66.145965625,
                      26.458386249999997,
                      13.229193124999998,
                      6.614596562499999,
                      2.645838625,
                      1.3229193125,
                      0.529167725,
                      0.2645838625,
                      0.13229193125],
                    tileSize: 512
                  }),
                  style: 'default',
                  wrapX: true
                })
              });


texts:
I have a fresh installation of GeoNode 2.6 on Ubuntu 16.04. The Web site is accessible, however, we are experiencing issues with the "template-overriding" feature.

I'm aware of GeoNode 2.6 problem template overriding.

I have read the doc about "theming" GeoNode. (I could have missed something though, feel free to point it out if that's the case)

From my understanding if I put my HTML sources in /etc/geonode/templates they should override the default sources (that's the behaviour I'm seeking) and the page should update automatically. However that does not seem to be the case. In fact nothing happens (apparently) when I do the latter. I even tried executing the following commands to no avail.

geonode collectstatic
service apache2 restart

My local_settings.py seem all right

GEONODE_ROOT = os.path.dirname(geonode.__file__)
TEMPLATE_DIRS = (
'/etc/geonode/templates',
os.path.join(GEONODE_ROOT, 'templates'),
)

# Additional directories which hold static files
STATICFILES_DIRS = [
    '/etc/geonode/media',
    os.path.join(GEONODE_ROOT, 'static'),
]

Is it possible use my templates by putting them in /etc/geonode/templates (as mentioned here) in order to theme GeoNode?

If possible how can that be done ?

If not possible, why is it not?


texts:
I am importing some layers in postgis. Now I will do a SQL-query and the result should be exported as a MapInfo-file:

SET PGCLIENTENCODING=LATIN1
set MAPSCHEMA=public
set PSQLPATH="C:\Program Files\PostgreSQL\9.5\bin"
set OGR2OGRPATH="C:\Program Files\QGIS 2.14\bin"
set ACDCONN="PG: host=localhost port=5432 dbname=test_mo user=postgres password=admin"
set INPATH_SCHABLONE=C:\Users\LiM\Desktop\batch_ExtractedOcean\schablone
set INPATH="R:\TomTom\Africa\2015_12\mn_inkl_POI\
set OUTPATH=C:\Users\LiM\Desktop\batch_ExtractedOcean\Extracted_MapInfo

echo importing data into POSTGIS

set sIU=reu_oin_oi2
%OGR2OGRPATH%\ogr2ogr -f "PostgreSQL" %ACDCONN% %INPATH_SCHABLONE%\CutAFROceanTT.TAB
%OGR2OGRPATH%\ogr2ogr -f "PostgreSQL" %ACDCONN% %INPATH%%sIU%\Mapinfo\%sIU%_area_orig.MIF

echo exporting  to MapInfo files (tab/mif/mid)...
%OGR2OGRPATH%\ogr2ogr.exe -f "MapInfo File" %OUTPATH%\%sIU%_area.tab %ACDCONN% -sql "SELECT * FROM CutAFROceanTT"

The import is working but the export is not working. The file has been created but the file is empty. I can Export a hole table inside PostGIS but I only want to export the result of the query.

Where is my mistake?


texts:
I have downloaded some DEM from Landsat 8 from usgs website.
when I opened the compressed file I saw some Tiff files and two text files.
I opened the txt file using ArcGIS 10.7 and it contains many other files such as Panchromatic, Water, Thermal,Landcover and etc as following:

I opened Panchromatic bcoz I have heard the Horizontal accuracy is 15 meter.
I want to report slope and elevation of some certain points.
but something is odd. when I read the elevation, the elevation is not correct I think these are in terms of binary values!

plus I downloaded this area from Landsat 7 as following:

the first one is Landsat 7 and the second one is Landsat 8.

*how should I convert Landsat 8 elevations?

*which layers in Landsat 8 I should use to report elevations and slope? (bcoz it contains many Tiff files) is it better than SRTM 1 arc-second?


texts:
Figure A Figure B

If I am generating line geometry with a marker to label parcels when I use a label call-out, (attached images) how can I use the automatically defined label placement point for my anchor?

QGIS draws my line to the centroid as soon as I move the label from it's original position.

Is there a variable that holds this original geometry position that QGIS put the label at?


texts:
I'm trying to get the value of a pixel from a coordinate.

The problem is when I use a large image, so this error appears

  Traceback (most recent call last):   File
  "D:\algoritmo\teste_read_2.py", line 61, in
      teste = band1.ReadAsArray()   File "C:\Python27\lib\site-packages\osgeo\gdal.py", line 2449, in
  ReadAsArray
      callback_data = callback_data)   File "C:\Python27\lib\site-packages\osgeo\gdal_array.py", line 339, in
  BandReadAsArray
      buf_obj = numpy.empty([buf_ysize,buf_xsize], dtype = typecode) ValueError: array is too big; arr.size *
  arr.dtype.itemsize is larger than the maximum possible size.

Follow my code below:

myshape1 = shapefile.Reader('TEMP_APP_BDG_TODAS_2018Point_selec_utm')
pontos = myshape1.shapes()
lon = []
lat = []
for p in pontos:
        xy = p.points[0]
        lon.append(xy[0])
        lat.append(xy[1])
 plt.figure(1)
 plt.plot(lon,lat,'r*',ms=15)
 alt_srtm = []
 arquivo = "mosaico_srtm_utm_5m"
 image = gdal.Open(arquivo+".tif")
 band1 = image.GetRasterBand(1)
 cols = image.RasterXSize
 rows = image.RasterYSize
 print rows,cols
 geo = image.GetGeoTransform()
 print geo
 lat1 = geo[3] # Pego a latitude do primeiro pixel
 lon1 = geo[0] # Pego a longitude do primeiro pixel
 pixel_i = geo[1]
 pixel_j = geo[5] * -1
 print lat[0],lon[0]
 teste = band1.ReadAsArray()
 for k in range(0,len(lat)):
        i = int(((lat1 - lat[k])/pixel_i))
        j = int(((lon[k] - lon1)/pixel_j))
        value = teste[i][j]
        alt_srtm.append(value)


texts: Does QGIS have the equivalent of ArcGIS's Layer (*.LYR) file?

texts:
Using ArcMap, it is possible to add a dataset (eg a shapefile) to the map, and configure the layer's properties (such as the symbology, table joins, switching some attributes off in the Identify results, etc).

This layer can then be saved as a Layer file:

The resulting *.lyr file references the data source (the shapefile) and all of the properties which have been configured for that layer. This makes it easy to add this same layer to multiple map documents, inheriting all of the properties which have been set.

Is there an equivalent in QGIS?


texts: is there any direct access to Mapbox data? (similar to CartoDB's sql api)

texts:
I've been asked to add a form to a map that will write data in the background (e.g., ajax to a magic file that does datawrites and some emailing, notifications etc) and a new marker will be added to the map being displayed.

I have done this before with CartoDB SQL api.  Sent all the form data to a magic ajax script that wrote the new info to the database, then the front-end would refresh the map to show the new marker.

So that's fine, and I had intended upon copying my code over wholesale, but then I found out this is in MapBox ... Is there any way to write to the MB database? or will i need to re-work the entire map storing the data elsewhere ?


texts: Inverse of Selected Records using ArcObjects?

texts: What is the Process to find the inverse of selected data from featureclass using vb.net arcobject?


texts: Road Miles Calculation without Google Maps

texts:
We have two ZIP codes we want to calculate the road mile difference.

Google maps states you cannot use their API if you do not show a map. Is there another service that provides this type of calculation? I have searched Google but nothing has turned up useful.

Currently we use lat/lon to get as the crow flies calculations.


texts: Python code to select streamlines within selected watershed

texts:
I'm trying to calculate total stream length for multiple, overlapping watersheds (i.e., each watershed contains the watersheds upstream of it). I have a flowline layer that I've run through the Intersect tool so that they end on the boundaries of the watersheds (not the overlapping watersheds because that created multiples of the same flowlines). My idea was to create a cursor that would iterate through the watershed table, select the individual watershed, then use the select-by-location tool to select the only the stream segments within the selected watershed. However, every time I try this, I get the error: "Error in executing tool."
Here's the code so far:

watershed = "Sampled_Watersheds"
flowline = "Flowline_Intersect2"

with arcpy.da.SearchCursor(watershed, "Stream_Length") as cursor:
  for row in cursor:
    arcpy.SelectLayerByAttribute_management(watershed, "NEW_SELECTION", [row])
    arcpy.SelectLayerByLocation_management(flowline, "COMPLETELY_WITHIN", watershed, "NEW_SELECTION")

After, I get the selection, I want to add the total stream length to the watershed table. Would something like the follow work?

arcpy.CalculateField_management(watershed, "Stream_Length", "SUM(!Flowline_Streams.Length_km!)","PYTHON")


texts: How to rename a file in the output (ArcGIS python)?

texts:
raslist=arcpy.ListRasters()
for ras in raslist:
    arcpy.ProjectRaster_management(ras,outfolder2+ras+'_Albers',spatialref2)
    print ras+' has been reprojected'

I was using the above code for batch processing in "Project Raster". I want to add "_Albers" in the output file name for each file. But the third line doesn't work. if I delete the "_Albers" in the third line, it works. But I want to rename the output file.

Thanks.


texts: Calculating radius of curves in a polyline with QGIS?

texts: I have a polyline made of curves and I would like to know if there is a way in QGIS to calculate the radius of the curves?


texts: Selecting Similarly Named Text Elements using ArcPy of ArcGIS Pro?

texts:
I have text elements I'm trying to align in ArcGIS Pro 1.4.1. I want to select all text elements with similar names, i.e., "latlonTableLinexx", where xx is the element number. I guess ideally I'd follow that with the "Align Center" command, aligning to the element at the top (usually element 01).

Probably a simple solution but I didn't find anything in a cursory search of ArcPy for selecting similarly named elements. I've certainly worked with elements in Python, but never asked it to select a group so I can work with them in the ArcGIS Pro layout.

Thoughts?


texts: Difference between Publish service and Save service definition file in ArcGIS for Server?

texts:
I have published few map services in which I directly publish on server and for remaining  services I have created service definition file.

All services are working fine but what exactly difference between "Publish a service" and "Save the service definition file" ? . In case "Publish a service " it also analyze the data and gives the warning message..

In case "Save the service definition file" we know the location of .sd file but in other case we do not.. so which one is best and why OR recommended ?

Software : ArcGIS Server 10.1 and Desktop 10.1


texts: Making a grid on a map and identifying points in each cell?

texts:
I have a map of a country with the location stores within the country (a store identifier and the coordinates of the store). I want to:

Create a grid over the map (say, 2 km by 2 km).
Create an index to identify all stores that are located in the same cell of the grid.
Export the resulting data (store identifier, coordinates, and cell identifier).

How do I do this in QGIS 2.18 on Ubuntu 16.04?


texts: GDAL doesn't work on Geoserver for ECW (and CentOS)

texts:
in a server with centos 7, tomcat8, java 1.8.0_121, geoserver v.2.10.1, i followed this link  to install the gdal plugin:http://docs.geoserver.org/latest/en/user/data/raster/gdal.html

i can see the list of the stores including the ecw store. but when i try to select the store geoserver crashes and restart.

I've tryed the gdalinfo command to check the file. With all the files i tested i have the same error:

terminate called after throwing an instance of 'std::length_error'
  what():  basic_string::_S_create
gdalinfo.sh: line 5: 19137 Aborted                 (core dumped) java -classpath "${CLASSPATH}" gdalinfo $input

i've searched in everywhere and it seems to me that this is a bug in the gdal plugin.

Is it possible to view an ECW file in a GeoServer installed on a centos server?


texts: Tutorial for WFS connection SQL query composer

texts:
Has someone a link to a tutorial or a manual that explain the functionalities of the 'SQL query composer' for QGIS?

I'm working with WFS vectors connected from the same server and I would like to download datas after querying that involve different layers to minimizing the downloading datas.


texts: Koppen-Geiger climate classification in R

texts: Are you aware of any R library that does the Koppen-Geiger climate classification, once the necessary input gridded (rasters, netcdf etc.) data are provided?


texts: Area of MultiPolygon including all the holes inside the polygon - Google Earth Engine

texts:
I am facing a very simple problem but my code is getting very complex. There must be some simple solution for this. That's why I am posting.

Problem:
I have a MultiPolygon. In this Multipolygon I have some polygons containing holes (some may not). I want to find the area of the holes and thus find the total area of MultiPolygon including the holes. (i.e.,  Area of outershell of all polygons.)
Is there any simple solution to get the area of Multipolygon including the hole area?

P.S.
This is a simplified version of my problem. In real problem threre are thousands of polygons inside the multipolygon which is the result of reduceToVectors. I know the first coordinate list item of a polygon is the outer shell and all the rest list items are inner hollow shells.

Link to my code
https://code.earthengine.google.com/fd058b8d3864dfec0d1efab741a24a15


texts: Why is ArcGis change precision of double field after copy?

texts:
I create the attribute table of a shapefile and while I copy it by using the "copy features" tool it change the precision from 11 to 14. I create two fields with the name poly_x and poly_y and they represent the centroid of the polygon.

I use ArcGis 10.2.1 and the attribute contain some text fields and two float fields which contain the X and Y coords of the cendroid of the polygon, with precision 11 and scale 3. After I use the above mentioned tool the precision is change to 14 and the scale remain the same.

I would like to mention also, when I insert the shapefiles in QGIS even if I define the precision 11 and scale 3 I see that the precision is 12 and the scale remain the same. The same thing happens when I add the data in Arcgis 9.3.1.


texts: Automatic cutoff circle placement in FIJI

texts: I am working with FIJI (imageJ). I have images of cells stained with DAPI. My goal is to fit a circle with only 12 surrounding cells. It is not too precise to do manually. Any suggestions? I should say that the cells are easily distinguished from one another on the images.


texts: Receiving value from QgsPropertyOverrideButton if data defined override is checked

texts:
I am trying to develop a plugin for QGIS 3. For me it is quite the first time really working with pyqgis or python, so there might be something simple I am missing.

I want to use data defined override button within this plugin. If it is checked, the used value should be taken from this button (expression or attribute) instead of a spinbox. This is my layout with the button:

Unfortunately I could not find any documentation or related questions, which explain how to actually use this QgsPropertyOverrideButton. So I am only able to grab the value from the spinbox. This is my relevant code snipped:

    #Check where to gather attributes from: GUI or Layer?
    if self.dlg.Isochrones_WalkSpeed_Override.isActive() == True:
        WalkSpeed = self.dlg.Isochrones_WalkSpeed_Override.value() #Receiving Value from GUI: DataDefinedOverride
    else:
        WalkSpeed = self.dlg.Isochrones_WalkSpeed.value() #Receiving Value from GUI: SpinBox

How do I grab the value from there either via an attribute from layer or via a given expression within edit field?


texts: WebGL Issues with popup in ArcGIS API 4.17 for JavaScript

texts:
In addition to the main map view in my app, I have a popup that has 2 distinct map views. The popup is launched when ever the use clicks on the map. I am using ArcGIS API 4.17 for JavaScript for my maps, views, popup, widgets, etc... I have noticed an odd bug that I believe is related to the popup or map views in the popup (and/or main view). After the user has clicked the map and launched the popup 6-7 times, on the next click, the base map and popup disappear and all that is left are my widgets and the navbar.

The only warning I get in the console is: WARNING: Too many active WebGL contexts. Oldest context will be lost.. I see that there is some esri documentation on this, but I am already using using 4.x ( as recommended in the article. This seems like a context issue; perhaps the old popup context is remaining. Is there a way to fix this or any recommendations on how to destroy the context on each popup close?


texts: How to deflate only red buildings in Leaflet (L.Deflate plugin)

texts:
I want the L.Deflate plug-in to be applied only to red buildings.

Created L.featureGroup() and put red buildings there, but Deflate does not work.

This code works for all buildings:

L.Deflate({minSize: 10).addTo(map);

If I try to put red buildings in this layer, Deflate does not work:

L.Deflate({minSize: 10, featureGroup: group_red}).AddTo(map);

how can I do this and is it possible, is it at all?

here JSfiddle, demonstrating the problem (uncomment line 15)

https://jsfiddle.net/petrovnn/65y1ed8w/


texts:
I have a huge list of parcel data that I need to assign land use values. The criteria is such that if more than one land use type overlays a parcel, then I need to assign the parcel with the land use that covers it the most (i.e., has the highest LU_Pct). Here’s an example of what my attribute table looks like:

I tried doing this by using Select by Attributes, but it didn’t work. I also tried to use the summarize function which allowed me to find the maximum LU_Pct for each Parcel, but it wouldn’t carry over the corresponding LU_Type record in the final output.

Does anyone know how I can find the maximum LU_Pct for each Parcel that will also carry over the LU_Type in the final selection or output feature (polygon)? I apologize if I'm overlooking something very basic!


texts: Unique constraint on geometry column in PostgreSQL/PostGIS

texts:
I have a table of points, and I set some time ago a simple unique constraint as such :

ALTER TABLE points ADD CONSTRAINT points_geom_key UNIQUE (geom);

It generally worked fine, until now that I face a weird behaviour. I tried to directly copy/paste some points through QGIS with the copy/paste of entities, as I am generally doing, but got an error saying the constraint is not being respected. Long story short : I have two points that are nowher near each other, for which the geom field looks totally different.

I guess this is not the way to check for the unicity of the geom, so what would be the proper way?

I run PostgreSQL 9.4 and PostGIS 2.2


texts: ArcObjects Convert GxDialog to Workspace

texts:
ArcGIS 10.1 with VB2010
I am using the GxDialog to get the users option to save a feature class. The target can be a shapefile, geodatabase, or Dataset inside a geodatabase.
I can’t wrap my head around how I am supposed to parse the Workspace, FeatureDataset, and FeatureClass name to use in creating the feature class.

So far I have

    'setup the save dialog
    Dim pGxDialog As IGxDialog = New GxDialog
    Dim pGxObjFilterSh As IGxObjectFilter = New GxFilterShapefiles          'shapefiles
    Dim pGxObjFilterPg As IGxObjectFilter = New GxFilterPGDBFeatureClasses  'personal geodb (Access db)
    'Dim pGxObjFilterFg As IGxObjectFilter = New GxFilterFGDBFeatureClasses  'file geodb (ESRI db)
    Dim pGxObjectFilterCollection As IGxObjectFilterCollection = CType(pGxDialog, IGxObjectFilterCollection)
    pGxObjectFilterCollection.AddFilter(pGxObjFilterSh, True)
    pGxObjectFilterCollection.AddFilter(pGxObjFilterPg, False) 'this one actually contains both file and personal geodb feature classes
    'pGxObjectFilterCollection.AddFilter(pGxObjFilterFg, False)

    'get the users selection (it can be a shapefile or a geodb). Catalog paths consist of two parts: the workspace and the base name.
    If pGxDialog.DoModalSave(MyApp.hWnd) Then
        Dim catalogPath As String = My.Computer.FileSystem.CombinePath(pGxDialog.FinalLocation.FullName, pGxDialog.Name)
        Dim workspaceCategory As String = pGxDialog.FinalLocation.Category
        Debug.Print(String.Format("catalogPath={0}, workspaceCategory={1}", catalogPath, workspaceCategory))
        Debug.Print("ClassID={0}", pGxDialog.FinalLocation.ClassID)

        'I can use the workspaceCategory to get the target type. But how do I split out either the entire string path or
        'the pGxDialog to create a Workspace object. How do I set up the objects to allow a Feature Dataset?
  Select Case workspaceCategory
      Case "Personal Geodatabase", "Personal Geodatabase Feature Dataset"
          'personal geodb (Access database C:\myData\mypGDB.mdb)
          wkspFactory = New AccessWorkspaceFactory
          pCLSID.Value = "esriGeodatabase.Feature"
      Case "File Geodatabase", "File Geodatabase Feature Dataset"
          'file geodb (File database C:\myData\myfGDB.gdb)
          wkspFactory = New FileGDBWorkspaceFactory
          pCLSID.Value = "esriGeodatabase.Feature"
      Case "Folder"
          'shapefile (c:\mydata\)
          wkspFactory = New ShapefileWorkspaceFactory
          pCLSID.Value = "esriCore.Feature"
      Case Else
          Throw New Exception("Unsupported workspace category: " & workspaceCategory)
  End Select
    End If


texts: Does the projection of a PostGIS table affect the rendering speed in QGIS?

texts:
I am interested to know whether changing the projection of a PostGIS table will affect the rending speed when viewing it in QGIS.
I have two layers coming from QGIS which are exactly the same and cover a geographical area covering England in the United Kingdom. The only difference is their projection:

The first is ESPG: 4326 (WGS 84)
The second is ESPG: 27700 (BNG)

The former (4326) appears to render in 2500ms and the second one in 7500ms.
Additional info

Both layers have a primary key on the same field
Both layers have a spatial index
The QGIS version is 3.4.15 and the project CRS is 27700


texts: PostGIS shp2pgsql-gui with Windows 10 Enterprise crash when running

texts:
Using shp2pgsql-gui on Windows 10 Enterprise to import shapefile data, crashes systematically on some files; for some files, this does not happen.

No error message is sent.

Does anyone know the reason?  Is there an alternate ways to basically import shapefiles?

The files come from our government and are widely used.


texts: Choosing Equal Area Projection for Central Ethiopia (37N)?

texts:
Can anybody please tell me what is the best equal area projection for central highlands of Ethiopia. I am assuming Albers Equal Area Projection.

Any suggestions?

Lambert Azimuthal Equal Area preserves the area of the land therefore minimal distortion in area because of projecting.


texts: Using ST_Value from raster to query second table in PostGIS?

texts:
This is a continuation of Returning raster values from multiple raster tables via Lat/Lon in PostGIS?, but feels sufficiently distinct enough to warrant its own post.

I am trying to select values from multiple rasters based on a single point, but need to match the value returned from one of the rasters to a separate table that holds attribute data.

The code as it exists is:

 WITH mypoint AS (SELECT ST_SetSRID(ST_MakePoint(lon, lat), 4326) geom)
               SELECT
                   ST_Value(a.raster, geom) AS value_a,
                   ST_Value(b.raster, geom) AS value_b
               FROM mypoint p
                   LEFT JOIN raster_a a ON (ST_Intersects(p.geom, a.raster))
                   LEFT JOIN raster_b b ON (ST_Intersects(p.geom, b.raster));

Which does the trick of pulling the multiple raster values. I then want to take one of these values (value_a) to use to filter out a value from a secondary table with two columns "value" and "label." This table provides the text land cover class for a given numeric value.

Unfortunately my PostGIS familiarity is not strong enough to determine how I would join the attribute table in this case. I tried adding one more join as below:

LEFT JOIN table_c c ON (ST_Value(ST_Intersects(p.geom, a.raster)) = c.value);

thinking that I would need to still use the mypoint p intersection, but with no success.

My question is similar to Adding attribute table to raster in PostGIS?, but is in reference to how the actual join is made.


texts: Checking if polylines are connecting using ArcPy?

texts:
The below picture shows the table of a shapefile, including fields of TotalCT, and coordinates of the ends of the line. I want to merge the lines if they are adjacent and only keep the field of TotalCT.

I have no idea about the logical flow of this arcpy code.

Any suggestions?


texts: Finding duplicate values that ALL share same attributes?

texts:
I am using Arcgis Desktop 10.2. Using python script i am trying to resolve the following :

I have 2 columns of data, one containing alphanumeric values, the other some strings like "-1", "-2", "-3" etc.

"-1" means that the data in the first column in no longer needed. The data set is large, millions of rows/ I am supposed to find duplicate string in the first columns. For every match , if all of the rows are marked with "-1" then that value will be removed, or marked by a commentary in another column. If for the match one value(or more) are not "-1" then that value will ignored

Looking at the picture:

"string_1" is found 5 times. Because the last instance is not marked "-1" all the 5  "string_1"'s will remain in the table.
"string_2" is found 4 time. The last values is not marked with "-1", but with "-4", so all the 4 values remain in the table

"string_3" is found 5 times, and all of this values are marked with "-1". For each of this values a commentary must be made in another column, so that in the future i can remove them form the table.

How can i mach the values in the first column (i need the script to find for itself the duplicate values, because i can't tell it what it is supposed to find, there are millions of values) and then find the duplicates that ALL have "-1" attribute?


texts: Download USGS satellite images of a city only

texts: I need to download the satellite images of a city for land change detection. Specifically I need to show urbanization trends and predict the future expansion. Now I can't seem to download just the required tile. Insted USGS Earth Explorer is showing the downloadable area that is exceeding my required area or is missing it.


texts: Problems serving GeoTiff using GeoServer and ImageMosaic

texts:
Goal:
I am trying to serve a large (1GB+) GeoTiff file using GeoServer and ImageMosaic.

Process:
iant have written a simple guide for this, which he proposed as an answer for a previous and similar question.

I have used his gdal_retile example and this gave me a folder with lots of tiff files and 8 subfolders (1-8, level 1 having the most zoomed tiles and level 8 a small overview tile).

I created an ImageMosaic store in GeoServer 2.1.3 and pointed my URL to the folder containing the tiff files and the 8 subfolders. GeoServer then put 8 files in the same folder.

I published my ImageMosaic store as a layer and when doing a preview inside GeoServer, Java just starts working and I don't know if it will ever finish.

What am I doing wrong here?


texts: Replacing all "NULL" values in multiple rows with "0"

texts:
I am looking for a way to replace all "NULL"-Values throughout multiple rows/ all "NULL"-Values in an attribute table. The table is from a SHP-File (shapefile).
Is there any easy way? Right now I am selecting then by hand and replace them row by row which takes quite some time. I have already tried to look up a solution, however as I want to replace them in the whole table it seems to be not that easy.


texts: ArcGIS method to partially erase or clip thousands of related polygons

texts:
I am looking for a way to remove part of a polygon that overlaps with a polygon in another layer.
Specifically I have two datasets, one is a set of square polygons and the other is a set of convex hull polygons which originates from each square polygon. I now need to remove from each convex hull that part which intersects with the square it originates from. The part of the convex hull not within the originating square should remain intact as it will be dealt with in a subsequent process.

The square polygons, each is used to create a convex hull.

The convex hulls. The only part to erase is the circular area within each square.

This is my result from Model Builder. Correctly removed part within the square but only one, I need over 50,000 removed.

Any help improving this Model or do I need a different approach?


texts: Inserting X,Y to Geometry in Spring Boot and PostgreSQL?

texts:
I use @query to Insert X,Y to Geometry in PostgreSQL.

@Modifying
@Query(value = " insert into mypoints (name,geom) values (:name,SRID=4326;POINT(:x :y))", nativeQuery = true)
@Transactional
void addSPCTY(@Param("name") String ten,@Param("x") String x,@Param("y") String y);

Right now it ends up with an error

  Error in Java: org.postgresql.util.PSQLException: ERROR: syntax
  error near ;


texts: Selecting/Exporting points from one layer that do not overlap with the points from another layer in Qgis

texts:
I am working with QGIS and have the following problem to solve:

I have 2 point layers. Layer 1 has more points than Layer 2 but all points of Layer 2 have the same position as Layer 1. I need to find those points that do not overlap.
Does anyone have an idea how to do this?


texts: Identifying convex and concave portions of polygon using ArcMap

texts:
I have a polygon shapefile depicting an island. For my study, I'm interested in two coastal features of the island: convex sections (portions of the coastline that extend into the ocean) and concave sections (portions of the coastline that recede into the island, like an embayment or lagoon).

To accomplish this, I have completed the following steps:

Convert the polygon into a polyline
Generate points along the island polyline boundary at 5km intervals
Split the polylines at the points, to create 5km polyline segments along the island
Used the Calculate Sinuosity tool in ArcMap, which measures deviation of a line from the shortest path, on all 5km polyline segments

However, the Calculate Sinuosity tool doesn't distinguish how the line deviates from the shortest path. The direction of the line deviation (convex: toward the ocean; concave: receding inland) is relevant.

Is there a tool to calculate whether the line deviation falls under a more concave or convex shape for the entire island ?


texts: How is "Burn Stream Slope" in ArcHydro working?

texts:
I have a DEM with only small relief and many ditches. Because of the small relief, I can't calculate flowdirections for some parts of my DEM. I know that in reality water is flowing from west to east (there are pumping stations guaranteeing this). This is a detail of the DEM with 3 points and the ground level from the raster:

To ensure that water flows from west to east in ArcMap as well I

1. added a line following the ditch

2. run "Generate From/To Node for Lines" from the "Attribute Tools" in the ArcHydro toolbar

3. run "Assign Stream Slope" from "DEM Manipulation" from the "Terrain Preprocessing" menu

This is the resulting line and attribute table containing the values entered in the ArcHydro tools:

Afterwards I used the "Burn Stream Slope" tool hoping to get a DEM with a ditch flowing from -1,8 m in the west to -2,3 m in the east.

Result:

The difference between the two ends of the ditch is correct, the elevations aren't. Especially the value of 10013,1 m for the point located north of the ditch make me wonder.

I tried reproducing the single steps mentioned in the help for Burn Stream Slope:

1. Retrieve the step size from the XML. It is defined as a fraction of cell size and default to 0.75.
Seems to be working, see step 3.

2. Retrieve the minimum elevation of the input DEM.
That is: -1,87214

3. Create Edit Points for each from node and to node, as well as for point located along each line using the previous step size as interval.
That's working, points are created

4. For each point stores the current DEM elevation of the underlying cell as well as the new elevation computed by linear interpolation of the FromElev and ToElev along this line feature.
That's working, values for the easternmost point:
ElevOld: -1,430757
ElevNew: -2,3

5. Set the elevation of each cell under an Edit Point feature to the Minimum elevation of the input DEM – Maximum FromElev from Stream feature class + New Elevation stored associated to the point. Basically, the elevations along the streams will be dropped below the minimum elevation of the input DEM and will be decreasing towards the outlet.
This is what I calculated for the easternmost point: -1,87214 - -1,8 + -2,3 = -2,37214
This is the value the tool assigned to the underlying cell: 10,700195

6. Raise the DEM by maximum drop value + 10 to ensure positive values in the DEM.
I have no clue what is meant by "maximum drop value" so I can't reproduce that but I think the problem starts at step 5.

What am I doing wrong?

Since the "Burn Stream Slope" tools seems to be unable to just modify the DEM in the stream regions without touching the surrounding area (see first answer): If anybody has an idea how to accomplish this with other tools I'm willing to assign the bounty to this person.


texts: How to create GeoJSON constructor in openlayers 3?

texts:
I wan't to add several GeoJSON-overlay layers from a Postgis-database to an Openlayers 3 map.
I want to define the source as well as the style, the layer is rendered in ol 3.
Actually I'm using a slightly modified php-script from Bryan McBride PHP PostGIS to GeoJSON, that uses a GET-call to pass the options of the different layers to the script to define the source of my layer.

n2k_do_l = new ol.layer.Vector({
          title: 'Title of Layer',
          visible: false,
          style: new ol.style.Style({
                fill: new ol.style.Fill({color:'rgba(100,250,0,0.1)'}),
                    stroke: new ol.style.Stroke({
                        color: 'magenta',
                        width: 2
                    })
              }),
      source: new ol.source.Vector({
          projection: 'EPSG:4326',
              attributions: [new ol.Attribution ({html: '<br />Attribution of Layer' })],
              strategy: ol.loadingstrategy.bbox,
          loader: function(extent, resolution, projection) {
          var extent = ol.proj.transformExtent(extent, projection.getCode(), ol.proj.get('EPSG:4326').getCode());
          $.ajax({
          type: 'GET',
          url: "./mapdata/get_geojson.php?"+
              "geotable=n2k_do" +
              "&fields=gid as id,sitecode,sitename" +
              "&bbox=" + extent.join(','),
              context: this
              }).done (function (data) {
                 var format = new ol.format.GeoJSON();
                 this.addFeatures (format.readFeatures(data,{
                 dataProjection: 'EPSG:4326',
                 featureProjection: 'EPSG:3857'
                 }));
          });
          }
             })

});

For one layer this is ok, but, I would like to use the code from above, to write an extended ol.layer.Vector class (myol3Vector(options)) to create several layers by passing a list of options:

n2k_bd_l = new myol3Vector({
            title: "Natura 2000 Birds Directive",
            attribution: "Attribution text",
            geotable: "n2k_bd",
            fields: "gid as id,sitecode,sitename,surfha",
            where: ""
});

n2k_hd_l = new myol3Vector({
            title: "Natura 2000 Habitats Directive",
            attribution: "Other Attribution text",
            geotable: "n2k_hd",
            fields: "gid as id,sitecode,sitename,surfha",
            where: "sitename ilike '%moselle%'"
});

How can I transform my working code to a constructor, that creates multiple layers just by passing a list of options?

Maybe someone can get me a hint, how to start.


texts:
I have been looking for it for a while but couldn't find what exactly I need.

I have WRF outputs in NETCDF format and I want to convert them to Geotiffs with gdalwarp.

I think my -s_srs should be this one because I used lambert : http://spatialreference.org/ref/sr-org/29/
and my -t_srs should be this one since I want albers : http://spatialreference.org/ref/sr-org/6630/

By the way I want the QFX subdataset.

When I do this I receive an error:

gdalwarp NETCDF:wrfout_d01_2007-06-01:QFX -s_srs "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -tr 12000 12000 -of GTIFF output.tiff

and the error is this:

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
ERROR 1: Unable to compute a transformation between pixel/line
and georeferenced coordinates for NETCDF:wrfout_d01_2007-06-01:QFX.
There is no affine transformation and no GCPs.

So it requires some GCP but I do not know how to provide them.
Also the latitude and longitude is in the NETCDF as subdatasets.


texts:
I'm using Googlemaps service to plot some marker on the map. All markers are stored in a postgresql table, named "locations" and that has a location column (data type to geometry).
When the user clicks on the markers a route is created between them.

My question is: how can I generate a buffer (with a limit of 10km along route) and select all my locations inside the buffer. I`m really new in postgis.

P.S. This link inspired me but without success
Thanks!


texts: Getting ArcGIS.com-created map to show up in sample app

texts:
I've been experimenting with the My First iOS app example, trying to get it to display a Bing map layer I created on ArcGIS.com.  However, the layer doesn't display.  In the sample app's code, when it says:

NSURL *mapUrl = [NSURL URLWithString:@"http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"];
 AGSTiledMapServiceLayer *tiledLyr = [AGSTiledMapServiceLayer tiledMapServiceLayerWithURL:mapUrl];
 [self.mapView addMapLayer:tiledLyr withName:@"Tiled Layer"];

what am I supposed to replace the mapUrl from?  The bit.ly link created through the share function on ArcGIS.com doesn't work, since it links directly to a map viewer displaying the map.  I'd like to get the map directly to put into my sample iOS app.


texts: Can GeoServer serve hundreds of thousands of geoTIFFs?

texts:
I'm working on a project that is creating several hundred thousand GeoTIFFs.

Each GeoTIFF is a projection of species survivability, for the intersection of a climate future and a species. Unlike many of the "lots of GeoTIFFs" problems people have, these are all independent layers that will be viewed one at a time, and each one spans the world... so they cannot be mosaic-ed into a single big layer.

I read on the GeoServer mailing list that having tens of thousands of layers makes GS go slow, so I suspect it will be impractical for me to add all 300,000 layers to the catalog.

For an earlier version of this project, I scripted a web system around MapServer such that it would dynamically grab the source GeoTIFF on first request, then cache it for a while for subsequent tile requests.

Can I do something similar with GeoServer, or should I go back to trying to compile a working recent version of MapServer and MapScript/Python?

Edit: apparently I can use a mosaic with an extra dimension, like "species id", to treat all my GeoTIFFs as a single catalog entry. Will this work for millions of GeoTIFFs?


texts:
I have an enterprise Geodatabase (SQL Server) on my ArcGIS server machine. I perform my computations on several rasters on another machine (80 *.tif s with mean size of 50kb). When I want to add the result rasters to a mosaic dataset on my ArcGIS server machine using python, it performs very slowly and I have to wait for about 20 minutes.

for filename in glob.iglob('*.tif'):
    raster_name = 'Rater-' + raster_num +'.tif'
    arcpy.AddRastersToMosaicDataset_management('D:/example.sde/SampleMosaicDataset', "raster dataset", 'D:/GRIDS/' + filename)

How can I do this task with more speed?


texts: Filter Image Collection by date with conditions

texts:
I want to filter my ImageCollection for my analysis with specif time ranges. I have an Image Collection with annual max NDVI images ranging from 1982 - 2015 and a conflict FeatureCollection for a period from 1989 - 2015. I'm calculating slopes with a linearFit function before and after a conflict date. Now I want to introduce some conditions for the time ranges before and after the conflict.

The slopes (NDVI trends) need to be calculated for at least 6 years before and 6 years after the conflict date (conflictDate) (to be able to calculate an accurate slope for the vegetation)
The time range needs to be the same, e.g 7 years before and 7 after the conflict (but min 6 years)
The time range should be as long as possbile (the maximum number of years). Because the longer the time range the more images used and the more accurate the slope. So basically it means if the conflict date is 1997. Then the beforeSlope time range is from 1982 to 1996 (15 years) and the afterSlope time range is from 1997 to 2011 (15 years).

//Dates for filtering the GIMMS NDVI Image Collection
var startDate = '1982-01-01'
var endDate = '2015-12-16'

//Loading and filtering (by date) the UCDP conflict data (CSV)  1989 - 2015
var conflicts = ee.FeatureCollection('users/welp1992/Conflicts25_Date')
  .filterMetadata('date_start', 'greater_than', '1988-01-01')//min 6 years, but as long as possible
  .filterMetadata('date_start', 'less_than', '2010-01-01')

//filter image GIMMS by date and band and add time 't' as a band
var timeSeries = ee.ImageCollection('users/welp1992/GIMMS_Max_Annual_NDVI')
    .filterDate(startDate, endDate)
    .select(['b1'])
    .map(function (image) {
      return image.addBands(
        ee.Image(image.getNumber('system:time_start')) //Time past in millis. since 1970 for the composite (15 days) start period
          .subtract(ee.Date(startDate).millis())
          .divide(1000*3600) // year (*24*365 = days)
          .int() //casts the input value to a signed 32-bit integer
          .rename('t')
      )
    })

print(timeSeries,'time');

//Call function for adding slopes to conflicts feature  collection
var conflictsWithSlopes = conflicts.map(addSlopes)
print(conflictsWithSlopes.limit(500)) //(5815 conflicts)
print('slope change stats', conflictsWithSlopes.aggregate_stats('slopeChange'))
print('histogram', conflictsWithSlopes.aggregate_histogram('slopeChange'))
// Other aggregate_*() functions...

//Function to add slopes (3) to feature collection, while calling getSlope function
function addSlopes(conflict) {
  var conflictDate = conflict.getString('date_start')
  var beforeSlope = getSlope(timeSeries, conflict, startDate, conflictDate)
  var afterSlope = getSlope(timeSeries, conflict, conflictDate, endDate)
  return conflict
    .set('beforeSlope', beforeSlope)
    .set('afterSlope', afterSlope)
    .set('slopeChange', afterSlope.subtract(beforeSlope))
}

//Function to calculate slopes using linearFit and reduce GIMMS to intersections with the conflict geometry
function getSlope(timeSeries, conflict, startDate, endDate) {
  return timeSeries
    .select(['t', 'b1'])
    .filterDate(startDate, endDate)
    .reduce(ee.Reducer.linearFit()) //try other reducers ee.Reducer.linearRegression()
    .select('scale')
    .reduceRegion({
      reducer: ee.Reducer.mean(), //different reducers possible (first, min, max, mean)
      geometry: conflict.geometry(),
      scale: 8000
    })
    .getNumber('scale')
}

I'm not sure where to add these conditions. But I think in the getSlope function (bottom) I should introduce them. Somehow the 't' (as a time range) needs to be filtered so that the 'beforeSlope'  time range equals the afterSlope time range. And the max. possible time range.
Code:
https://code.earthengine.google.com/046eaa48b06f647f53f061583d00a0b3
Asset NDVI and Conflict:
https://code.earthengine.google.com/?asset=users/welp1992/NDVI_Annual_Max_81_2015
https://code.earthengine.google.com/?asset=users/welp1992/Conflicts25_Date


texts: Using R to sum raster pixels of specific values

texts:
I have a land use raster (lu_raster) and several thousand shapefiles (counties). Within each county, I'd like to calculate the proportion of agricultural land use. So, I think I need to find a way to 1). sum all the raster pixels that are either value 81 or 82 (see here for reference) within each of the shapefiles, and 2). divide that by the total number of pixels in each of the county shapefiles.
I've been using exact_extract to do this kind of stuff. So I think the second part of my code will look something like this:
Results <- exact_extract(lu_raster, counties, 'count')

But I can't figure out how to do those counts conditionally (only counting the 81s and 82s).


texts: I have a polygon (pond) and I need to ensure it is at one elevation. I have converted all vertices of the polygon to points, used summary statistics to determine the range in elevations and if the range is greater than 0, it indicates there is a discrepancy in elevation values. I do not have spatial analyst or 3D analyst. Is there an easier way or code to do this in ArcGIS?  I received this layer from a fellow employee that extracted the data from a LiDAR point cloud using MicroStation. Therefore, I do not have the original elevation (raster) layer, just a shapefile that has elevation values attached to the vertices. I am looking for a way to identify a polygon whose vertices have different elevation values. Each polygon should have one elevation value associated with each vertice. For example, Polygon A should have a 120m elevation value for each vertice in the polygon. I am looking for something to check if the vertices are different.


texts:
I'm trying to use Qfield for tree surveying. I am testing it on a Galaxy S7 running Android 7. I have installed from Google Play. I can't seem to get QFIELD working. So far, the stumbling blocks are:

When I open QField, I get a alarm icon in the top right corner. On clicking, the message reads "Can't open database: unable to open database file".
I have created a portable project using the plugin to QGIS. The portable project is made up of two files, one with a .qgs extension, the other with  a .qgs~ extension. I have placed these in the QField directory, specifically MyFiles>SD Card>Android>data>ch.opengis.qfield. When I navigate to this location using the Qfield_Open Project, the file is not visible.

It seems like I'm missing something very basic, Can anyone help...

Regards
O


texts:
How can I hide a column in qgis3 using pyqgis:

In Qgis2
it was quite simple:

        editFormConf = layer.editFormConfig()
        idx=layer.fields().indexFromName(columnName)
        editFormConf.setWidgetType(idx, "Hidden")

I investigated, and now in Qgis3 it has changed to something like

        idx=layer.fields().indexFromName(columnName)
        editor_widget_setup = QgsEditorWidgetSetup(??????)
        layer.setEditorWidgetSetup( idx, editor_widget_setup )

The problem is the use of the constructor QgsEditorWidgetSetup(??????)

How to tell the constructor that this field is hidden


texts: Using Data Driven Pages for Page Number with Count?

texts:
I have a set of 100 data driven pages showing different roads.  I would like to export the map book by road.  E.g. one road will be one pdf, another road will be exported to another pdf.

The issue I have is with the page numbering.  I have set each page to have a dynamic page number with count.  At the moment, when I export Road A, the page numbers show 'Page 1 of 100', 'Page 2 of 100' etc.

I want Road A to have 'Page 1 of 5', 'Page 2 of 5', etc. then Road B will have 'Page 1 of 2', 'Page 2 of 2'.

I can't seem to work out how to do this, I'm sure I have done it before, but my terrible memory has failed me!


texts: Calculating new field works inside QGIS Python console but not in standalone program

texts:
My code in the standalone program

import sys

sys.path.append("C:/OSGeo4W64/apps/qgis-ltr/python/plugins")

from qgis.core import *

from PyQt4.QtCore import *
from PyQt4.QtGui import *

import processing
from processing.core.Processing import Processing

QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis-ltr", True)
qgs = QgsApplication([], True)
qgs.initQgis()

Processing.initialize()
Processing.updateAlgsList()

# Here on
shp1 = QgsVectorLayer("C:/Users/Documents/shape.shp", "shp1", "ogr")

processing.runalg('qgis:fieldcalculator', shp1, "new_field", 0, 1, 2, True, "$Area", "C:/Users/Documents/shape2.shp")

qgs.exitQgis()

The exact same code after the label # Here on above works when using the Python console within QGIS. Whereas for the standalone program I get the field names but no values. I'm using version 2.18.17. Is this a bug with the software or am I doing something wrong?

Bug report submitted at: https://issues.qgis.org/issues/18395


texts: Displaying 2 different color lines based on 2 different table attributes in QGIS?

texts:
I have a line shapefile that has a bunch of data and it is broken down into different features which is assigned a different number. It is then styled on the rule based approach so each feature gets a different symbology.

What I am wanting to do is to take one feature (skid trail - "feature" = '11') and have 2 or 3 different color of lines based on what the corresponding column "label" says. In the picture you can see my attempt in the expression string builder but I'm not very good at it. I also want to make sure that if there is no value or another word in the "label" field that it will still display a line.

I know that I could add another rule with a different feature number (ex. feature = 44) but because this data is shared between others we keep the data info the same.

As I was typing this I was wondering if I could even just do this under the rule based symbology style and not even use the expression builder.


texts: Blank map when working with OpenLayers and Geoserver

texts:
I'm trying to develop a map with 4 layers: two shapefiles (layers: "roads" and "hospitals") and two World Image Files (layers: "hypsometric","shadows"). I installed OpenLayers and tried to write correct code, but something is wrong, because I get only blank map when I enter address: http://localhost:8080/geoserver/www/wms.html.

This is my code:

<script src="http://localhost:8080/geoserver/openlayers/OpenLayers.js"></script>
<script type="text/javascript">
    var map;
var layer1,options;
     function setHTML(response) {
        document.getElementById('nodelist').innerHTML = response.responseText;
    };

    OpenLayers.IMAGE_RELOAD_ATTEMPTS = 5;
    OpenLayers.DOTS_PER_INCH = 25.4 / 0.28;

    function init(){
        var bounds = new OpenLayers.Bounds(
            564840.1, 5075229.7,
            588160.6, 5093252.2
        );
        var options = {
            controls: [],
            maxExtent: bounds,
            maxResolution: 0.025439315679878505,
            projection: "EPSG:4326",
            units: 'degrees'
        };
        map = new OpenLayers.Map('map', options);

        layer1 = new OpenLayers.Layer.WMS(
            "Roads",
            "http://localhost:8080/geoserver/wms",
            {
                width: '800',
                srs: 'EPSG:4326',
                layers: 'roads',
                height: '450',
                styles: '',
                transparent: true,
                format: 'image/png'
            },
            {isBaseLayer: false,singleTile: true, ratio: 1}
        );
        layer2 = new OpenLayers.Layer.WMS(
            "hospitals",
            "http://localhost:8080/geoserver/wms",
            {
                width: '800',
                srs: 'EPSG:4326',
                layers: 'hospitals',
                height: '450',
                styles: '',
                transparent: true,
                format: 'image/png'
            },
            {isBaseLayer:false,singleTile: true, ratio: 1}
        );
        layer3 = new OpenLayers.Layer.WMS(
            "hypsometric", "http://localhost:8080/geoserver/wms",
            {
                width: '800',
                srs: 'EPSG:4326',
                layers: 'hypsometric',
                height: '450',
                styles: '',
                format: 'image/png'
            },
            {singleTile: true, ratio: 1}
        );
        layer4 = new OpenLayers.Layer.WMS(
            "shadows", "http://localhost:8080/geoserver/wms",
            {
                width: '800',
                srs: 'EPSG:4326',
                layers: 'shadows',
                height: '450',
                styles: '',
                transparent: true,
                format: 'image/png'
            },
            {isBaseLayer:false,singleTile: true, ratio: 1}
        );

        map.addLayers([layer1,layer2,layer3,layer4]);

       // setup controls and initial zooms
        map.addControl(new OpenLayers.Control.PanZoomBar());
        map.addControl(new OpenLayers.Control.Navigation());
        map.addControl(new OpenLayers.Control.Scale($('scale')));
        map.addControl(new OpenLayers.Control.MousePosition({element: $('location')}));
        map.addControl(new OpenLayers.Control.LayerSwitcher());
        //map.addControl(new OpenLayers.Control.OverviewMap());
        map.zoomToExtent(bounds);

        // support GetFeatureInfo

    }
</script>`

I'm working with OpenLayers for the first time.


texts: Filter Leaflet markers by input selectionin R Flexdashboard

texts:
I am trying to filter the markers that are drawn on a Leaflet map based on what is selected in a user defined date range. The below codes is returning an error that it is unable to find hfxLoc object. I have used the below approach (filtering within the render  function in other elements (i.e. renderPlot, renderValueBox, etc.). Is my approach to filtering correct, or does it need to occur in a different place in the Leaflet code block?
The inputSelect code is as follows:
selectizeInput("shiptypeInput", "Vessel type: ", choices = hfxVessels, selected = "Cargo ships", multiple = TRUE)

The leaflet code is as follows:
renderLeaflet({
  hfxLoc <- hfxETA %>%
    dplyr::filter(as.Date(eta_date) >= input$dateRange[1] & as.Date(eta_date) <=  input$dateRange[2]) %>%
  leaflet(data = hfxLoc) %>%
  setView(lng = -60.25, lat = 46, zoom = 6) %>%
  addProviderTiles(providers$CartoDB.Positron,
                   options = providerTileOptions(minZoom = 2, maxZoom = 16)) %>%
  addMarkers(lng = ~lon, lat = ~lat, clusterOptions = markerClusterOptions())


texts: Is it possible to use DotNet to show a QGIS map?

texts:
Has anybody ever tried to show and interact with a QGIS map from a DotNet language like C#? What I have in mind is, to write a desktop application with C#, and show an embedded map with the QGIS-API on a form.

The most useful document I found, is this Coding Guide, and it seems that QGIS supports exclusively C++ and Python to write applications/plugins. So my question, is it doable to create DotNet wrappers in C++, and use them in a managed language (with a reasonable amount of work)? Any positive or negative experience?


texts: offline map tiles in standalone desktop app

texts:
Which mapping API can be used (either in a C# desktop app or as a standalone mapping app) that accesses an offline tile collection and requires only a minimal amount of additional servers / other applications?

We are currently developing a desktop app in C#, which enables the user to do a fuzzy search on a csv and outputs a place name, for which I have generated two columns of lat/long values.

Now we would like to add map functionality with mouseover labels and give location/district/province, next to a table view of the search result (name/district/province). This map can either be implemented as a window in the search app or as another standalone. I know how to generate a collection of tiles and use these in a web app, but I am stuck on how to bring a map window into a desktop app without the need of having to run all sorts of servers. Our users will not be able to get online and as they will not be very IT literate, we are trying to minimise the amount that they need to install (bundled for Windows) or start up before using the app.

I've so far had a look at jxMapViewer (which we probably can't use as we are developing in C#): https://www.java.net//pub/a/today/2007/10/30/building-maps-into-swing-app-with-jxmapviewer.html


texts:
Till now, I was using PostGIS-pgrouting pgr_tsp function to solve Travelling Salesman Problems. This function can give us the best possible route depending upon the cost parameter when the starting, finishing and via points are given.

But soon I realized to make my function more flexible. For example, in pgr_tsp function you provide the starting and finishing points along with the via points. But what if I am not sure about the exact via point? If I want to calculate the route starting from my home and going to the office via McDonald Restaurant and Shopping Mall, I am not explicitly defining which McDonald or which Shopping Mall? All I want is my route going from the Home to the Office via atleast one McDonald and one Shopping Mall; and has to be the shortest one depending upon the cost parameter I am using (cost parameter could be road length or travel time).

Literature reviews suggested me that Traveling Purchaser Problem is what I am looking for and Travelling Salesman Problem is a special case of this problem.

How can I solve the above explained Traveling Purchaser Problem in PostGIS? I came across the jsprit open source toolkit (a java based product) and MATsim and these other projects, but I have no idea how to use them? Am I supposed to write a PL/Java procedural function to implement this jsprit toolkit in my PostGIS Db or what?

Any initial hint could jet propel from progress.


texts:
The value part of a GP Result can be any type of GP object (or literal).

Does anyone know of a general way to deserialize this, preferably with either WCF or Json.NET?

I'd prefer not to have to write a deserializer for each different GP tool.

{
   "results":[
      {
         "paramName":"Output",
         "dataType":"GPFeatureRecordSetLayer",
         "value":{
            "geometryType":"esriGeometryPolyline",
            "spatialReference":{
               "wkid":4326
            },
            "features":[
               {
                  "geometry":{
                     "paths":[
                        [
                           [
                              -91.6849975582907,
                              24.206096649241
                           ],
                           [
                              -91.8398132323948,
                              24.4056720735074
                           ]
                        ]
                     ]
                  },
                  "attributes":{
                     "FID":1,
                     "FNODE_":0,
                     "Shape_Length":0.25258314060805
                  }
               }
            ],
            "exceededTransferLimit":"false"
         }
      }
   ],
   "messages":[

   ]
}

I never figured out a way to do this.  If Esri wants to encourage consumption of GP services by a wider variety of clients they should consider serializing the .NET type name, similar to what can be done with Json.NET by setting TypeNameHandling = Objects.


texts:
Does the ArcGIS Pro feature "least cost corridor" have an equivalent module or workflow in QGIS?

http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/creating-a-least-cost-corridor.htm


texts: Determining US zip code location information from first three zip code digits?

texts:
I am working on a dataset that comes from survey results.  The only location identifying information we have from those who took the survey is the first three digits of the U.S. zip code.

How detailed of location information does this give?


texts:
My legend is too big and it is blocking the view of my map.

The size of the map is fixed, so the only solution is to make the legend smaller? Is there any way to do this?

    Cases reported are decreasing or circulation is ongoing

    Cases reported are increasing or peaking

    Risk of spreading

    Awaiting Confirmation

New code:

  <div.cartodb-legend.custom {
width: auto; height: 50px;
  }

 div.cartodb-legend.custom.bkg {
font-size: 50%;
  }
 div.cartodb-legend.custom ul li.bkg .bullet {
 background-size: 15px 15px !important;
 background-repeat: no-repeat !important;
  }

    Cases reported are decreasing or circulation is ongoing

    Cases reported are increasing or peaking

    Risk of spreading

    Awaiting Confirmation


texts: QGIS Field calculator / Multiple CASE operator and % expression

texts:
I can't seem to be able to troubleshoot my expression here. I'm trying to calculate a new field. Column "FOO_U_CODE" contains numerous 4 digit codes, the new column I want to create is a simplified version based on the first character of the code in "FOO_U_CODE".

I saw somewhere that you can use % to specify anything after one or multiple characters. Is it right?

What is wrong here?

   CASE
    WHEN "F00_U_CODE" IS'1%' THEN "Utilisatio"  = 'Residentiel'
    WHEN "F00_U_CODE" IS'2%' THEN "Utilisatio"  = 'Transports communications et services publics'
    WHEN "F00_U_CODE" IS'3%' THEN "Utilisatio"  = 'Industrielle'
    WHEN "F00_U_CODE" IS'4%' THEN "Utilisatio"  = 'Transports, communications et services publics'
    WHEN "F00_U_CODE" IS'5%' THEN "Utilisatio" =  'Commerciale'
    WHEN "F00_U_CODE" IS'6%' THEN "Utilisatio"  = 'Services'
    WHEN "F00_U_CODE" IS'7%' THEN "Utilisatio"  = 'Culturelle, recreative et loisirs'
    WHEN "F00_U_CODE" IS'8%' THEN "Utilisatio"  = 'Production et extraction des richesses naturelles'
    WHEN "F00_U_CODE" IS'9%' THEN "Utilisatio" = 'Immeubles non exploites et etendues deau'
   END


texts: Is there a way to create layer summary Pie Charts?

texts:
In Qgis2.6. I want to create a pie chart to summarize the layer data. I have a layer of land use polygons (woodland, agricultural, urban etc.). I want a pie chart to show the % of each land use in the layer.
I see the diagrams section of layer properties will provide pie charts for each polygon, and the 'Group Stats' plugin will calculate the totals I need. All tantalizingly close.
Does anyone know of a way this can be done?


texts: QGIS plugin 'Module' object has no attribute 'classFactory'

texts:
I am trying to learn how to make QGIS plugins and have been trying to follow tutorials such as this one;
    https://anitagraser.com/2014/04/26/getting-started-writing-qgis-2-x-plugins/
    but I keep getting this error right in the first few steps

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISES~1/apps/qgis/./python\qgis\utils.py", line 306, in startPlugin
    plugins[packageName] = package.classFactory(iface)
AttributeError: 'module' object has no attribute 'classFactory'

I did not change any code in the utilis.py so I am wondering if anyone knows what might be the issue and how to fix it?


texts: Need help making polygons from lines in QGIS

texts: I am trying make a layer of buildings, divided into apartments, by making polygons from crossing & intersecting lines in QGIS.  When I use the vector/geometry tools/lines to polygons tool, I get just one polygon with no divisions at the inner lines, or otherwise odd polygons.  The attached screenshot shows a random example of three polygons which were created, with the associated lines I hoped to use for the divisions.   I also then tried the single part to multipart, but could not make that divide the polygon by the separate line file.  I don't know if I am simply using the tools incorrectly, or if there is some setting I am missing.


texts: Scheduled reports in ArcGIS Enterprise (emailed; condition-based)

texts:
ArcGIS Enterprise 10.7.1; Oracle 18c

I want to set up scheduled tabular reports based on Oracle GDB data:

The reports would run weekly on a server (not an individual's computer)
Emailed to a list of recipients (Microsoft Exchange; no Office 365/Power Automate)
The report filetype can be anything; the important part is that people get notified of records in a resultset by email — so that they can follow-up as needed.
The reports would only be emailed out if the resultset/count meets certain conditions:

Example: email report A if the result set has any records.
Or email report B if the result set does NOT have any records.

if recordcount in resultset > 0 then email user@municipality.ca

--selects construction projects that are out of date
select * from ROAD_CONSTRUCTION
where status='FUTURE' and YEAR  <  extract(year from sysdate)

Question:
As a GIS administrator, are there any mechanisms in ArcGIS Enterprise  10.7.1 that would allow me to create and manage scheduled reports?


texts: Find all side-roads of a road via Overpass

texts:
I am curious about how roads spread out – in particular how far you can get, if you start at one road and are only allowed to use that road and all its side-roads.

Right now I use Overpass-Turbo with following query to e.g. find "Frankfurter Allee" in Berlin:

// gather results
(
  node["name"="Frankfurter Allee"](52.451172, 13.291531,52.587381, 13.518124);
  way["name"="Frankfurter Allee"](52.451172, 13.291531,52.587381, 13.518124);
  relation["name"="Frankfurter Allee"](52.451172, 13.291531,52.587381, 13.518124);
);

// print results
out body;
>;
out skel qt;

So what I'd like to do is query for a road by name (like above), but receive all connected side-roads.

Is something like that possible? If so, how?


texts: Best design for Open Source Python / PostGIS prototype

texts:
I am writing a data intensive web application that is delivered through apache. My question is about how to best arrange processing given that there are multiple options.

I have at my disposal OpenLayers/JQuery/Javascript, PostGIS/Postgresql (with pgsql), python/psycopg2, php.

The database contains about 3 million rows and the prototype currently runs as follows:

User clicks on a point on the OpenLayers window
The coordinate is sent as an AJAX request through to a python function  on the server
Currently my application is stateless
Python's psycopg2 is used to call a pgsql stored procedure and a largish set of WKT values  (and a data field) are returned back to the python module
The data field is used to categorize the WKT records in python as follows:  all WKT values are categorized into one of 5 groups. About 1% of the WKT values are actually modified.
The five sets/groups of WKT are buffered to create five distinct polygons.  I currently call a stored procedure in the database to do this. This in turn just uses ST_BUFFER. (I've considered using Shapely but am not sure there will be a performance advantage since the GEOS library is used in either case...)
Finally the 5 WKT text values are wrapped up in a JSON string and sent back to  OpenLayers for rendering as five layers.

I am finding that the bottlenecks are the initial spatial search and the final buffering stage.

I guess the Question is:

Is there a better way to arrange things? For example, should ALL of the data processing be done in PostgreSQL (e.g. with cursors) and is would this be a good thing in terms of maintenance and performance?   Would it be better to use a tile server to avoid passing long WKT strings over to the web client?  How would you address it?


texts: What does buffer mean in QGIS Voronoi dialog?

texts: Could you explain what "buffer" does in the Voronoi dialog i.e. what does the number stand for?


texts: Is this distribution of GPS positions caused by multipath distortion?

texts:
I have some GPS points that are close to where they should be, but not nearly close enough (horizontal precision is around 10 meters). Looking at the positions which comprise these points (there are 10-15 positions per point), the positions are in a nearly straight line. Normally, I would expect to see a more or less round grouping of positions with the point somewhere in the middle, which I do see for most points in this dataset.

The points were collected on a 10-20 foot wide road through a forest, with handheld Trimble units (the antenna was in the handheld, not on a pole), and I was told that in most cases the points were taken by someone who had their back against a tree. The area is in the middle of a bayou, there are no nearby overhead obstructions or buildings nearby, overhead or otherwise. I've also differentially corrected the locations with base data from a nearby CORS base station.

One obvious explanation is that the person was actually walking and not standing still when taking these points, but I've been told that wasn't the case. My only other guess as to what's causing this pattern of positions is multipath distortion. I don't think these trees are big or dense enough for that, and the fact that the positions are in a nearly straight line seems too convenient.

So does this look like multipath distortion, or something else? The red point in this screenshot is the worst of the bunch, but there are others with the same problem. In this spot, I would have expected the red point to be pretty close to the southernmost (bottom-most in this screenshot) GPS position (the small black dots). As you can see, points 740, 736, and 737 have pretty tight clusters of positions, the blue points are ones with much higher precision. 738 and to a lesser extent 739 have a bit of the straight line problem but not as bad as the red point.


texts: How do I call ogr2ogr from a C# application?

texts:
I have an ogr2ogr (GDAL) command that I run to export points from a database to an XML file.

ogr2ogr-f GPX c:\temp\points.gpx PG:"host=localhost port=5432
dbname=SpatialPlayground schemas=public user=postgres password=password"
-sql "SELECT * FROM PointsForGpxExport"

How do I run this from a C# console application and get the resulting file (or a stream) so I can perform additional processing?

One thought is to start a new process from within the console app then grab the exported file and continue processing. Is this the best approach?


texts: Converting ArcGIS layer style to Styled Layer Descriptor (SLD) file?

texts:
I am putting together a GIS server using GeoServer and need some help with importing my layer styles from ArcGIS 10 to SLD format.

I have one layer in ArcMap properly symbolized and labeled, how can I get this layers symbolization and labeling into a SLD file?


texts: Draw a Line between two points and zoom on the new layer - PyQGIS

texts:
In my QGIS3 Plug-in I have this code to draw a line between two points. The code works, but after running it, QGIS zooms somewhere in the middle of Pacific Ocean, when I'd like to zoom on the new layer ("line"). The problem occurs at the last line of this code when I call addMapLayers().

def drawLine(self, line_start, line_end):
  start_point = QgsPoint(line_start.asPoint())
  end_point = QgsPoint(line_end.asPoint())
  v_layer = QgsVectorLayer("LineString", "line", "memory")
  pr = v_layer.dataProvider()
  seg = QgsFeature()
  seg.setGeometry(QgsGeometry.fromPolyline([start_point, end_point]))
  pr.addFeatures( [ seg ] )
  v_layer.updateExtents()
  crs = v_layer.crs()
  crs.createFromId(4326)
  v_layer.setCrs(crs)
  QgsProject.instance().addMapLayers([v_layer])

I tried to add iface.zoomToActiveLayer() but nothing changed. How can I zoom on this new layer avoiding the Pacific Ocean issue?


texts: Mapbasic v15 open wor within existing mapview

texts:
Using Run Application "G:\System\SSA\wor\Current\Townplanning_2018_v1\Zoning.wor" that's called from within a mapbasic program I want to load the workspace in the current map view rather than in a new one.

If I have

Dim nWindowID As Integer
nWindowID = FrontWindow()

How can I tell the run application to open in this window? This has to be in v15 32bit and not 64bit.


texts: Error in image.visualize: Creating NDWI animation series for Landsat 8

texts:
I'm currently working on QGIS with GEE plugin, making NDWI animation series for Landsat 8. Im getting a error in image.visualize: Expected a string or list of strings for field 'bands'... at the end of the script. I would some advises how to solve this or how to end the script and get the animation link in GEE.
L8_comp_1 = ee.ImageCollection('LANDSAT/LC08/C01/T1').filterDate('2013-02-01', '2014-12-01').select(['B3','B5']).sort('system:time_start').filter(ee.Filter.lt('CLOUD_COVER',Cloud_cover))

def function(img):
  doy = ee.Date(img.get('system:time_start')).getRelative('day', 'year')
  return img.set('doy', doy)

col = L8_comp_1.map(function)
distinctDOY = col.filterDate('2013-01-01', '2014-01-01')
filter = ee.Filter.equals(leftField = 'doy', rightField = 'doy')
join = ee.Join.saveAll('doy_matches')
joincol = ee.ImageCollection(join.apply(distinctDOY, col, filter))

def addNDWI(img):
    NDWI = img.normalizedDifference(['B3','B5']).rename('NDWI').copyProperties(img,['bands'])
    return img.addBands(NDWI)

withNDWI = joincol.map(addNDWI).select('NDWI')

def function2(img):
    doycol = ee.ImageCollection.fromImages(img.get('doy_matches'))
    return doycol.reduce(ee.Reducer.median())

comp = withNDWI.map(function2)

def function3(img):
  return img.visualize({'min': 0.5,'max': 1,'palette': ['00FFFF', '0000FF']})

rgbVis = comp.map(function3)

gifParams = {
  'region': region,
  'dimensions': 600,
  'crs': 'EPSG:3857',
  'framesPerSecond': 10
}
print(rgbVis.getVideoThumbURL(gifParams))


texts: Creating choropleth map using OpenLayers?

texts:
Is there any way to generate a chloropleth map using OpenLayers 3?

The map will contain static data that is stocked in a JSON variable.

The map will be published like the picture below:


texts: Dealing with layers that sometimes load and sometimes do not via ArcGIS API for JavaScript?

texts:
I have a map with 3 layers: 2 featureLayers and one ArcGISDynamicMapServiceLayer.

sometimes it loads perfectly, sometimes(50% or 60% of the time) it doesn't.

When the layers do not load, a very cryptic message appears:

Error: Unable to complete operation.
...,"https:"));k=h.content;var x=h.url,y=s&&h.form,F=D; v=t.isDefined(v)?v:F.useCors...

At first I believed this error is caused by network latency. That is, sometimes the network is too slow and when the js map performs some operation, the data hasn't arrived yet and the error is thrown, however, i've noticed that when data does not load a call to /ArcGIS/rest/info?f=json fails.

What will be the best way to handle this situation?

Code:

require(["esri/map", "esri/geometry/Extent", "esri/layers/FeatureLayer", "esri/layers/ArcGISDynamicMapServiceLayer"],
function(Map, Extent, FeatureLayer, ArcGISDynamicMapServiceLayer) {
    var bounds = new Extent({
        "xmin" : -79.50,
        "ymin" : -4.60,
        "xmax" : -66.50,
        "ymax" : 12.8,
        "spatialReference" : {
            "wkid" : 4686
        }
    });

    var theMap = new Map(idDiv, {
        extent : bounds,
        slider : false,
        logo : false,
        showAttribution : false
    });

    //window.setTimeout(function() {
        var urlFeatureServer = "xxx/FeatureServer";

        var fl = new FeatureLayer(urlFeatureServer+"/0", {
            id: "fl1",
            mode: FeatureLayer.MODE_SNAPSHOT,
            outFields:["*"]
        });

        var  myId = xxxx;
        fl.setDefinitionExpression("OBJECTID=" + myId);

        theMap.addLayer(fl);
    //}, 1500);

    theMap.addLayer(new ArcGISDynamicMapServiceLayer("xxx/MapServer", {id: "myBaseLayer"}));

    fl = new FeatureLayer(urlFeatureServer+"/1", {
        id: "fl2",
        mode: FeatureLayer.MODE_SNAPSHOT,
        outFields:["*"]
    });
    fl.setDefinitionExpression("id_xxx =" + myId);

    theMap.addLayer(fl);
});


texts: Ordinary kriging tools available in ArcGIS?

texts:
I have been using Geosoft's Target to interpolate X, Y coordinates then use them to create grids using ordinary kriging.

I am now working with ArcGIS which I am new to and can't seem to be able to interpolate my coordinates and create the grids. Any ideas?


texts: Convert Geometry.boundingbox() coordinates to geometry PyQGIS

texts:
I have a shapefile of 50 line segments and I want to create a separate shapefile containing a bounding box (envelope) for each segment.
for feature in roads.getFeatures():
    #we need to get an envelope around each feature
    inAttr = feature.attributes() # Input attributes
    inGeom = feature.geometry() # Input geometry
    bf_inGeom = inGeom.buffer(25,1).boundingBox()
    xmax = bf_inGeom.xMaximum()
    xmin = bf_inGeom.xMinimum()
    ymax = bf_inGeom.yMaximum()
    ymin = bf_inGeom.yMinimum()
    cordlist = [QgsPoint(xmax, ymin), QgsPoint(xmax, ymax),
                QgsPoint(xmin, ymax), QgsPoint(xmin,ymin)]
    gPolygon = QgsGeometry.fromPolygon(cordlist)

I am just learning the API and while I've successfully loaded the shapefiles, I can't seem to make and display a bounding box polygon.
I'm intermediate on Python. I'm pouring over the docs, but any tips about learning this API would help answer this question too.
How can I convert QgsRectangle to QgsGeometry?
I'm on QGIS 2.14.11 using Debian 9 Stretch.


texts:
I have to connect points to a polygon. I tried to use "points to path" and "point connector" plugins in QGIS, but neither of them are working.
For testing, I just drew 4 points and used "points to path" and "point connector". The created layer does not contain any features.
I did the same with importing 4 points from a txt file, the result is the same.
Is there some trick on how to use those plugins?


texts: Why is the city of San Francsico floating over Point Richmond?

texts:
I have a shapefile that is supposed to be the city of San Francisco. I got it from the good folks at Transbase. On their tiles, it seems to align with the city itself, but when I uploaded it to CartoDB, it wound up hovering about ten miles north of where it belongs:

You can see the map on CartoDB at http://cdb.io/1Lv2U4o.

I did use QGIS to merge the census tracts into a single polygon, but I don't see how that would shove everything north.

I'm not even sure how to troubleshoot this! I assume it's a projection issue, however?

If I do SELECT ST_srid(the_geom) AS geom_srid, ST_srid(the_geom_webmercator) AS webmerc_srid FROM san_francisco_california I get two different values:

geom_srid: 4326
webmerc_srid: 3857

though I see those same two values in other, perfectly functional datasets.


texts: Graphic clipping of Data Driven Pages to index feature that defines each page?

texts: I wonder if there is a way to cut all the data in each polygon separately when using data driven pages?


texts:
I am facing issue while loading the KML file on the mapbox map. For loading KML, I am using leaflet-omnivore https://github.com/mapbox/leaflet-omnivore library.

For reference purpose I am using this example https://www.mapbox.com/mapbox.js/example/v1.0.0/omnivore-kml/

It is working fine for the KML who has the single Placemark tag. But is not working for the multiple Placemark tags.

Can any please let me know how can I load the KML with multiple placemarks?

Or is there any other issue.

Here is my sample code

var runLayer = omnivore.kml('kmlfiles/sample.kml', null, customLayer)
            .on('ready', function() {
                runLayer.eachLayer(function(layer) {
                    layer.bindPopup(layer.feature.properties.description);
                });
            })
            .on('error', function(e) {
                console.log("Error loading KML: ", e);
            })
            .addTo(map);


texts: Prevent row deletion in versioned tables using database privileges

texts:
I have versioned tables in an Oracle geodatabase (versioned with the option to move edits to base).
I would like to prevent users from deleting rows from the tables when using ArcGIS Desktop.
(With that said, the users still need to be able to select, update, and insert rows into the tables.)

The reason for preventing the deletion of rows is:
The asset records in the tables need to be preserved/decommissioned
so that they can be integrated to a work order management
system. We cannot allow rows to be deleted from the GIS,
because the GIS records would become out-of-sync with the records in
the work order management system (in which deletes are not possible, and
historical records are preserved).

ArcGIS functionality:
It is my understanding that preventing deletes is possible with out-of-the-box ArcGIS functionality on unversioned
tables:

Grant and revoke dataset privileges:
If the dataset is not registered as versioned, you can grant and
revoke the update, insert, and delete privileges individually using
the Privileges dialog box. For example, you can grant a user select
and update privileges, which allows the user to connect to the dataset
and alter existing features but does not allow the user to add new
features or delete existing features.
If the dataset is registered as
versioned, the privileges that allow a user to modify a dataset
(update, insert, and delete) must be granted and revoked as a group.

Question:
I understand that preventing deletes is not possible on versioned tables using out-of-the-box functionality in ArcGIS.
However, is there a reason why I shouldn't prevent deletes via database privileges?
REVOKE DELETE ON ROADS FROM ROLE_EDITORS_GROUP_A


texts:
Even though this is a single line, it is stored as a two-part polyline. I can split it and rejoin at the break, but whatever I do it won't join the two into one line. (This is most obvious if you turn on labels.) What's going on here?

http://pastelink.me/dl/a8ca3e


texts: Convert .mxd with embedded VBA command button to ArcMap 10

texts:
I'm trying to convert an .mxd file from an older version to ArcMap 10.  The old version had two simple VBA scripts attached to toolbar buttons.  I am trying to convert this .mxd to ArcMap 10.  For deployment reasons I would like to maintain everything in a single file.  I would also like to continue using toolbar buttons to run the scripts.

I thought all I would have to do is convert the VBA code to Python.  It's apparently not that simple.  Can someone please explain whether the following is possible in ArcMap 10 (seems like it should be since it was simple in previous versions):

single .mxd (no extra .tbx or .py files)
two custom toolbar buttons
click button1 and it runs some_python_function_1
click button2 and it runs some_python_function_2


texts: Get Modified Feature from Modify Interaction

texts:
THIS QUESTION IS STILL OPEN AND OF INTEREST

Openlayers: 4.6.5

I have a Modify Interaction and two active Select Interactions. The Modify Interaction is instantiated like this:

   var interaction = new ol.interaction.Modify({
     source: myLayer.getSource()
   });

When I click the edge of a polygon feature, a vertex of this polygon is highlighted - when I draw, the vertex is moved and when I release, the MODIFYEND is triggered. There is obviously exactly one feature affected.

After Modifying a Feature Geometry with the Modify Interaction, I expect the triggered modifyend event to provide the feature or the geometry that has been affected by the modification. However, it only provides me an ol.Collection containing ALL Features of the interaction's source.

How can I obtain the feature then?

interaction.on(ol.interaction.ModifyEventType.MODIFYEND,function(event) {

    // event does not have a feature property

});

Edit: I totally want to avoid this:
I found this as a working solution, but I am not happy about it and I don't think this is the appropriate way:

interaction.on(ol.interaction.ModifyEventType.MODIFYEND,function(event) {

       var featureModified = interaction.dragSegments_[0][0].feature;

    });</s>


texts: How to paint the color of the road between points geolocation?

texts: We have sensors set on the map, each of them has  coordinates, we made a rectangle using Graham scan and for the moment we need to paint up all roads which are included into this rectangle. Can you advise which method we should use to make it? We have not found any example in ARCGIS API (version 3.18).


texts: Is there a way in QGIS to cluster rather than overlap SVG icons from all layers

texts:
I have many layers in QGIS containing points. All layers were added as delimited/csv.

PROBLEM: These are points from different layers, I'd like them clustered around the area instead:


texts: Automatically Generate Attribute Table in Qgis

texts: I am working with historical parcel data in Qgis and have the following fields in my attribute table, (name, grantee date, grantor date, parcel ID, acreage, previous owner). Unfortunately, I have around 10-15 already created parcels that does not have these fields in the table. Is there an easy way to have those fields automatically generated in those attribute tables or will I have to manually enter them?


texts: QGIS Bonn - Trying to use query builder to exclude multiple test strings

texts: I am entering "textString" != 'Hilltop' OR "textString" != 'Amberley' in to the Provider specific filter expression in order to EXCLUDE these two names in OS Cartographic text. The query correctly shows 1353 rows, but these no longer show on my OS map. What am I doing wrong?


texts: Using matplotlib to create and export graphs automatically from ArcMap

texts:
I am going to start a project soon where, after the data has been processed in various ways, an output jpg showing a graph of the data should be produced. The workflow is as follows:

1. A feature class with several fields containing various values is created.
2. A radar/bar chart based on several fields within the feature class should be produced
3. Colour the chart according to whether the field value is within a certain threshold e.g. if the threshold for the field CHEM1 is 4, a feature with CHEM1 = 2 will be coloured red, while a feature with CHEM1 = 6 will be yellow.
4. Each unique feature within the feature class should have its own graph i.e. 10 features = 10 graphs
5. The graphs should be saved as jpgs
6. The jpgs should be added as attachments to the geodatabase the feature class is in - not needed at this point but should be kept in mind

I will be automating this process after the initial run, as the new data will be processed monthly and the jpgs generated automatically, so a python script will definitely used. Before I begin though, I am unsure as to which graphing method I should use. I have researched the following options for my approach:

1. Python + ArcPy graphing
2. Python + matplotlib
3. Python + Excel

My pros and cons for each option:

I am used to ArcMap's graphing options, and have been using arcpy
and ModelBuilder daily for the past year, so it will be relatively
simple for me to use the inbuilt arcpy functionality.
I have never used matplotlib before, but my brief research indicates that it is quite simple to use for basic charts, and it doesn't seem to be that much more difficult to create more advanced charts using it. Also, since it's included with ArcGIS, it seems like a logical choice.
I have worked heavily with Excel for many years, and am very familiar with creating charts in it. I also did VBA, and it is my understanding that after looking at the VBA code for an action, one can translate that into Python quite quickly.

I am leaning towards the third option, since I feel most comfortable with that combination. However, since this process will have to be automated completely and will eventually run on the server, I do not like the idea of having to install Excel on the server so that this process can be run. I also fear that all the arcpy code will slow the process down (from a combination of not having enough experience with the graphing methods, not being able to optimise my code, slow computers in general, network issues), and that even though matplotlib appears to be the simplest/fastest approach, I will end up taking longer writing the code because I have not used it before.

My question is not which combination works the best overall (because that would be subjective/too localised), but which option would be the best combination to achieve the outcome for this particular task?


texts: How would one geocode historical US census records?

texts:
I want to geocode census records from the U.S.  Here is one example from 1870.  This particular census page is from the 6th election district, 8th ward of New York County.  Here's one from the 1910 Census.  This one mentions the county, township, supervisor district, and enumeration district.

I can geocode the county, state.  I'm wondering if there are any resources which include latitude/longitude for districts, wards, supervisor districts, and enumeration districts from the past to help me geocode census records.  I did find NHGIS, but their data is quite sparse for anything more fine grained than the county level before 1940.


texts: Identifying Open Areas in OpenStreetMap?

texts:
I am currently creating a system using OpenStreetMap (OSM) and I need to identify open areas in the map.

I know there are existing open areas on OSM but is there a way to identify more open areas that are not pre-identified using OSM data?
Or can I do this using a different approach?

By open areas I mean, open fields, parks, or any areas away from structures like buildings


texts: Python code does not work with RasterCalculator after exported from ModelBuilder?

texts:
I have a model authored in ModelBuilder that uses "RasterCalculator" to process input raster datasets.

 arcpy.gp.RasterCalculator_sa("OutRas = Con(\"%STATIC_forestland%\"  ==  0, -255, Con(\"%mapvictoria20121210%\" >= 115, -255, Con((\"%mapvictoria20121210%\" >100) & (\"%mapvictoria20121210%\" < 115), 100, Con(\"%mapvictoria20121210%\" <= -15, -255, Con((\"%mapvictoria20121210%\" < 0) & (\"%mapvictoria20121210%\" > -15), 0, \"%mapvictoria20121210%\")))))", mapvic1)

It works well as a model. However, when exported to a Python script and saved in the same Toolbox, the script does not work.

The error message is as below

RuntimeError: ERROR 000732: Input Raster: Dataset %mapvictoria20121210% does not exist or is not supported

Failed to execute (RasterCalculator).

I am sure the input raster are there. Can anyone give some help? Is it to do with raster naming issues?


texts: Using pyqgis logic to stop process if output has no geometry

texts:
I have a pyqgis script which customises the functionality of the grass algorithm r.lake.coords:

##Seed Points Lake=name
##DEM_input=raster
##Seed_vector_points=vector point
##Overspill_vector_boundary=vector polygon
##output=folder

# Import classes
from qgis.core import QgsRasterLayer
from qgis.utils import iface
from PyQt4.QtCore import QFileInfo
import processing

# Define layers
point_layer = processing.getObject(Seed_vector_points)
layer = processing.getObject(DEM_input)
boundary_layer = processing.getObject(Overspill_vector_boundary)
fileName = layer.source()
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)

# Get extent of the region
extent = iface.mapCanvas().extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

for feat in point_layer.getFeatures():
    geom = feat.geometry()
    attr = feat.attribute("Z_ADD_1")
    i = 0
    while i < 5:
        output_1 = processing.runalg('grass:r.lake.coords', layer, "%f"%(attr), "%f,%f"% (geom.asPoint().x(), geom.asPoint().y()) , False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), 0.00, None)
        output_2 = processing.runalg('saga:rastercalculator', output_1['lake'], '','a>0', False, 3, None)
        output_3 = processing.runalg('gdalogr:polygonize', output_2['RESULT'], "DN", None)
        output_4 = processing.runalg('qgis:dissolve',output_3['OUTPUT'], True, "DN", None)
        output_5 = processing.runalg("qgis:extractbylocation", output_4['OUTPUT'], boundary_layer, u'within', output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + "_%.2f"%(attr) + 'm.shp')
        i += 1
        attr += 0.5

My problem is that the final algorithm which the script calls upon produces additional blank shapefiles where there is no geometry, this is technically correct if the outputs from output_4 are not within the overspill boundary.

However, I would ideally like to have the script set up using some sort of python logic where it avoids producing final shapefiles with no geometry.


texts:
I have "unsorted" points from various sources  - steming from a hiking track / path.

How can I create one line from the closest points (...following the track). I´ve already imported all points to one layer, but I´m not able to create a proper line.


texts:
I got a table with the name Predios

like this:

OBJECTID   PRM_IDENTIFICA

1          4030200606130

2          403020060615

3          40302006011700

Can I use search cursor to find the records when the PRM_IDENTIFICA len its <> 14.???

Right now i got:

import arcpy

tabla="C:/Users/jaguero/Desktop/validacion tabla/Valid.gdb/Predios"

rows=arcpy.SearchCursor(tabla)

row=rows.next()

while row:

print len(row.PRM_IDENTIFICA)

row=rows.next()

giving me this:

13

12

14

What if I want only as a result the records when the len is <> from 14 ?


texts: How to create a smallest quadrilateral concave polygon around GPS points?

texts:
I have this GPS track in a metric grid (Symbology by elevation, if anyone wonders):

I need to define it's general area - for registrations and permits, with just 4 points. Therefore, I assume that the corners of a circumscribed quadrilateral polygon (4 lines) will do, as that will be the smallest possible polygon area-wise.

I don't know which Arcmap tools I can use to do it, though I imagine the result should look something like this:

With a polygon I can easily know the corners' coordinates, but if there's another way to extract them, I'm all ears.

Thanks in advance.

Edit:
This is the resulting convex hull tool: