texts:
I downloaded a .tif raster file from opentopography.org. It's an Alos (AW3D30
) topography file showing mountainous terrain in Sweden.
On the very right the file contains terrain, and on the left a lake surface.
The lake surfaces have values of -9999 for each cell.

Is there some way to interpolate these -9999 values into values which would somehow correspond to the lake surface elevation? For example, the elevation of the nearby shore/coast is around 300 something meters.
Is there a way to replace all the -9999 values with 300 values in order to get a somewhat acceptable transition from shore/coast to lake surface?

The mentioned file can be downloaded from here:
https://www.dropbox.com/s/kai6202uewas47j/no_value.tif?dl=0


texts: OpenLayers 3: Dynamically or batch add layers from user folder

texts:
I would like to dynamically load all layers (KML, for instance) from a user's folder to a group map. This likely requires an iteration that searches the folder and adds the layers to the users' map.

I've been to here and here, but neither of these discuss iterative  layers from a given location, nor how to parse n layer names, presumably by Document.getElementById().

Whatever the route taken, list needs to pass to a function I have in a table for creating toggle buttons on each layer.  This code works, just curious as to the iteration's effects.

function toggleControl(element){
     if(element.value == "dynamicallyaddedlayer1"){
      dynamicallyaddedlayer1.setVisibility(element.checked);
     }
     if(element.value == "dynamicallyaddedlayer2"){
      dynamicallyaddedlayer2.setVisibility(element.checked);
     }}
//...... and so on for as many n layers present

    <fieldset data-role="controlgroup" style="margin-top:0px">
        <input type="checkbox" name=“dynamiclayer1" value="dynamiclayer1btn" id="dynamiclayer1btn" class=".ui-corner-all" onclick="toggleControl(this);" />
        Name-of-Dynamically-Added-Layer-#1
//... and so on for as many n layers present
        </fieldset>


texts:
I have two airborne lidar datasets and, for purposes of comparing change between the two, would like to standardize the vertical and horizontal coordinate information for each. The plan is to transform the horizontal/vertical system of the 1st dataset to match that of the 2nd dataset. Here are systems for each:

1st dataset:
Year - 2010-2011
Horizontal - NAD83 UTM 4N
Vertical - GRS80 Ellipsoid

2nd dataset:
Year - 2017
Horizontal - ITRF00 projected into UTM 4N
Vertical - NAVD88 (Geoid12A)

I simply am not well versed in geodesy enough to know how to proceed in completing this transformation. I have access to Blue Marble Geographic Calculator and NOAA's Vertical Datum Transformation tool.

How should I proceed here?

I should clarify my issue is not with the usability of any of the software packages. Rather, I'm confused about the actual best practices for the conversions to make (i.e. converting from lidar shot at an ellipsoid to those referenced to NAVD88 with a particular geoid12A).


texts:
I have a shapefile that contains multiple polygons. I want to convert it to a kml but also keep the face color. From other answers here, I have found that I can add a field called OGR_STYLE to the shapefile and specify the fill parameters as a string:
OGR_STYLE: BRUSH(fc:#0000FF);PEN(c:#000000)

I then use ogr2ogr to convert the shapefile to a kml using the kml and libkml drivers:
ogr2ogr -f kml -f libkml output.kml input.shp

This appears to run okay, however my kml does not have a blue face color as expected. The outline is black, as expected, but there is no face color.
Am I overlooking something?


texts: Adding vector layer and newly created layer always in Top- Down Order Panel using Pyqgis 3

texts:
I have two layers in my project

Point Layer
Building Layer

I recently created Voronoi Polygon on Point Layer say (vor-points layer)
But as soon this Happen, the newly created vor-points layer appears on the top of Layer Panel List.

How should i make sure that my Layers (by uploading or newly created ones ) always appear in top to down order and my canvas should show all the layers simultaneously ?


texts: Third party base layer other than Google map

texts: Can we use third party base map (images tiles) for displaying vector data?, like we use google map as base layer for displaying vector data in OpenLayers application.


texts:
I have several polygons from a drawing tool (eDraw) exported from the Web AppBuilder dev edition into a directory. They are exported in JSON format with a EPSG:3857 projection.

Through a script, I want to transform these JSON files to shapefiles, and store them in a different directory. This works perfectly, however there are three files that keep getting the error that they are out of bounds.

  arcgisscripting.ExecuteError: ERROR 001558: Error parsing json file
  ''. The coordinates or measures are out of bounds.

I have tried to reproject them to no avail, I have tried to trace the failing polygons in the WAB and re-export them, but no luck.
If I compare these three files to other files that do work I see no significant difference. Also, the two example files are right next to each other on the map, and I have no problem uploading the 'broken' JSON files to the drawing tool and display them in the right location on the map, so I do not understand how one can be passed with JSONToFeatures, and the other is out of bounds.

Working file:

{"features":[{"geometry":{"rings":[[[-8734229.349438265,-1044267.8697553115],[-8733964.208496368,-1044023.0323991263],[-8733867.46788243,-1044123.9531630479],[-8733481.69975528,-1043799.095792834],[-8733485.282740982,-1044816.6637318594],[-8733914.046696614,-1044820.8438817654],[-8734001.232681995,-1044775.4593962248],[-8733921.212668015,-1044651.249225272],[-8734229.349438265,-1044267.8697553115]]],"spatialReference":{"wkid":102100}},"attributes":{"name":"Scallops Peru 1 SAC ","description":"","objectid":1},"symbol":{"color":[255,199,198,99],"outline":{"color":[207,101,99,255],"width":0.75,"type":"esriSLS","style":"esriSLSSolid"},"type":"esriSFS","style":"esriSFSSolid"}}],"displayFieldName":"","fieldAliases":{},"spatialReference":{"wkid":102100},"fields":[{"name":"objectid","alias":"objectid","type":"esriFieldTypeOID"},{"name":"name","alias":"name","type":"esriFieldTypeString"},{"name":"description","alias":"description","type":"esriFieldTypeString"}]}

Not working file:

{"features":[{"geometry":{"rings":[[[-48806317.836461626,-1048103.3855233326],[-48806231.844804816,-1048077.1102948618],[-48805966.70386297,-1048325.5306367677],[-48806058.66716262,-1048642.0277069844],[-48806058.66716262,-1048788.9301207075],[-48806317.836461626,-1048790.1244492745],[-48806317.836461626,-1048103.3855233326]]],"spatialReference":{"wkid":102100}},"attributes":{"name":"Cultivos Casma SAC","description":"","objectid":1},"symbol":{"color":[255,199,198,99],"outline":{"color":[207,101,99,255],"width":0.75,"type":"esriSLS","style":"esriSLSSolid"},"type":"esriSFS","style":"esriSFSSolid"}},"displayFieldName":"","fieldAliases":{},"spatialReference":{"wkid":102100},"fields":[{"name":"objectid","alias":"objectid","type":"esriFieldTypeOID"},{"name":"name","alias":"name","type":"esriFieldTypeString"},{"name":"description","alias":"description","type":"esriFieldTypeString"}]}

This is part of the script:

import arcpy

import os

import json, sys

dirOfJSON = '<path>'

outputDIrectory = '<path>'

for f in arcpy.ListFiles('*.json'):
    j = os.path.join(dirOfJSON, f)
    outputShp = os.path.join(outputDIrectory, os.path.splitext(f)[0] + '.shp')
    print("Processing : {}".format(outputShp))
    arcpy.JSONToFeatures_conversion(j, outputShp)


texts: Geoprocessing Tools on ArcGIS 10.2 and Visual Studio 2012 -issues

texts:
I am migrating one of the existing project (custom Geo-processing tool- class library project) source code which runs on ArcGIS 10 SP4 and Visual studio 2010. This is working fine. I am could run the project and able to add this tool in ArcTool box. Publish this as GP service and consumed in application.

I am moving this code to new machine which Windows 2012 server and  ArcGIS 10.2 and Visual studio 2012. The source code is fetched from TFS and and code is built successfully after updating ArcGIS libraries to 10.2 When I run the project in VS2012 - class library (starting - ArcMap.exe with specific mxd document  as input) break points are not hit.

After doing some research I have found that ESRI Link here I have to update ArcMap.exe.config file. After updating config file to use Framework 4.0 break points are set but not getting hit. Whenever I run the class library project MXD is opens, but the break points are not getting hit and not able to debug the code. Project properties in VS2012 is same as I used in VS2010 (checked Register COM Exception) and enabled debugging.

Question:
a) Why VS2012 is not hitting breaking points?
b) Dll is successfully built, but I am not seeing this custom tool box in "Add tool " box in ArcTool box.

I am not sure what am missing or some changes in ArcGIS 10.2. Please help.


texts: Get Isochrones for Public Transport in QGIS

texts:
I'm using QGIS and trying to get public transportation isochrones in Montreal for 20, 40, and 60 minutes.

I was able to use the ORS plugin to do the same for cars, walking and cycling in the area. I tried to use the TravelTime plugin for public transit but the results came out very inaccurate for the Montreal area (it showed a very small area accessible by public transportation within 60 minutes and I live here so I know its way off)

I put together .csv files for all of the bus and train transit lines in the area and brought them into QGIS as delimited text layers then merged them using MMQGIS plugin

Basically I have a sample size of roughly 4000 postal codes from residents in the greater Montreal area and I want to see how many of them live within 20, 40 and 60 minutes from the Concordia University campuses (which I have points for).


texts: Line and Polygons Intersect Results QGIS 3.0

texts:
Dark Green - Line

Light Green - Buffer created from Line

Pink Polygons - Parts of the black polygons that fall within the buffer zone

I have this line going though a bunch of polygons.
The result I need is the overall footage of the Line that potentially intersects the black polygons. The buffer is what defines the potential intersects.

I have come up with a few methods, but they are very complicated with many steps.
I was hoping someone knows a simple or simpler solution to do this.
My current solution goes like this Buffer, Intersect, Dissolve, Multipart to Singlepart, Bounding Boxes, Expand Bounding Boxes, Dissolve, Intersect, Multipart to Singlepart, Add $length attribute.


texts: How to use Hub distance in a processing model?

texts:
I am trying to make line layer (hub distance layer ) between two point layer in QGIS 2.16.

It works fine when I make layers step by step but, I want to make MODEL in processing modeler and automaticly make final layer.

I got an error message like this.

This is a result of what I want.


texts: Data interpreted incorrectly when opening NetCDF file in ArcMap

texts:
I'm trying to open a NetCDF file with some current data on a polar stereographic projection in ArcMap, but the data are not displayed correctly. The only obvious problem I can see is that ArcMap seems to be interepreting some of the global attributes incorrectly. For example, in ncdump it says

:false_easting = 2384800.f ;

Whereas in ArcMap it says

False_Easting 1.890518526342814e-273

Is there anything I can do to ensure the attributes are read correctly? I thought the point of a format like NetCDF was that it would specify both the data and the data type in a standard way.


texts: Trying to label graphic using feature attribute (ArcGIS Javascript 3.7 API)

texts:
I have a Dialog DIV containing an HTML list of buttons with a value for each school district. I want the user to be able to click on a button and have a polygon graphic of the district placed on the map, zoom to its location, and label the graphic with the district name.

Using A HREF (javascript:findSchoolDist('district_name')) I call the function and zoom to the extent of the returned graphic. To my understanding, when using FindTask, all attributes are returned along with the graphic.

The FindTask successfully finds the feature, returns the graphic, and zooms to its extent. The problem is I cannot get the label to appear. While I am using the Javascript 3.7 API, I am not using the AMD format so the code I am using is borrowed from a 3.5 API sample as they are coded using the legacy format.

For finding the center of the polygon so the app knows where to place the label, the require "esri.tasks.geometry" is included and I have included the Geometry Service URL under the layers as "var geometryService = new esri.tasks.GeometryService("http://gis.abag.ca.gov/arcgis/rest/services/Geometry/GeometryServer");"

The two functions that make up the tool are as follows:

function findSchoolDist(searchText) {
   map.infoWindow.hide();
   map.graphics.clear();

   //Create find task
   var schoolFindTask = new esri.tasks.FindTask("http://gis.abag.ca.gov/arcgis/rest/services/admin/abag_District_School/MapServer");

   //Build school district filter
   var schParams = new esri.tasks.FindParameters();
       schParams.returnGeometry = true;
       schParams.outSpatialReference = map.spatialReference;
       schParams.layerIds = [0];
       schParams.searchFields = ["name"];

   //Set search text and find school district
   schParams.searchText = searchText;
   schoolFindTask.execute(schParams, schShowResults);
}

function schShowResults(results) {
   var schoolSymbol = new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_DIAGONAL_CROSS, new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID, new dojo.Color([255,255,0]), 4), new dojo.Color([0,0,255,0.25]));

   map.graphics.clear();

   dojo.forEach(results, function (result) {
      var schGraphic = result.feature;
      var schAttrib = schGraphic.attributes.name;
      schGraphic.setSymbol(schoolSymbol);
      map.graphics.add(schGraphic);
   });

   var poly = dojo.filter(results, function(r) {
      return r.feature.geometry.type === "polygon";
   });
   if ( poly.length > 0 ) {
      map.setExtent(poly[0].feature.geometry.getExtent().expand(1.5));
   } else {
      alert("No polygon results...");
   }   //WORKS UP TO THIS POINT, FAILS AFTER

   //Add district name to map
   var schFont = new esri.symbol.Font("14pt", esri.symbol.Font.STYLE_NORMAL, esri.symbol.Font.VARIANT_NORMAL, esri.symbol.Font.WEIGHT_BOLD, "Arial");

   var schTextSymbol = new esri.symbol.TextSymbol(schAttrib);
      schTextSymbol.setColor(new dojo.Color([255,0,0]));
      schTextSymbol.setAlign(esri.symbol.TextSymbol.ALIGN_MIDDLE);
      schTextSymbol.setFont(schFont);

   var pt;
      pt = schGraphic.getExtent().getCenter();
   var gra = new esri.Graphic(pt, schTextSymbol);

   map.graphics.add(gra);
}

I am hoping that someone is still using the legacy format who has experience labeling graphics.


texts: Outside ArcMap session, script fails to add/save layer to mxd file

texts:
The following code generates an layer output, but fails to add and save it to the specified mxd.  There are no other layers in the mxd, and is just a 'ghost' / placeholder for running the script outside an ArcMap session.

import arcpy,os

fc_path = "C:\\TEMP\\Footprint.shp"
layer_name = "Footprint"

mxd = arcpy.mapping.MapDocument("C:\\TEMP\\Parser.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
layerfile = os.path.join(arcpy.env.scratchFolder, layer_name + ".lyr")
arcpy.MakeFeatureLayer_management(fc_path, layer_name)
arcpy.SaveToLayerFile_management(layer_name, layerfile, "ABSOLUTE")
add_layer = arcpy.mapping.Layer(layerfile)
arcpy.mapping.AddLayer(df, add_layer)
arcpy.RefreshTOC()
mxd.save()


texts:
I am using a SQL Server 2012 database. I have noticed that geometries and geographies interpret WKTs differently. For instance, if I construct a square polygon, clockwise:

DECLARE @geom Geometry=geometry::STGeomFromText('POLYGON((-30 -45, -30 45, 30 45, 30 -45, -30 -45))', 4326)
DECLARE @geog Geography=geography::STGeomFromText('POLYGON((-30 -45, -30 45, 30 45, 30 -45, -30 -45))', 4326)

The geometry will be interpreted (as intended) as the original square (inwards):

But the geography will be interpreted as the complement of the square (outwards):

Is this standard (and can you find a reference to why this is so)? And if so, how is it possible to construct a geography from a WKT?


texts: OSM available data extract for addresses only?

texts:
Is there a dataextract available for all addresses without the street network?

I would be interested in something like this:

[streetNumber, street, postcode, city, county, state, country, lat, lon]

I would need this because I build a custom geocoder (similar to nominatim) and would need to minimize the dataset size.


texts: Refresh leaflet marker data in place ? (on insert to sql db)

texts:
I have a form that sends data via ajax to be geocoded and inserted into a mysql db.

Once that data is there, can I force the marker layer to refresh so that it displays the new data?

The SQL API for Cartodb includes a layer method .invalidate() which refreshes the data if there has been a change.  I don't care about detecting a change, but I do hope to find a similar refreshing capability in pure leaflet.

The map will not be resized or zoomed when the new data is pulled in, so something like map._onResize() won't have any effect.

There will be reference URLs etc that need to be pulled from the db to add to a popup, so just putting on the marker on the front end while pushing to the db via ajax won't work.


texts: How to use arpy.mapping to change 'out jpeg' name to reflect page name

texts:
I am an experienced Arc user but new to python. I can not figure out where my syntax is wrong. I'm trying to use data driven pages to export to jpeg where the name of the file is the page name. It works fine through the usual export to pdf in File, Export, Pages dialog so I know that the data driven pages are set up correctly. The first line references the field name in the attribute table of the layer file that I have set to be the page name. The script works if I use the page number.

I'm using ArcGIS 10.2.2 for Desktop with a basic licence.

arcpy.mapping.ExportToJPEG(mxd, r"C:\Temp\Linn\Linn_KCH" + str(pageNum) + ".jpg")

This is the script I'm trying to run now.

field_name = "Page_Name"

mxd = arcpy.mapping.MapDocument("CURRENT")
for i in range(1, mxd.dataDrivenPages.pageCount + 1):
mxd.dataDrivenPages.currentPageID = i
row = mxd.dataDrivenPages.pageRow

pageName = mxd.dataDrivenPages.pageRow.Name

print row.getValue(field_name)
arcpy.mapping.ExportToJPEG(mxd, r"C:\Temp\Linn\Linn_KCH" + "." + row.getValue(field_name) + str(pageName) +".jpg")
del mxd


texts: Summary statistics (sum of raster cell values) are equal in R and QGIS, but different in ArcGIS. Why?

texts:
I'm new to GIS and I have a question related to the use of different programs (R, QGIS and ArcGIS) to do summary statistics. I've used R and QGIS to sum the cell values of my raster (using respectively the cellStats function from R (as follows) and the Raster layer statistics from the QGIS toolbox) and in both cases I get the same result.

Following is a link to the data I'm using: https://1drv.ms/u/s!Aod2icrdkQPXgSDjgZg3O_NIO71Y

In R, I do (with raster package):

my_raster <- raster("path_to_my_raster")
sum_raster_cells <-cellStats(my_raster, 'sum')

In QGIS, I do:

However, when I use the arcpy.Statistics_analysis from ArcGIS (as follows), I get a different result. The later actually seems to give a more accurate result than the one that I found with R and QGIS, but I wanted to understand why these results are different. Why isn't the total sum equal when using the same data?

import arcpy
arcpy.env.workspace = "path_where_gdb_is"
arcpy.Statistics_analysis("GPW_ac_ad1_15","output_path",[["SUM_E_ATOT","SUM"]])


texts:
I am trying to smooth 8 days LAI data from MODIS and managed to work around it with several algorithms such as loess, Savitzky–Golay, smooth.spline etc. This is pretty much going well if the level of noise is not really high (please see below plot) for a target feature.

However, for some pixels (i.e. features), the level of noise is so high that I am not sure if the fitted smooth line is realistic or not (please see below plot)?!

I have already filtered low-quality MODIS pixels using QA dataset but still no improvement.

Is there any suggestion for me how to smooth values of such pixels with high levels of noise?


texts: Sending long JSON objects (polygon geometry, table rows) in POST request to Geoprocessing Service?

texts:
I am seeking to understand my options for sending long JSON objects in a request to a Geoprocessing Service.

In particular, I am looking for some example code to illustrate how to include two long JSON objects (one a polygon geometry as a coordinate string, the other a table as rows of ASCII text with the first row being the column headings) into a POST request to my Geoprocessing Service, and then hints as to how best to read them as parameters in the Python script that I've published as my Geoprocessing Service.

At http://help.arcgis.com/en/arcgisserver/10.0/apis/rest/gettingstarted.html it says:

  When using the REST API, you will normally use an HTML GET method in a
  form. When you use GET, the entire request is encoded in the URL. This
  is the preferred method to use whenever possible. However, in this
  mode, a URL is limited to as few as 1024 characters depending on the
  browser. Thus, if you have a long JSON object to include in the
  request, you will need to use POST.

  A second option, when using certain Geometry Service and Geoprocessing
  Service operations, is to continue to use GET and to specify a URL to
  the input JSON object contained in a file on a public server.

I'm familiar with using a GET method when the parameters can all fit into a URL that I can paste into a browser.

I'm also comfortable with the idea of using GET to specify a URL to input JSON objects that are contained in files on a public server.  However, I am prevented from using this method this because my client does not want files placed on a public server.

That leaves me with having two long JSON objects to include in the request and having to use POST which is a method beyond my current level of understanding.

Is any example code readily available?


texts: Deciding whether to Select Transformation in QGIS

texts:
I've installed QGIS and downloaded the sample files (I'm on a Mac). Following the documentation for getting started I double-click on the landcover.img file which is supposed to load in the background. But I get a "Select Transformation for landcover" dialog box which is not mentioned in the documentation and I have no idea what to select.

This is probably a good learning opportunity but I fear that I am on a steep curve.


texts: Implementing WPS client in OpenLayers 3?

texts:
Is it possible (and what is the recommended way) to use OpenLayers 3 to do the work of theOpenLayers.WPSClient?

As shown in the old WPS client example, OpenLayers 2 can be used to operate WPS.

There are OpenLayers 3 examples for using other web services like WFS and WMS, but not WPS (as of 3.2.1).

Is this temporary or by design?


texts:
I am working within a local city government on creating a routing network and started last week with little to no knowledge of pgRouting. I want the routing network to be used by our fire department and public schools. A series of travel-time buffers radiating from each fire-station is my goal. Example: how far can we travel in 1 min? 5min? 10min? from each station.

I found what visual I am looking for via:https://anitagraser.com/2011/02/12/drive-time-isochrones/

So far I have successfully created a database for the routing, created both extensions (PostGIS & pgRouting), imported my shapefile road centerlines successfully into my database, and lastly I ran pgr_createtopology. The results are great! My pgRouting plugin within QGIS works flawlessly. The next step is where I need help!

Should I create nodes? Or, alternatively, can I just use the vertices from the ""_vertices_pgr table generated from function: pgr_createTopology?

The ""_vertices_pgr table contains an "id" column but I can't see any column with any geographic data within the table. Is this why we can't use the vertices table? Is there something I am misunderstanding about creating a node network when I already have vertices created from pgr_createTopology?


texts: Rapid digitizing of Point Data using ArcGIS Desktop?

texts:
I have several big maps with a lot of points that I have to digitize. All these Points have different Values (depth in Meter, double) and will be stored in 1 Shapefile.

Is there a way to speed up my digitization process?
e.g. Create Point with my Mouse in right Hand and directly Type the Value into a NumPad in my left Hand?

Right now I always have to select the Item in the Attribute Table, where i can edit the Value. (Too much clicks)


texts: Improving select and edit features performance of QGIS plugin?

texts:
I wrote a plugin for selecting and instant updating the value of the selected feature in the attribute table. I tested with a big shapefile and to select and edit 5 features takes ~5 seconds. Which isn't really fast.

Can someone look at the piece of code and give me suggestions for how I can do for faster editing?

        lRect = self.canvas.mapRenderer().mapToLayerCoordinates(layer, r)
        layer.select(lRect, False)
        layer.startEditing()
        #for feat in layer.selectedFeatures():
        for feat in layer.selectedFeaturesIterator():
            feat["foo"] = my_value
            layer.updateFeature(feat)
        layer.commitChanges()
        layer.removeSelection()

I could use selectedFeatures() or selectedFeaturesIterator() but I tested it with both and I got nearly the same result (time).


texts: Getting layer from multivalue input in python toolbox

texts:
param5 = arcpy.Parameter(
        name="C",
        displayName="C",
        direction="Input",
        datatype="DEFeatureClass",
        parameterType="Required",
        multiValue=True)

    parameters = [param0, param1, param2, param3, param4, param5]

    return parameters

And now I want to get each of layer from param5. I must do this such as here?

c = parameters[5]
 c1 = c[0].valueAsText
 c2 = c[1].valueAsText

Multivalue parameter is a list? What is the way that get each of leyer? If I would like to do Merge (http://pro.arcgis.com/en/pro-app/tool-reference/data-management/merge.htm). That's enough pass parameters[5] as first argument to Merge_management??


texts: 6S atmospheric correction method for Landsat 8

texts:
I m not clear with 6S atmospheric correction. which image we use to run 6S atmosphere as an input image to calculate coefficient xa,xb and xc
top of atmosphere radiance or top of atmosphere reflectance image?


texts: Bulk load multiple shapefiles into PostGIS

texts:
I've got over a hundred data layers in shapefile format that I want to upload into a PostGIS database.  They are all in the same projection, but represent different data layers so they do not have the same schema.

What is the most time-effective way to bulk convert all of these files into my PostGIS database (on windows 7 OS)?

UPDATE: Bulk importing in the pgAdmin3 "PostGIS shapefile uploader" (as noted below) is now available by default with PostGIS 2.0.


texts: How to turn off OpenLayers 4 smooth zooming

texts:
After upgrading from OpenLayers3 to OpenLayers4, the zooming has changed. Now the user can zoom to zoomlevels between integers. It does not look okay for some raster layers.

Good (same as source raster map):

Little bit zoomed by OpenLayers (looks blurred):

How can I turn off this zoom option?


texts: Openlayers 2 - zoom to multiple coordinates

texts:
I have multiple asset icons on a map

I want to be able to zoom to show all icons on the map
I have the Lat and Long of each asset

Currently, i have it working with 1 asset as I can use this code to zoom to the asset
map.setCenter(new OpenLayers.LonLat(obj.d[k].vLongitude, obj.d[k].vLatitude).transform(WGS84, map.getProjectionObject()), map.getZoom());

Do I need to pass in an array of co-ordinates?


texts: Grouping counties into provinces, st_exteriorring(MultiPolygon)

texts:
I loaded into PostGIS (to be used in QGIS and web) a shapefile containing areas of a country. These areas are grouped in zones. I'm trying to create the shapes of these zones (and later an additional level).

The naive approach was:

create table eric_z as
  select
    z.id,
    z.name,
    st_multi(st_union(a.geom)) as geom,
    count(*)
  from geo_a a
    inner join geo_z z on a."Z_id" = z.id
  group by z.id, z.name;

This works but for some areas where the borders are not exactly matching, I'm getting inner points like this:

I then thought about using st_exteriorring. Since this produces outer LineString, I used st_makepolygon to get full shapes. The issue would not be fully visible in the screenshot below with just lines.

create table eric_z2 as
  select
    z.id,
    z.name,
    st_makepolygon(ST_ExteriorRing(st_union(a.geom))) as geom,
    count(*)
  from geo_a a
    inner join geo_z z on a."ZS_id" = z.id
  group by z.id, z.name

This seems to work at first but as predicted in comments by @ThingumaBob, it can be tricky because st_exteriorring can have issues with multipolygons. The result is that some areas suddenly become null:

@Thingumabob suggested to use st_dump to get the individual lines, exteriorring them and then turn it back to polygons. However, I can't figure out how to do that.

Something that does work but that I find awfully dirty is to enlarge the polygons by a short distance and then shrink them again. This actually does what I want but I'd feel dirty committing that code into posterity... :(

st_buffer(st_buffer(st_multi(st_union(a.geom)), 0.00001), -0.00001)

Any better suggestion ?


texts: Showing only selected features from a a layer in the QgsComposerLegend?

texts:
I have developed a custom qgis plugin. I have a button to export the map canvas in the form of image along with the legend. I came up with the following code for the legend:

    legend = QgsComposerLegend(c)
    legend.model().setLayerSet(mapRenderer.layerSet())
    c.addItem(legend)

The problem is that it adds all the features in the map legend whereas I only want to include the selected features in the map legend. I also tried:

    legend.setLegendFilterByMapEnabled(True)

But this does not seem to work. Any ideas on how to solve this issue?


texts: What might cause gdalwarp to lose CRS information?

texts:
I asked a question on StackOverflow focusing on the Docker aspect of my issue, here I would like to focus on the GDAL matter.

I have a Docker container in which gdalwarp -of GTiff in.JP2 out.tif is executed via Python's subprocess.run().

When I run the container in interactive mode, with docker run -it <image_name> bash and then pyhon script.py, I get the expected output. However, if I run the script with docker run <image_name> python script.py, the CRS information is missing from the output.

Below are gdalinfo outputs in both cases, see the missing parts in the second one.

Expected:

Driver: GTiff/GeoTIFF
Files: file.tif
       file.XML
Size is 19771, 19771
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (16.295439814814817,48.338634259259265)
Pixel Size = (0.000004629629630,-0.000004629629630)
Metadata:
  AREA_OR_POINT=Area
  METADATATYPE=DIMAP
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  16.2954398,  48.3386343) ( 16d17'43.58"E, 48d20'19.08"N)
Lower Left  (  16.2954398,  48.2471019) ( 16d17'43.58"E, 48d14'49.57"N)
Upper Right (  16.3869722,  48.3386343) ( 16d23'13.10"E, 48d20'19.08"N)
Lower Right (  16.3869722,  48.2471019) ( 16d23'13.10"E, 48d14'49.57"N)
Center      (  16.3412060,  48.2928681) ( 16d20'28.34"E, 48d17'34.32"N)
Band 1 Block=19771x1 Type=UInt16, ColorInterp=Red
Band 2 Block=19771x1 Type=UInt16, ColorInterp=Green
Band 3 Block=19771x1 Type=UInt16, ColorInterp=Blue
Band 4 Block=19771x1 Type=UInt16, ColorInterp=Undefined

Faulty:

Driver: GTiff/GeoTIFF
Files: file.tif
       file.XML
Size is 19771, 19771
Origin = (48.247101851844533,16.386972222229545)
Pixel Size = (0.000004629629630,-0.000004629629630)
Metadata:
  METADATATYPE=DIMAP
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  48.2471019,  16.3869722)
Lower Left  (  48.2471019,  16.2954398)
Upper Right (  48.3386343,  16.3869722)
Lower Right (  48.3386343,  16.2954398)
Center      (  48.2928681,  16.3412060)
Band 1 Block=19771x1 Type=UInt16, ColorInterp=Red
Band 2 Block=19771x1 Type=UInt16, ColorInterp=Green
Band 3 Block=19771x1 Type=UInt16, ColorInterp=Blue
Band 4 Block=19771x1 Type=UInt16, ColorInterp=Undefined

I'm curious as to what might cause this loss of information.


texts: Labeling selection of OSM layer in QGIS depending on certain tags?

texts:
I just imported original OSM data to export a pdf-map by assigning graphic attributes to the features I need.

The label tool in QGIS 1.8 allows to select certain features, like I can do with the graphic style. But I do not know and can not find how to label a feature like ref (street-nr) of a certain tag like highway = motorway (especially if all other line elements are labeled by name).

tags LIKE '%"highway"="motorway"%'

and then? ref is a feature with an attribute in the column tags.

Another Example: Can I achieve several labeling sizes for different tags like city, town or suburb?

The old labeling menu allows only one label attribution for each layer. Can I select labeling by any attribute?

Thanks for any help and tips concerning the syntax in QGIS attribution and labeling. (The answer: "How to label features in QGIS using SQL statement?" did not help at least)


texts: Get all points from multiline in shapefile using C# and GDAL

texts:
I'm trying to read a shapefile and draw it using DirectX and GDAL library. So first I will read the shapefile and get a list of all the points in the shapefile. It works with linestring using GetPoint().

var polylines = new ConcurrentDictionary<string, List<Point3D>>();
Feature feat;
while ((feat = layer.GetNextFeature()) != null) {
    var geom = feat.GetGeometryRef();
    var points = new List<Point3D>();
    for (var i = 0; i < geom.GetPointCount(); i++)
        {
            double[] geopoint = {0, 0, 0};
            geom.GetPoint(i, geopoint);
            points.Add(new Point3D(geopoint[0], geopoint[1], geopoint[2]));
        }

    polylines.TryAdd(feat.GetFID().ToString(), points);
    feat.Dispose();
}

But when the geometry is MultiLineString then GetPoint(), it always returns 0. Does anyone know if there is a way to unpack MultilLineString to LineString or read MultiLineString using GDAL?


texts: Simplify OS MasterMap Data without Losing Key Details

texts:
I am building a flood model (Infoworks ICM), using OS MasterMap Topography data supplied in shapefile format to define land use such as roads, buildings, fields etc across the river's catchment. The model uses a triangular flexible mesh, the vertices of which are partially influenced by the vertices of the various layers used in the model build process. Many of the "road" polygons within the OS MasterMap data contain a significant number of vertices around bends for example. Buildings often contain minor details such as porches or steps. These details mean that very small triangles are generated due to the high concentration of vertices within a small area of the OS MasterMap data. Smaller triangles often lead to significantly longer computational times or model instability/innaccuracy. For the purposes of the study, the buildings could be simplified heavily to rectangles or similar to remove the minor details.
I've tried the built in simplify tools in QGIS to reduce the vertex count within the OS MasterMap data, but with mixed success. Whilst I am able to alter the simplification method and tolerance to reduce the number of vertices, the polygons within the output often overlap each other which would in turn produce further small triangles within the model mesh amongst other issues.
Is there a way to simplify the OS MasterMap data by removing unnecessary vertices around curves, trimming minor details from buildings etc, but keeping the essence of the original data in tact and resulting in no overlaps or gaps? Ideally using QGIS 3.x only, but I have access to MapInfo 10.5, ArcMap 10.4 and Global Mapper 18 if necessary.
I've attached a few images showing how a the bend of a road has been simplified to reduce the number of vertices from the red points to green, however the algorithm seems to retain different vertices within each polygon, causing the overlaps/gaps.
Image 1: Before simplification - vertices shown as red points
Image 2: After simplification - vertices shown as green points with overlaps/gaps shown due to mismatched retained vertices
Image 3: Mesh layout showing smaller elements within buildings (green polygons). Issues highlighted in yellow. Ideal, simple polygons circled in black.


texts: Moving SpatiaLite Table into PostgreSQL using ogr2ogr?

texts: How can I move a Spatialite table from a Spatialite db to a PostgreSQL db using ogr2ogr?


texts: Converting Road shapefile with line geometries to another shapefile with polygon geometries by adding buffer to line on both sides using QGIS

texts:
I have Roads shapefile which are denoted by lines.

I want to create another shapefile of roads as polygons, by adding 10-meter buffer on line.

How do I do this using QGIS?


texts: Converting KML into a vector graphic... and back

texts:
Is there a simple means to convert, say, an SVG or other vector image format file to KML and vise versa?

There are incredibly powerful vector editing tools which attempt to replicate the way reality works, or at least the way it presents itself to the artist's eye. GIS manipulation tools are much more focused on the data itself, as opposed to the "whys" of the data.

There are cases where I would like to import a KML file into GIMP or Fireworks, re-touch it, and convert it back. Anyone?


texts: Comparing overlap with different SRIDs in MSSQL

texts: I am sourcing our GIS data (ESRI ArcGIS) into a data warehouse for reporting. I am trying to join our sites in Table_A (with a point geometry column) to our markets in Table_B (with a polygon geometry column) by using the STContains() function. However, I found that these two columns are using different SRIDs (4326 and 4269 respectively) so the function will not work. I've done a lot of reading, and found that SQL Server has no transform function. So I used the .STAsText() function to store the text and STGeomFromText function to recreate the geometry, but this time matching the SRIDs. Now the STContains() function returns 1 that they do overlap as expected. My concern is I don't fully understand the impact of changing this SRID. By doing so, did I shift the point/polygon from where it should be on the map? Is what I'm doing a valid solution?


texts: How to detect Query Layer for toolvalidation using arcpy?

texts:
For an arcpy toolvalidation script, I have been trying to figure out a good way to detect ArcGIS layer files (*.lyr) based on "Query Layers", and to distinguish them from "normal" feature layers pointing to a File or Enterprise Geodatabase. I haven been searching fruitlessly through the ArcGIS for Desktop Help and especially the 'describe' object properties, and Google to see if there was any clear and easy way to do this, but I haven't found anything yet.

I am likely to have overlooked something obvious, but does anyone have ideas of how to do this?

The ultimate goal is to be able to restrict layer files selected by a user for a certain parameter of the script to those that are not based on a Query Layer, so to exclude them.


texts: Converting ArcGIS style (symbology) to font for MapInfo/OpenGeo?

texts:
This is related to Retain ESRI symbology in Mapinfo TAB

As a workaround I am wondering if I can save all my ESRI symbology from the lyr files created using the script below to an ESRI style sheet and then convert that to a font set which can be opened in mapinfo/opengeo?

The image above shows some of the [complex] symbology that we have in our ESRI mxd's that need to go into Mapinfo and/or Opengeo to be available via WMS.

Tests show that KML created from the lyr files handle the symbology better than tab but not every time.

Code so far...

import arcpy, string, datetime, shutil, os
import arcpy.mapping as MAP

#Read input parameters from script tool
MXDList = string.split(arcpy.GetParameterAsText(0), ";")
ProjectPath = arcpy.GetParameterAsText(1)

#MXDList=r'P:\2012\183_TownPlanning_Symbology\Working\test.mxd'
#ProjectPath=r'P:\2012\183_TownPlanning_Symbology\Working'

count=0
#Loop through each MXD and print
#for MXDPath in MXDList:
#MXDFile=r'P:\2012\183_TownPlanning_Symbology\Working\Strategic_Plan_Ravenshoe_8.mxd'
count=count+1
print count
##try:
mxd=arcpy.mapping.MapDocument(MXDFile)
outPath = ProjectPath+'\\'+mxd.filePath[mxd.filePath.rfind('\\')+1:mxd.filePath.rfind('.')]
os.mkdir(outPath)
msd = outPath+'.msd'
dfList = arcpy.mapping.ListDataFrames(mxd)
for df in dfList:
#arcpy.mapping.ConvertToMSD(mxd, msd, df, "NORMAL", "NORMAL")
#arcpy.AddMessage(str(count)+'...' +str(outPath)+'...'+'\n')
    lyrList=arcpy.mapping.ListLayers(mxd, "", df)
    outPath = ProjectPath+'\\'+mxd.filePath[mxd.filePath.rfind('\\')+1:mxd.filePath.rfind('.')]+'\\'+str(df.name)
    os.mkdir(outPath)
    for lyrFile in lyrList:
        if lyrFile.isFeatureLayer == True:
            if lyrFile.visible == True:
                arcpy.AddMessage(str(lyrFile)+' is visible')
                print str(lyrFile)+' is visible...exporting dataset'
                outFileN=str(arcpy.ValidateTableName(lyrFile.longName[lyrFile.longName.rfind('\\')+1:]))
                try:
                    arcpy.FeatureClassToFeatureClass_conversion(lyrFile.dataSource, outPath, outFileN)
                    arcpy.SaveToLayerFile_management(lyrFile,outPath+'\\'+outFileN+'.lyr', "ABSOLUTE")
                    arcpy.LayerToKML_conversion(outPath+'\\'+outFileN+'.lyr', outPath+'\\'+outFileN+'.kmz')
                except:
                    print 'Exists...continuing'

        arcpy.AddMessage(str(lyrFile)+' is visible')
        print str(lyrFile)+' is not visible'
##except:
##    errorM=arcpy.GetMessages()
##    arcpy.AddMessage(str(count)+str(errorM)+'\n continuing...\n')
##    print errorM
##    #MAP.PrintMap(MXD, printer)
#Remove variable reference to file
del mxd, msd

BTW// I tried the msd route (change to .zip and open the xml) but it is very convoluted.


texts: Is there an easy generic way to cluster Markers with CartoDB?


texts:
in the new version of qgis 2.8, in the c++ API and in the class QgsSymbolV2, a lot of function has been deleted, by exemple setfillcolor, setcustomtexture, setnamedpoint, setpointsize

i would like to know if these function are useless in the 2.8 version ? or they had been moved to another classes ?

please help


texts: Avoiding arcgisscripting RuntimeError: Not signed into Portal from ArcPy with ArcGIS Pro?

texts:
Sometimes when I run a Python script from IDLE that uses the ArcPy installed with ArcGIS Pro 1.4.1 the script will abort with RuntimeError: Not signed into Portal:

RESTART: C:\temp\test.py
Traceback (most recent call last):
  File "C:\temp\test.py", line 9, in <module>
    import arcpy
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\__init__.py", line 66, in <module>
    from arcpy.geoprocessing import gp
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\__init__.py", line 14, in <module>
    from ._base import *
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 14, in <module>
    import arcgisscripting
RuntimeError: Not signed into Portal.

If I do nothing other than run the script again, in exactly the same way, using the same Python shell and script windows, then it will usually just work.  Very occasionally I may have to try running it a third time.  At a very rough estimate it works 90% first time, 9% second time and 1% third time.

I have seen the same error on two machines running two versions and two different license sources:

ArcGIS Pro 1.4.1 sourcing its license from ArcGIS for Portal
ArcGIS Pro 1.3.1 sourcing its license from ArcGIS Online

In both configurations I have not tested using Authorize ArcGIS Pro to work offline so all tests have been run while working online.

Does anyone know a likely cause for this error, and any way to try and avoid it occurring?

The same error seems to have been seen deep in a GeoNet forum thread and suggests that it is related to logging out of ArcGIS Online while running the script.  However, I am running my script while ArcGIS Pro is open (but not being used), and I am not logged in via a browser to ArcGIS.com so I think the advice offered there does not apply.


texts: ArcPy Python Script To Update Shapefile Records From CSVs

texts:
I am fairly new to python scripting. Most of the basics I have down however still learning.

I am trying to create a python script to update a shapefile with information from two separate CSV files. There is a common attribute that will link the CSV records to the Shapefile. I want to update the records of the shapefile with information from provided in the CSV files.

I found a script from a website that appears to do what I am looking for (link below), however when I run the script it it errors on line 11: for row2 in cursor2 and states: "Runtime Error; Invalid SQL Statement". I have tried several different attempts to negate the issue to no avail. If anyone could point out what is wrong with the script or if you have a more efficient way to create a script to accomplish my goal I am all ears.

import arcpy
from arcpy import env
env.workspace = "C:/Data"

print 'Processing...'

fc = "parcels.shp"

cursor = arcpy.da.SearchCursor(fc, ["APN"])
for row in cursor:
  SQL_stat= "APN = "+ str(row[0])
  fc2 = "Numerical_CSV.csv"

cursor2 = arcpy.da.SearchCursor(fc, ["ASSESSMENT"], SQL_stat)
for row2 in cursor2:
  UpdatedValue = row2[0]
  cursor3 = arcpy.da.UpdateCursor(fc2, ["TAXYEAR"],SQL_stat)
    for row3 in cursor3:
        row3[0] = UpdatedValue
        cursor3.updateRow(row3)
del row
del cursor
del row2
del cursor2
del row3
del cursor3

print "Done"


texts:
Is there any of you ever work to develop a web-mapping project for a plantation company? Or is there any such project published? If there is then :

Can you point where I can see an example of the published project?
What was the main goal of the application? is it only for viewing purpose or is it combined with reporting system?
What platform do you use? Free or Paid?

Best Regards,


texts: Where is official documentation of gdal/ogr?

texts:
Where can I find the official docs for gdal/ogr python API?

I am looking for something comprehensive which shows all the classes with descriptions of their attributes, methods and input/output types. So far I have been using https://www.gdal.org/annotated.html but this appears to have been generated for an older version of gdal/ogr. For instance I am using gdal version 2.2.2 and it appears that in the docs the OGRDriver class is not shown, only the GDALdriver class is shown. At the bottom of the link the following is stated:

  Created for gdal by doxygen 1.8.8

I'm not sure if 1.8.8 refers to gdal version or doxygen version, in any case I would like to find the proper documentation for the version of gdal I am using?


texts: Changing feature attribute by selecting feature in vector layer using QGIS

texts:
I wanted to change a feature attribute by selecting a feature in a vector layer and opening the feature form. The following QGIS online course does only show the procedure for an older version of QGIS: https://docs.qgis.org/2.2/en/docs/training_manual/create_vector_data/forms.html.
Is there a way to accomplish the same in QGIS 2.18.9?


texts: difference in elevation between closest pixel of another layer or 3d distance

texts:
I want to find the closest point of water to each forest pixel, and the difference in elevation, to see how deep should a channel be dug in order to get water to the mangrove.

My aim is to obtain a raster with the values of the subtraction between a layer and its closest point in another elevation layer

I want to keep it raster (finding closest distance between polygons may take hours, while it needs just a couple of minutes in raster). It can be in GRASS, R, QGIS or python. I've been trying with r.distance and library(raster) without success.
It may be a 3d distance I guess.


texts: Is there a minimum number of required samples/ data points for IDW Interpolation in ArcGIS for Desktop?

texts: I am interpolating air pollution values from 6 data locations using IDW in ArcGIS. I know that more the number of points, more reliable the results are. But, I am not sure if a minimal number of data points is necessary.


texts: Adding dynamic text to layout from attribute table with Data Driven Pages

texts:
I was trying to figure out how to use dynamic text from an attribute table in a collection of Data Driven Pages in ArcGIS 10.

I've looked in the help files, but can't find anything on dynamic text from attribute tables.


texts:
I have a data set of ~3500 points located across the Pacific Ocean, and am looking to measure the distance from each point to it's nearest coastline to identify the "host" country. I had originally used Plate Caree, but now believe that it doesn't maintain distances accurately?
I am hoping not to need a country by country approach, but could someone confirm Plate Caree does / does not maintain accurate (reasonably accurate) distance, or suggest another?
Also, I am fully a user of ArcMap, not QGIS, etc. so 'how to' answers outside ESRI aren't helpful in this instance.


texts:
A question regarding mapnik:

Is it possible to configure a custom max. zoom level? I am trying to find a way to get mapnik to render tiles for zoom levels 20+ (till 28 would be nice).

PS: The solution according to this link does not seem to exist anymore.


texts: How to create a irregular buffer around a point?

texts:
I need to create an buffer around a point which is 3.58km x 3.0km. i.e I need two sides of the buffer to be 3.58km long and the other two sides to be 3.0km --> making a rectangular shape. I am finding this task difficult as the buffer area is not an equal area. does anyone know how to do this?

I have tried to, however There are a few issues I face with this.

Firstly, my coordinate reference system is in degrees not meters. Would it be accurate if I was to convert the CRS from WGS84 to once which uses meters instead of degrees. then type in the specific dimensions I want --> save the shapefile and then convert it to WGS84.

Secondly The buffer tool does not have a viable option to input two data points. Would I need to save the data points (i.e 3.58km x 3.0km) In Excel and save it as a text file then upload it and enter that data set into the buffer tool?


texts: Loop points in polygon based on subset of two columns in R

texts:
I have a spatial point data frame with two columns with attributes, one for months, another for species (and off-course the lat, long for those points).

long=c(-1.747429,-1.519554,0.430455,-2.049941,-1.288655,1.311327,-2.217094,0.718583,-1.646491)

lat= c(54.407632,53.369744,51.519982,52.347591,54.530768,51.12631,53.430349,51.552895,55.028696)

species = c('d','c','d','d','d','c','c','c','c')

month = c(1,2,3,2,3,4,3,2,3)

I also have a spatial polygon of my country projected in the same CRS as my points.

I want to create a loop that will give me the count of points per month and species and that output to be in separate data frames.

But so far I only know how to do count points in polygon overall

Df2<-over(Df, poly)

df2Freq<-count(Df2, "NAME_1")

What I want is a loop that will do by month and species and produce multiple data frames with outputs


texts: Converting ArcMap Cartographic Representations to ArcGIS Pro and exporting stylx

texts:
Here I have a Feature Class named ContactsAndFaults, which I have imported into ArcGIS Pro from a File Geodatabase (you can find this same geodatabase here if you want to reproduce this analysis). The geodatabase was created in ArcMap and Cartographic Representations (CartoReps) were used for the styling. CartoReps are not supported in ArcGIS Pro. However, you can draw your symbology from CartoReps by using Appearances -> Symbology ->  Symbolize Your Layer Using Representations.

This works and correctly renders the lines (you can see the difference by comparing the pictured screen shots as some lines have now become dashed etc.) It does this by converting the CartoReps to Unique Values.

Now, would like to BATCH export this styling information as an ArcGIS Pro .stylx file, but how do we do this?


texts: Label dynamic text but also shorten it to only 50 characters

texts:
I am using dynamic text and labeling this:

<dyn type="page" property="number"/>

but I only want it to label the first 50 characters and than if it has to shorten the name to than put et al. at the end.

Thoughts?


texts: How to interpret ArcGIS kernel density units when a population parameter is specified?

texts:
Possible Duplicate:
  How to interpret GRASS v.kernel results?

The attached image shows trees in a study area based on their actual height in the field.  The trees are overlaying a raster that was produced in ArcGIS using the Kernel Density tool.  A population field was used for the kernel density estimate which included the tree heights in meters.  The following descriptive statistics are for the displayed trees in the image:

n = 117
Min height = 0.2 m
Max height = 8.68 m
Mean height = 0.76 m

How should I interpret the results of the kernel density estimate as displayed in the legend?
You can see the largest KDE value is approximately half the average tree height.  What am I missing here?


texts: Re-projection of a shape (with a dbf) and converting the original X Y data in ArcMap

texts: In my database, I have a dataset (shp, dbf etc) which has X Y coordinate data that was taken in NAD27 from a GPS (a long time ago). When importing it in ArcMap desktop I can easily project the data in NAD83 with the accurate transformation (NtV2) taken for the region. And thus, it shows me my data correctly projected in NAD83 by creating a .prj file if I understand it well. HOWEVER, it does not change the original NAD 27 X Y coord data in the dbf. Is there any way/tool to extract the NAD83 X Y data resulting from the reprojection?


texts: Change Defination query of a layer in ArcGISDynamicMapServiceLayer

texts: Is there any solid example to change the definition query at run time of the ArcGISDynamicMapServiceLayer using javascript API 3.5 or higher version and ArcGIS server 10.1. I found API documentation setdynamiclayerinfos but was not able to find example to find how it exactly works.


texts: Using arcpy search cursor to relate a .csv and a .shp?

texts:
I'm trying to write a script that relates data from a .csv file to a utility pole shape file. Basically, we want the script to iterate through the pole # fields in each table and determine which pole numbers are missing from the .shp. Right now, the script iterates through every row of the field and appears to not be reading the pole # value inputted by the user, since every print out reads the pole # and NO (ex: 111A0111, NO) no matter if the pole # is matching or not. We want the program to read "YES" if the pole number is in both tables, and "NO" if it is only in the .csv file and not the .shp.

import arcpy
import csv

#define other variable
fc = "P:\Scripts\data_check\prd_poles.shp"
field = "GISID"

#Create a list to be populated by the search shapefile search cursor
gis_poles_assets =[]

cursor = arcpy.SearchCursor(fc)

#Populate the list
for row in cursor:
  gis_poles_assets.append(row.getValue(field))
#print gis_poles_assets

#define the csv to be read by the reader and make that a variable
f = open(r'P:\Scripts\data_check\discoverer_poles_snip.csv')
csv_f = csv.reader(f)

#search the shapefile based on fields in a .csv
for row in csv_f:
  if row [1] in gis_poles_assets:
   pass
  elif "111A0111" in gis_poles_assets:
    print row [1], "YES"

  else:
    print row [1], "NO"


texts: Conversion from CGIAR GEOTIFFs to SRTM3 like HGTs

texts:
Im writing a C program for querying height information from SRTM data. Since reading directly GeoTIFF format is tough job, I decided to load HGT format, which serves great for the purpose.

Since SRTM3 data contains void areas, I found the CGIAR preprocesed data which looks much better. But they only offer geotiff download for 5×5° altogether saved in continous numbered filenames (thus not compatible with degree number naming).

I would like to convert these TIF files (ie srtm_39_03.tif) to HGT and name them properly. Or alternatively please suggest another download site.

So far I realised, that I can convert 1° TIF to HGT easily, but 5° fails:

$ gdal_translate -of SRTMHGT srtm_39_02.tif srtm_39_02.hgt
Input file size is 6001, 6001
0Warning 1: The corner coordinates of the source are not properly aligned on plain latitude/longitude boundaries.
ERROR 1: Image dimensions should be 1201x1201 or 3601x3601.

As of the splitting, it could be probably achieved with gdal_translate -srcwin 0 0 1201 1201 srtm_39_02.tif out.tif, but Im really confused how to split 6001 in 5 parts each 1201 pixels :-)

Thanks,
zbycz

Update: according to SRTM manual, HGT files overlap with the 1 upmost/rightmost line, so the splitting is clear now. But still would be happy for consultation.

Update2: finally I used SRTM hgt data and wrote my own ANSI-C SRTM parser.


texts: Seting directory automatically in ArcGIS GUI Batch Commands?

texts:
When I do a batch command in ArcGIS e.g. Projection transformation the batch dialog has a field for the "new file name".

However, the default Directory is always the same folder as the original files.

If I want to change the Directory but keep the old names, then there is a strange syntax which I can use, something like %n%, but I would still have to copy and paste the complete path to the new folder.

Is there an option to simply pick a new folder and have the rest complete automatically?


texts: Conditional to update field not working properly - Arcpy

texts:
I am attempting to update a text field with a shortened version of what is currently there.  Currently, there are values from 1 to 4 with a description of each.  I would like to return only the value.  I tried a few different snippets of code to achieve this.  However, when I run the code, only some of the values change as I would like.  My most current code is:

try:
    with arcpy.da.UpdateCursor(protectTable, "gap_stat") as cursor:
        for row in cursor:
            for n in row[0].split():
                if n[0] == "1":
                    row[0] = "1"
                elif n[0] == "2":
                    row[0] = "2"
                elif n[0] == "3":
                    row[0] = "3"
                elif n[0] == "4":
                    row[0] = "4"
                else:
                    pass
            cursor.updateRow(row)
except AttributeError:
    pass

Here is an example of the result where some of the values change and some do not:


texts: Publishing tools fails to start: Error: service failed to start System.PublishingTools.GPServer:  Service's containing process has crashed

texts:
My question is about PublishingTools failing to start.

The publishing tools service fails with the error message:

  "service failed to start System.PublishingTools.GPServer:
  Service's containing process has crashed"

I have uninstaled the Spectre and Meltdown security patches as well but the problem still exists.

The server logs give an error code 9003 with message

  Cannot obtain a free instance.; nested exception is:
  com.esri.arcgis.discovery.ejb.ArcGISServerEJBException:
  java.lang.Exception: Could not initialize service
  'System/PublishingTools.GPServer'.

How do I fix this issue?

I have ArcGIS Server 10.1 on Windows Server 2008 R2.


texts: Seeking QGIS plugin to create Isopleth map

texts:
I'm looking for a plugin for QGIS that would allow me to create an effective isopleth map.

Does one exist?

I haven't been able to find any, only independent software that would allow this type of mapping.


texts: Copying new data into shared Geodatabase in ArcGIS for Server using ArcObjects?

texts:
How do I need to in order to publish a geodatabase into a server or just to copy new data into the server using ArcGis Server 10.1, ModelBuilder and .NET 2010.

Note : I have already a geodatabase published on the server but each time I have a new added data on the geodatabase they must be copied into the server. to be able to consult them on the map.

Do I need perhaps to a script ArcObjects, ModelBuilder or an addtional tool?


texts: Rasterfile clipping error: Cannot compute bounding box of cutline

texts:
I have a raster file and want to clip it, a vector file is supposed to be the mask layer. I use the most obvious method (Raster->Extraction->Clipper).
No matter what I do, I get this error (see picture below):

  Cannot compute bounding box of cutline

I tried it with changing the raster format, changing the projection etc.
Of course I made research in older posts/Google too. I have no Idea how to solve this problem.  If I clip by extend it does work.

Weird thing is, I remember doing this three months ago, same method and I did not get this error.

I also would (instead of a solution for this error) welcome an alternative way to perform this action. I just need to clip a rasterfile as described. Can someone help?


texts: Down sample with imagemagick then import worldfile?

texts:
I'm trying to get a better understanding of the limitations of world files and how I can use them in a current project.

If I produce an ESRI world file (tfw) with listgeo:

listgeo -tfw clip.tif

Then downsample the image with ImageMagick and use in geotifcp:

geotifcp -e clip.tfw output.tif final.tif

Will the georeferencing be accurate on the downsampled image? If not, how do I address this situation?


texts: Add Point to Layer WITHOUT Typing in Coordinates

texts:
I have a point layer and a line layer in my QGIS project (QGIS 3.0.3). I want to add a point to the point layer at the intersection of two lines in my line layer. I toggle editing for the point layer, select "add point feature", locate the cursor on the intersection and left click.

I get a dialog box for the point with three fields: EAST, NORTH, and ID. I type in the ID, without typing in the EAST and NORTH coordinates because they are 12 significant figures and therefor are HARD to type in. Furthermore, the numbers in the "coordinate" display at the bottom of the screen move around constantly even though I have already fixed the point by clicking on the position I want it at. SO, I put in the ID (simple) and click OK. The point, with the ID, displays right where I want it to be. Great!

Having created my point, I then SAVE LAYER EDITS. Immediately my point DISAPPEARS! When I open the attribute table for the layer, the point is listed with NULL / NULL for X and Y. And there is no point on the layer display, even though there was one before I clicked SAVE LAYER EDITS.

How do I get QGIS "Add Point Feature" to AUTOMATICALLY post the X and Y values automatically into the attribute table since I have already indicated the coordinates by clicking on the position where the point is to be placed? Or automatically into the dialog box, instead of leaving them NULL?

The answer to the question "Update point coordinates in QGIS after manually moving locations in edit mode" does not answer this question. At the time that a point is CREATED (not moved) a dialog box is presented that requests that the X and Y (EAST and NORTH) coordinates be entered. Since QGIS knows the X and Y (EAST and NORTH) coordinates at that point, I want to find an easy way to put those values into the dialog. That is not what the other question is about, nor does it have an answer that would work in this context.


texts: How do I mirror a vector layer along a line in QGIS?

texts:
ArcGIS (Map & Pro) has a nifty utility that allows you to mirror a vector layer along a line segment you give it.  Before I spend time trying to write something from scratch, I wanted to see if there was a workflow where I could accomplish the same in QGIS.  I need actual features here, not just displayed graphics.
Here's an example of some dots and the line along which I'd like to mirror them:

I've played around with the affine and v.transform (grass) tools in QGIS 3.14.  Neither appeared to allow me to issue the tool a line segment to use as the reference. I installed Vector Bender as well, but it doesn't seem to work for this use case (though the documentation is too slender to tell).


texts: What is my coordinate system preserving?

texts: I'm trying to get a better grasp of coordinate systems, datums, projections etc. I understand that certain systems preserve different things; distance, shape, angles, area. I'm wondering how you can tell what a certain system is preserving. What portion of the CRS name indicates this? Specifically, I'm wondering about a UTM system in Alberta, Canada. NAD83(CSRS)/Alberta 10-TM (Forest)


texts: How to prevent transparent borders (collars) of a tif wipe out overlapping regions of adjecent tif while using gdalbuildvrt

texts:
I have a set of Geotifs that have got transparent borders (collars) of irregular shapes & dimensions around them. First I created a VRT over this set of Geotifs through gdalbuildvrt. Then I generated tiles through gdal2tiles.py using this VRT.

I noticed that the transparent borders (collars) of some Geotifs wiped out the overlapping regions of the adjacent Geotifs creating holes in the output.

Is there any setting that I should use while creating VRT through gdalbuildvrt or while generating tiles through gdal2tiles.py and the above created VRT, so that I can get rid of the holes in the output caused by the transparent borders (collars) of some Geotifs wiping out the overlapping regions of adjacent Geotifs.


texts: I have openstreetmap tiles for offline (png images), how can I display these images (tiles) on my application with zooming and rotation functionality as on openstreetmap.org using c++?


texts: Making multiple fields NULL using Python in ArcGIS?

texts:
I have a repetitive task here at work. It requires that I make multiple fields NULL as I update the geodatabase. Is there a script I can write that I can input the fields and the feature class. Then make all those fields NULL. This script must not make everything NULL in the feature class, only the fields designated.
    import arcpy

fc = "N:\Updates\CarsonCityNV\CarsonCityNV.gdb\J_POC_S_Studies_Ln"

with arcpy.da.UpdateCursor(fc, ("C1_GAGE", "C2_DISCH", "C3_MODEL", "C4_FCSTR", "C5_CHANN", "C6_HSTR, "C7_SCOUR", "S1_REGEQ", "S2_REPLO",     "S3_IMPAR", "S4_HSTR", "S5_CHIMP", "S6_TOPO", "S7_VEGLU", "S9_HWMS", "S10_REGEQ", "CE_TOTAL", "SE_TOTAL", "A1_TOPO", "A2_HYDRO", "A3_IMPAR", "A4_TECH", "A5_FOAPASS")) as cursor:
 for row in cursor:
    row[0] = None
    row[1] = None
    row[2] = None
    row[3] = None
    row[4] = None
    row[5] = None
    row[6] = None
    row[7] = None
    row[8] = None
    row[9] = None
    row[10] = None
    row[11] = None
    row[12] = None
    row[13] = None
    row[14] = None
    row[15] = None
    row[16] = None
    row[17] = None
    row[18] = None
    row[19] = None
    row[20] = None
    row[21] = None
    row[22] = None
    cursor.updateRow(row)

I'm getting a parsing error on line 5.


texts: Iterating using search cursor and layer selection performance?

texts:
I'm trying to use the results of a search cursor to iteratively make selections on a feature class and then create new feature classes (so you get a new feature class or shp for each unique value in a field... kind of an explode by attribute tool).

As an example - if you had a feature class containing 14763 records and one of the fields had 18 unique values in it, after running the script you would have 18 new feature classes in your gdb.

So far I have this:

def unique_values(table, field):
with arcpy.da.SearchCursor(table, [field]) as cursor:
    dict = sorted({x[0] for x in cursor})
    print type(dict)
    for i in range(len(dict)):

        whereD = dict[i]
        fname = whereD[0:5] + "_" + whereD[-4:]
        arcpy.MakeFeatureLayer_management(table, fname)
        query = "relatedBirds_Species = '" + whereD + "'"
        arcpy.SelectLayerByAttribute_management(fname, "NEW_SELECTION", query)
        arcpy.CopyFeatures_management(fname,fname+ "_c")

This works but it's incredibly slow when it runs.  Does anyone know of a quicker way to get the same results?

Requested edits below

The Data

Geodatabase workspace
The feature class I'm running this with has 18 unique values in the field in question.
Total number of records in the original feature class is 14,673

Timings

As the code is above = 2' 35.04"

With the {where clause} added to MakeFeatureLayer (w/out copy feature either)= 1' 27.75" (thanks Stephen)

With Select_analysis = 2' 10.58"

I think most of the time was being taken to actually draw these in Arc as I was running this directly from the python window so I ran it again from cmd with arc closed.

Stephen wins at 57.93" (with copy feature put back in) running from cmd.


texts:
New-ish to GIS.

Using Arc 10.3 and experiencing an annoying problem. When attempting to resize a neatline (in order to use it as a locator box on my inset reference map) the lines disappears when reduced in size too quickly, or after getting below ~0.6x0.6". The blue vertices remain visible afterwards, but cannot be selected/moved/brought forward.

Any thoughts or work arounds?


texts: Setting PostGIS default sequential value in a QGIS "editable" view

texts:
Introduction

In a QGIS project I have a VIEW, that I want to became editable (I'm trying to  replicate the versioning system used DB Manager with some tweeks). Therefore I created a few  rules to make is so (see below). I'm able to delete and update the features without problems, but while inserting new features, I'm not being able to use PostGIS default value (sequence) for gid (that in this case is not a primary key) but has a NOT NULL constraint. So, if the user don't fill the gid value, QGIS is not be able to commit the changes, if it does... we never know where thing going to get wrong.

The question

Is there a way I can prepare my view or set QGIS form to read the default value?

What I got

-- The original table
CREATE TABLE "PGHP".infraestruturas
(
  gid serial,
  nome character varying(40),
  classe character varying(40),
  geom geometry(POINT, 3763),
  "id_hist" serial primary key,
  "time_start" timestamp,
  "time_end" timestamp,
  "user_update" character varying(40)
);

-- Creating the view
CREATE OR REPLACE VIEW "PGHP".infra_portoes AS
  SELECT "gid", "nome", "geom"::Geometry(POINT, 3763), "time_start", "user_update"
  FROM "PGHP"."infraestruturas_pontos"
  WHERE "classe" = 'portoes' AND "time_end" IS NULL;

-- Make the view editable
CREATE OR REPLACE RULE "_DELETE" AS ON DELETE TO "PGHP".infra_portoes DO INSTEAD
  DELETE FROM "PGHP"."infraestruturas_pontos" WHERE gid = old.gid;

CREATE OR REPLACE RULE "_INSERT" AS ON INSERT TO "PGHP".infra_portoes DO INSTEAD
  INSERT INTO "PGHP"."infraestruturas_pontos" ("gid","classe","nome","geom")
    VALUES (NEW."gid",'portoes',NEW."nome",NEW."geom");

CREATE OR REPLACE RULE "_UPDATE" AS ON UPDATE TO "PGHP".infra_portoes DO INSTEAD
  UPDATE "PGHP"."infraestruturas" SET
    "gid" = NEW."gid",
    "classe" = 'portoes',
    "nome" = NEW."nome",
    "geom" = NEW."geom"
  WHERE gid = NEW.gid;


texts: QGIS 3: How to iterate in a QgsPrintLayer with pyQGIS

texts:
I'm looking for a way to iterate through items from a composer in QGIS 3, with pyQgis.
First I create a QgsPrintLayout and add it to the layouts manager in my project.
    self.newLayout = QgsPrintLayout(projet)
    self.newLayout.initializeDefaults()
    manager.addLayout(self.newLayout)

Then I load a template to populate the QgsPrintLayout
    doc=QtXml.QDomDocument()
    doc.setContent(file1, False)
    self.newLayout.loadFromTemplate(doc,QgsReadWriteContext())

And I want a list of all the QgsPrintLayout items:
    labels = [item for item in self.newLayout.layoutItems()\
    if item.type() == QgsComposerItem.ComposerLabel]  ...

But Qgis answers:
    labels = [item for item in self.newLayout.layoutItems()\
    AttributeError: 'QgsPrintLayout' object has no attribute 'layoutItems'

But in doc for QgsLayout (from which QgsPrintLayout takes a part of its methods), I can read:

void  layoutItems (QList< T *> &itemList) const

Returns a list of layout items of a specific type. More...

And I really can't find another Public Member function to list all items from a QgsPrintLayer.
What am I missing here?


texts: In transferring a family cemetery layout to a digital format, I need a software program

texts: We have a large historic family cemetery. I've had surveying experience, drafting in the days when it was done by hand, and a lot of real estate experience with deed descriptions and layouts on maps. What I don't have is experience with this in the CAD digital world although I use the computer daily for other work. We want to transfer the information from our paper layout (some 2 acres with 300+ burial sites as well as additional information) to one that's digital to make it easier to make changes. There are many programs out there for commercial use - expensive - and many that are free or inexpensive. Can somone provide any advice on which of the latter to choose?


texts: Alignment and distribution tools for shapefiles?

texts:
Does anyone know of a native Esri tool, Python, or R package that would allow one to emulate arcmap's graphic element toolbar functionality, but for shapefiles rather than graphic elements?

For example, the image below shows a layer containing multiple rectangles of varying size arranged in a quasi-gridlike pattern.  Is there an existing tool that would allow me to nudge each rectangle such rectangles within some tolerance distance become perfectly aligned (using either an edge or the centroids) and distributed such that the spaces between them are uniform?  Note that I'm looking to shift the entire polygons rather than moving only some of the vertices (e.g. by using using ESRI's snap tool).

If I were to do it manually, I could envision creating a 2nd layer composed of one or more rectangles and snapping features in the 1st layer to the nearest feature in the 2nd layer.  However, this wouldn't necessarily handle the spacing between them properly.

The suggestion from @klewis would work well for a one-off type situation.  I'm more interested in an automated and reproducible/versatile strategy though.  I'm using these "bounding boxes" as an intermediate step in the creation of "callout tables" (image below) based on data collected at points from which the leader lines emanate.  I'm using the R arcgisbinding package to build the "tables" from the bounding boxes joined with the source data.  As you mention, conversion to graphic elements would eliminate any possibility of performing table attribute joins.

I can envision adding a grouping field to the attributes of my rectangular table outlines that will allow me to generate a bounding box from a union by group.  These group bounding coordinates could probably be used to do the alignment and distribution programmatically.  I just wanted to avoid recreating something that already exists if at all possible.


texts: Statistic for measuring (local) autocorrelation of ordinal data

texts:
I have a spatial data set with a bunch of variables in nine-level ordinal scale.

I want to calculate (local) spatial autocorrelation for each of these ordinal scale variables.

I've been using spdep-package for R. Join count -statistic seems to be the way to go for categorical data, but I'm unsure whether it works with multi-level and ordinal data.

What is the proper way to measure local autocorrelation or clustering of ordinal data?


texts: Can't extract chlorophyll a values from NetCDF file using R

texts:
I'm trying to get the chlorophyll a values from a NetCDF file using R software, but all I'm getting are missing values, NA. I would like to know if I'm doing something wrong or if the file really just has missing chlorophyll a values. I can get the Longitude and Latitude values with this method.

The file I'm using is from here https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/Mapped/Monthly/4km/chlor_a/ and I get missing values from any file I tried not only the one showed on the script.

require(rgdal)
require(maptools)
require(raster)
require(sp)
require(rorwr)
require(RNetCDF)

clorofila<- "C:\\Users\\User\\Desktop\\files\\A20172132017243.L3m_MO_CHL_chlor_a_4km.nc"

cla <- open.nc(clorofila)

print.nc(cla)
file.inq.nc(cla)

clor <- var.get.nc(cla,"chlor_a",start=c(1,1),count=c(8640,4320))
Long <- var.get.nc(cla,"lon")
Lat <- var.get.nc(cla, "lat")

With ncdf4 and raster I got the same results

require(ncdf4)

clorofila10<- "C:\\Users\\User\\Desktop\\files\\A20172132017243.L3m_MO_CHL_chlor_a_4km.nc"

nc <- nc_open(clorofila10)

val <- ncvar_get(nc, "chlor_a")
nc_close(nc)

raster

require(raster)
clorofila10<- "C:\\Users\\User\\Desktop\\files\\A20172132017243.L3m_MO_CHL_chlor_a_4km.nc"
clacla<-raster(clorofila10)
CHL1 <- raster(clorofila10, varname="chlor_a")
names(CHL1) <- 'chlor_a'

z <- getValues(CHL1)


texts: qgis_process run script from command line

texts:
qgis_process allows us to run scripts already installed in QGIS from the command line. But is there a way to run qgis scripts not yet installed?
e.g. if I have a file script.py
from qgis.core import QgsProcessing
# (...etc...)
class Foo(QgsProcessingAlgorithm):
    # defines my algorithm, which works fine in qgis

Can qgis_process or some other tool run this from the command line without me having to add script.py to the toolbox in qgis first?
The command-line help to qgis_process is encouraging in that it says

(...) run     runs an algorithm. The algorithm id or a path to a model
file and parameter values must be specified. Parameter values are
specified via the --PARAMETER=VALUE syntax

But if I provide a path to script.py (with or without the .py) it says not found.


texts: Group legend in ArcMap basing on field value

texts:
I'm using ArcMap 9.3.1 and have a point feature with an attribute table like this:

What would be the easiest way to generate a legend like this?

Generating different tables (feature classes) for each of the 3 states is no solution since the table is growing over the next month and it would be nice to use the same layout for other data as well.

I thought of creating 3 legends (one for each category) but wasn't able to display entries with a specific state.


texts: QGIS, displaying attribute values over geographies

texts: Is there a way to use graduated symbols that display a value off the attribute table? For example, I have a map of Vermont towns, and I would like to have an industry code floating above each town (I am using a color gradient to show unemployment rate, so I need another method to compare industry). I'm guessing I cannot, because I cannot reference to this anywhere or previous questions on the issue.


texts:
I am trying to concatenate two fields into one new field in GRASS, but I am unable to get the syntax to work. Ideally, I would like "FIELD1" and "FIELD2" concatenated together in "FIELD3".  All fields are VARCHARs.  I have tried:

v.db.update layername col=FIELD2 value=FIELD1+FIELD2

using many combinations of +, ||, &, and wrapping in single and double quotes all to no avail. Can someone provide some guidance on the required syntax?


texts: prevent ogr2ogr from generating new fid attribute field

texts:
I am exporting data from PostGIS to gml using python script with ogr2ogr. But the output generates unwanted new field  "fid" that contains the following string: "sql_statement.0" , "sql_statement.1", ... and so on. Is there a way to prevent this? I already have a field with unique feature identifiers.

cmd = ogr2ogr + " -f  \"GML\" " + filename + " PG:" + db_connect + " -sql " + query


texts:
I have a Site reference feature service data in polygons served up by ArcGIS Server.  I need to get the centroid coordinate of a polygon by the Id but can’t find a method out of the box.

I’m working on a GeoEvent Processor Service that will use this Site service to grab a centroid from a site Id.  So far I believe my options are:

Custom GP service, python – takes in polygon rings, returns point
Server Object Extension – same as above
Sql view, query layer served up as a feature service that uses STCentroid (The database uses a SQL GEOGRAPHY type) on the Site layer.  Would query the layer by ID to return the geometry.

The first two options keep me on track with a service orientated architecture however it certainly seems overkill to be writing an SOE for one function.  The final option could be a ‘goer’ but it means my logic is coupled to SQL Server.


texts:
I am getting the

  Unable to synchronize replica, sync replica internal

error in the collector app. We have interns out in the field and this just started happening after the second week, so I am confused. We can't just reload the map because then we could lose 50 data points that they have collected throughout our county.

Am I missing something simple?


texts: GeoFence: How to authenticate to WFS service?

texts:
I have Geoserver and Geofence properly installed and configured. I have rules in place for different users and their access to various layers. HOWEVER, when I send a WFS request for a layer, how can I add authentication info for that user?

Ex:
http://myserver:8080/geoserver/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=rapid%3Aecologically_significant_areas&maxFeatures=200&outputFormat=text%2Fjavascript&format_options=callback%3A%20getJson&callback=getJson&_=1463949789269

I tried putting the authorization in a header to no avail. On the server logs I see the following when I make the request:

Requesting access for RuleFilter[user:""+ role:ANY inst:name+:default-gs ip:"---.---.---.---"+ serv:"WFS"+ req:"GETFEATURE"+ ws:"rapid"+ layer:"ecologically_significant_areas"+]

How do I tell Geofence which user to use?


texts: ArcPy: (SearchCursor/export to tiff ) with no output

texts:
I'm trying to write a script that uses a search cursor to loop through a field in a  .shp file (Incident_Z) find every value that meets a specific condition (zip code=10004) and then at a scale of 1000, the layout at that scale to a pdf.

With advice from various folks I've tried two approaches:

Using a search cursor. However I've run into problems with it:
The problem is that there is no output, and interestingly, no errors. The script runs, and when I include a print "Done" statement at the end, prints correctly.

import arcpy

mxd = arcpy.mapping.MapDocument("C:\temp\Graffiti_DM1.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
expdir= "C:\temp"

cursor=arcpy.SearchCursor("C:\temp\itscomplicated.shp",'"Incident_Z" = ' + "'10004'", "", "", "")

for cur in cursor:
    cur.rotation = 0
    cur.scale = 1000
    outFile = expdir + "\" + cur.name + ".tif"
    cur.panToExtent(lyr.getSelectedExtent())
    arcpy.mapping.ExportToTIFF(mxd, outFile, df)

del cursor
Additionally, I've tried as per @PolyGeo's suggestion, the definition query method, below, which runs and produces an output, but not exactly what.I'm looking for.. The problem is that the output is a pdf (I'm actually preferential to pdf than tiff) of the current extent of the layout, and only the current layout. What I'm looking to do is to have a page for each instance of the where_clause, i.e. each record where '"Incident_Z" = ' + "'10033'", at a specified scale of 1000. Is it possible to modify the existing definition_query script to achieve that? As always, any advice is greatly apprecaited.


texts: intersection between two multipolygons yielding anomalous GeometryCollection object full of LineString's & Polygon's (trying to get intersect area)


texts: How to crop a river from a raster layer?

texts:
So I need to include the extent of a river into my map, but only have a raster file that shows the outline of the river with a blank (white) background:

Is there any way to crop this file so that my result is a raster which includes only the river without the white background so that I can put the layer above for example a ortho and still see the ortho with the river on top?

Alternatively, a raster to vector conversion should work as well, if possible?


texts: Would it be possible to construct a stand-alone grass.script python package?

texts: I would like to migrate dependency on ArcPy to grass.script without requiring the full installation of Grass. Would it be technically possible to have a portable directory containing grass.script and all the libraries it depends upon?


texts: QGIS - Creating a batch of drawings by automatically swapping layers

texts:
I have a project where I need to create several drawings showing model results in QGIS. I have 60 drawings to create where everything is the same apart from the model results themselves and also the drawing title.

I have anonymized these maps as an example, both show model results of flooding on a river and in this example the results are polygons and the drawing title changes. I need to do something similar, 60+ times, where only the results and drawing title change.

Is there an efficient way of doing this?


texts: Share ArcGIS Desktop and Pro 2.7 with FME Desktop

texts:
Actually, I use FME 2020.2, ArcGIS Desktop and ArcGIS Pro on the same computer.
FME uses the ArcGIS products to read/write the data in a Enterprise geodatabase.
Since I installed ArcGIS Pro, I cannot read/write SDE anymore
C:\Users\Stanusjl\AppData\Roaming\Esri\ArcGISPro\Favorites


texts: Creating random Image in Google Earth Engine

texts:
The final input of my work will be quintile ee.Image objects. To test my display I would like to create an ee.Image with the value in each pixel set randomly in [1, 5] (only integer).
I have already an ee.Image that could be used as a base:
var data = ee.Image('Oxford/MAP/accessibility_to_cities_2015_v1_0')

I just want to rewrite every pixel in it.
Is there a simple way to do that in Google Earth Engine ?


texts: AGOL & SQL database: Is ArcGIS for Server necessary?

texts:
A client is requesting we come up with a Python script to have an AGOL feature layer and a SQL Server database containing similar information be able to read/write to each other without having to use ArcGIS for Server. Their concern is the amount of traffic on the in-house server and would like to bypass this. The client is also looking into ArcGIS Portal but currently doesn't have it. Also, while AGOL is clearly public, their SQL database is behind a firewall.

Doing some research, it doesn't look likely that this is possible. However, I figured I'd ask here in case anyone has tried to get around this before.

Is there any way around ArcGIS for Server I'm not finding through my own research?

[Update]: With much more clarification, I believe we have come up with a solution that would work. We know Python can interact with AGOL, and we know Python can interact with SQL Server DBs. From here, we just want the two tables to verify each other. Seems much simpler than was previously discussed.


texts:
When segmentizing a geography line in Postgis, the max segment length does not seem to be respected when the line is more than 3000km. The error grows when the length of the line increases.

Below the code to demonstrate it:

-) The printLineSegmentLength method logs the length of each segment of the line.

-) The select create a geography line, then segmentize it and finally call the printLineSegmentLength  method.

SELECT printLineSegmentLength(ST_Segmentize(ST_GeographyFromText('LINESTRING(20 20, 60 60)'), 50000), true);

CREATE OR REPLACE FUNCTION printLineSegmentLength(line geography, use_spheroid boolean)
RETURNS VOID AS $$
DECLARE
    points geometry_dump[] := array(SELECT ST_DumpPoints(line::geometry));
    point_count INTEGER := array_length(points, 1);
    segmentLength NUMERIC;
BEGIN
    FOR i IN 1..point_count - 1 LOOP
        segmentLength := ST_Distance(geography(points[i].geom), geography(points[i + 1].geom), use_spheroid);
        RAISE INFO 'Length: %, %, %', segmentLength, ST_AsText(points[i].geom), ST_AsText(points[i + 1].geom);
        END LOOP;
END;
 $$ LANGUAGE 'plpgsql' IMMUTABLE;

When executing this function, using ST_Distance both with the spheroid or not, the error grows to reach 53263m whereas it shouldn't be more than 50000m.

Changing the length of the segment does not make the error disappear.

Is there anything wrong with my test code?


texts: I know that KML can contain symbology in it. I understand that GML is the open version of it. Can GML contain symbology too?


texts: PostGIS - st_intersects using polygons that cross antimeridan/equator

texts:
New to using PostGIS and have two situations that I'm working with.

In the one case, I have a polygon that crossed the antemeridian and was detecting whether a two separate points located within that polygon one on the west and one on the east side of the antimeridian intersects with the polygon (data types used are geometries). I found I had to use the st_shift_longitude on the geometries passed to the st_intersects function.

The other situation is a polygon crossing the equator with points above and below the equator within the polygon. st_intersects doesn't detect both points intersect the polygon. Do the latitudes need to be shifted in this case? I didn't see a shift latitude function available.

I'm using the following points and polygon at the equator:
st_point(170.0,45.0)
st_point(170.0,-45.0)
and the polygon 'POLYGON((165 40,175 40,175 -50,165 -50,165 40))'

I'd expect that the st_intersect function return true for both points, however, it's only returning true for st_point(170.0,-45.0)


texts: How to change Openlayers MousePosition control as mouse up?

texts:
I am using openlayers 6.x and Mouse position control.

this.control = new MousePosition({
  coordinateFormat: createStringXY(this.precision),
  className: this.className || "mouse-position",
  target: this.target,
  projection: this.projection,
  undefinedHTML: this.undefinedHTML || "&#160;"
});

But this works on mouse move event. But I want to get mouse position on mouse click. So I want to change it.

I read that MousePosition had a function named read(mapEvent). Description is: "Update the projection. Rendering of the coordinates is done in handleMouseMove and handleMouseUp."

But how can I change it?


texts: GeoServer - stores intermittently becoming disabled

texts:
Recently in my company, we are having issues with our GIS web services. The problem is due to the PostGIS Stores in GeoServer have been intermittently been disabled. When we manually enable them again, everything works again like it should.

Does anyone know what could be causing the stores to become disabled?
Is there a solution / work around to either
(a) prevent this from happening; or
(b) fix it automatically when it does happen; or
(c) is there an option to set up an email notification - highlighting that the stores have been disabled and need to be re-enabled?


texts: Road design using Least cost path by longitudinal slope in ArcGIS Desktop?

texts:
I have previously asked for help regarding this issue Finding least cost path in QGIS?, but obviously people do not know exactly the problem related to road design
I used the least cost path, this time in ArcGIS, and got the result as shown below. Other lines are automatically generated with the module least cost path, where the cost raster use slope raster.You can see that the obtained trajectory can not possibly be as roads.

Red lines represent the paths that I hand-drew under rule: In my case a roads may not have longitudinal slope greater than 10%, can be between 1-10%.

What is the longitudinal slope of the terrain?

If the height difference between two contour lines 10m, longitudinal gradient of 10% is the line that connects the two contour lines at 100m, for 8% distance is 125m ...
the formula looks like this: K = E / L-slope

K- step
e - equidistance
L - longitudinal slope of road

Is it possible to use longitudinal slope as raster cost for road design?

If we solve this, we finally resolve the issue of the optimal network paths.


texts: Leaflet controls are not appearing with Promise function

texts:
I am not receiving any errors in the console either.
I am attempting to use the Promise function instead of callback.

var map;
var centerlatlng = L.latLng(38, -98);
var aDiv;
var geojsonfeatures;

$(document).ready(function() {

  var aLayerOne = L.tileLayer('http://{s}.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png', {
    attribution: '&copy; OpenStreetMap, Tiles courtesy of Humanitarian OpenStreetMap Team'
  });

  var aLayerTwo = L.tileLayer('http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}', {
    attribution: 'Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community'
  });

  var aLayerThree = L.tileLayer('http://{s}.tile.openweathermap.org/map/precipitation_cls/{z}/{x}/{y}.png', {
    attribution: 'Map data &copy; OpenWeatherMap',
    opacity: 0.5
  });

  var aLayerFour = L.tileLayer('http://{s}.tile.openweathermap.org/map/temp/{z}/{x}/{y}.png', {
    attribution: 'Map data &copy; OpenWeatherMap',
    opacity: 0.5
  });

  var aLayerFive = L.tileLayer('http://{s}.tile.openweathermap.org/map/pressure_cntr/{z}/{x}/{y}.png', {
    attribution: 'Map data &copy; OpenWeatherMap',
    opacity: 0.5
  });

  var points = L.geoJson(geojsonfeatures, {
    filter: function(feature, layer) {
      return feature.properties.COUNTRY != "UNITED STATES";
    },
    pointToLayer: function(feature, latLng) {
        return L.marker(latLng)
          .bindPopup(feature.properties.NAME);
      } // end of function
  });

  map = L.map('myMap', {
    center: centerlatlng,
    zoom: 4,
    layers: [aLayerOne]
  });

  var promise = $.getJSON("js/weatherstations.js");
  promise.then(function(geojsonfeatures) {
    var weatherstations = L.geoJson(geojsonfeatures);
    points.addTo(map)
    // Add a Layers Control
    var baseMaps = {
      "Open Street Map HOT": aLayerOne,
      "Esri Imagery": aLayerTwo
    };

    var overLays = {
      "Precipitation": aLayerThree,
      "Temperature": aLayerFour,
      "Barometric Pressure": aLayerFive,
      "California Weather Stations": weatherstations
    };

    L.control.layers(baseMaps, overLays).addTo(map);

    // Add a Scale Control
    L.control.scale().addTo(map);

    // Add a Full Screen Control
    // use  "new" in its construct
    new L.Control.FullScreen().addTo(map);

    // Add a Locate Control
    // images folder must be within the same folder as L.Control.Locate.css file
    L.control.locate().addTo(map);

    ///Creating custom control here
    // part 1/3 : Creates a control with the given position
    var aControl = L.control({
      position: 'bottomright'
    });

    // part 2/3 : Should contain code that creates all the neccessary DOM elements for the control
    aControl.onAdd = function() {
      aDiv = L.DomUtil.create('div', 'aCustomC'); // create a div with a class "aCustomC"
      aDiv.innerHTML = "Open-Source Weather Map"; // Each HTML element has an innerHTML property that defines both the HTML code and the text that occurs between that element's opening and closing tag.
      return aDiv;
    }; // end function onAdd

    // part 3/3 : Add the control to the map
    aControl.addTo(map);
    ////////
    var searchControl = new L.esri.Controls.Geosearch().addTo(map);
    var results = new L.LayerGroup().addTo(map);
    searchControl.on('results', function(data) {
      results.clearLayers();
      for (var i = data.results.length - 1; i >= 0; i--) {
        results.addLayer(L.marker(data.results[i].latlng));
      }
    });

    $('.aToolTip').tooltip({
      'container': 'body'
    }); // this works

    $('.aPopOver').popover({
      'container': 'body'
    }); // this works

  });

});


texts: Schema for offline map editing and pass to production environment

texts:
I'm working with QuantumGIs and PostgreSQL/Postgis. So when I edit the layers, the changes are reflected on the DB.

The thing is I want to do some "offline" editing before committing the changes to the DB, as the DB is in production and serves Geoserver for displaying the maps on a website.

Anyone know a good schema to accomplish this? Thanks.


texts: Method to generate points in any projection

texts:
I am writing an application in ArcMap to generate a grid of points. These points have a very defined layout and I have worked out an algorithm that works great under a projected coordinate system. However, as soon as I try it under a geographic coordinate system, my math doesn't stand up as well and I am not an expert in this type of math (but will to learn if pointed in the right direction)

Are there any suggested methods, libraries, or anything else that might make this easier?

A little bit more information about the dots themselves: In a circle of 144 foot radius, the dots will be generated in a hex type pattern. The hex is made up of an equilateral triangle with sides of 9 meters (in meters because ArcMap IPoint X,Y properties are in meters).