texts: How to retain domain descriptions from selected features in database to new database?

texts:
I have the following script that selects features from an enterprise database and then exports them to a new database. However, two columns in the database contain domains. How can I keep the descriptions of the domains in the new database that is created rather than the domain numbers?

import arcpy
from arcpy import env

env.workspace = r'\\path\'
stateFeature = r'\\test.shp'
parish = arcpy.GetParameterAsText(0)
output = arcpy.GetParameterAsText(1)
nameField = "PARISH"

arcpy.CreateFileGDB_management(output, parish)

arcpy.Delete_management('SelectionStateLayer')
arcpy.Delete_management('lyr')

arcpy.MakeFeatureLayer_management(stateFeature, 'SelectionStateLayer', '"' +
str(nameField) + '" = ' + "'" + str(parish) + "'")

for fc in arcpy.ListFeatureClasses():

    if arcpy.Exists('lyr'):
        arcpy.Delete_management('lyr')

    arcpy.MakeFeatureLayer_management(fc, 'lyr')
    arcpy.SelectLayerByLocation_management('lyr', "INTERSECT", 'SelectionStateLayer')

    outname = fc.replace(".","_")
    outpath = output + r'/' + parish + '.gdb' + r'\{0}'.format(outname))
    arcpy.CopyFeatures_management('lyr', outpath)

arcpy.Delete_management('lyr')


texts: Data Source (s) for names and polygon shapes for the world

texts:
I am looking to find out where I can ideally get every country, district, city, town, postcode in the form of polygons.

Is there a central place for this? What are the licences for using the data?


texts:
So I already have the code working to initiate a map and to populate it with data points with a dynamic slider. What I am trying to do now is to have a user type in a query and when they submit that query, the map changes its center lat,lng to the new one.

So it'll pretty much be like

Search Rutgers Club
Get Lat & Lon
Change map center to new Lat & Lon

You can access my live site at parkru.com/v2.html

Also, here is what my js file looks like


texts: Transfer from Android phone to Windows OS truncates lat long

texts: I am clicking GPS enabled pictures of sites which are less than 500 metres apart through Android based Samsung ACE Duos phone. I later use Google Earth on my Laptop  to view these geo tagged pictures. I am facing a problem that Windows Vista OS on my laptop truncates the lat long to six places after decimal. This I realised after seeing the Properties of the transferred JPG file in Laptop. However, in the original picture  the details of the jpg when checked in Gallery of the phone shows 12/15 digits after decimal.


texts: Access MDB connection file result does not have an object ID in ArcGIS 10.8

texts:
I have made an MDB Access Connection by linking the Excel file with the Access one, but when I looked at the result in ArcMap the feature class table does not have an Object ID on it. I tried to export the file in shapefile in order to get Object ID on it.
How can we keep the connection with the MDB Access connection with the new feature class in order to update the data directly from Excel table to be displayed and updated directly on ArcMap?


texts: Calculating features from one vector layer in another layer (grid)?

texts:
I have two vector layers

a landcover map showing the different habitats in a region
a grid of 4 x 4 km squares that covers the region

I have to calculate the percentage of landcover in each grid square, so the percentage of what is in the landcover layer for each individual 4 x 4 km square.

How to achieve this in QGIS?


texts: Is the half installation issue just related to my desktop?

texts: I have tried to install QGIS version 3.16 on to my work computer (I am working remotely from home) and it asked me to uninstall a programme which was apparently already on the system (I think it may have been QGIS 3.16). I clicked yes, but halfway through I was concerned that I may be deleting colleague files, so I stopped the process. Now I am unable to open some QGIS based files (I think they may have only been half installed) and I'm just wondering if there's a way to rectify this, and whether this will have affected my colleagues QGIS software too?


texts: Looking for geodata of Europe's major roads

texts: As stated in the title, I'm currently looking for a shapefile/other geodata source of Europe's major roads (highways). If anybody knows a free datasource (or a possibility to extract that data from Open Street Map without processing the world file), it would be greatly appreciated.


texts: I have combined MODIS blue marble images into a mosaic, a screen shot of the mosaic raster image with its source info opened in ArcMap is attached. The  file source shows that it has angular coordinates defined. I tried to define projection/project but nothing seems to be working. My problem is that this image does not align with other layers having projections defined (e.g. basemap). So I need to align it with other global maps (e.g. basemap) so that it overlaps with them. I know the coordinates of all the four corners of the map (e.g. bottom left corner is -180, -90 or 180W, 90S). I thought it should be a simple problem of georeferencing but it does not seem so after searching all the information online. How can I project or georeference this image in a simple way using the tools in ArcGIS? Please help.


texts: I want to have a raster representation of the river, so that I could further analyse combining other raster data. I tried extract the DEM, polygon to raster and resample. But the results are not good enough: missing raster cells in some narrow parts. Is there a way to fix it? I have spatial analyst in Arcgis 10.


texts:
has anyone a suggestion which segmentation-tool is best suited for this task.
I'm looking for a fast tool to do this approx 1km²/day the faster the better of course.

I've tried Orfeo's "meanshift" and "exact large meanshift" as well as "watershed", "meanshift" has offered me the best results so far but the whole process is about 5 days for 0.6km² so i'm looking for a faster way.
spatial resolution is 5-10cm

It doesn't matter if it's with orfeo, saga, grass or just python.


texts: Why is controlling shell script not working?

texts:
I updated my post after I resolved the previous errors, but additional errors were produced.

I want to create a controlling shell script that contains water.mapcalc within it. The problem is, I do not know how to save a code in a file called water.mapcalc. I'm trying out the hydrologic simulation.

It is stated in page 12:
"For those wishing to repeat the experiment with GRASS enter the model
code in a file called water.mapcalc then create a controlling shell script."

I understood the instruction as I should save the following in a raster map called water.mapcalc:

r.mapcalc water.mapcalc = "water = water + eval(x=elev + water,
    if(x>(y=elev[-1,0]+water[-1,0]),
        -.15*if(elev>y, water, x-y),
        .15*if(elev[-1,0]>x, water[-1,0], y-x))+

    if(x>(y=elev[1,0]+water[1,0]),
        -.15*if(elev>y, water, x-y),
        .15*if(elev[1,0]>x, water[1,0], y-x))+

    if(x>(y=elev[0,-1]+water[0,-1]),
        -.15*if(elev>y, water, x-y),
        .15*if(elev[0,-1]>x, water[0,-1], y-x))+

    if(x>(y=elev[0,1]+water[0,1]),
        -.15*if(elev>y, water, x-y),
        .15*if(elev[0,1]>x, water[0,1], y-x))+

    if(x>(y=elev[-1,1]+water[-1,1]),
        -.10*if(elev>y, water, x-y),
        .10*if(elev[-1,1]>x, water[-1,1], y-x))+

    if(x>(y=elev[1,1]+water[1,1]),
        -.10*if(elev>y, water, x-y),
        .10*if(elev[1,1]>x, water[1,1], y-x))+

    if(x>(y=elev[1,-1]+water[1,-1]),
        -.10*if(elev>y, water, x-y),
        .10*if(elev[1,-1]>x, water[1,-1], y-x))+

    if(x>(y=elev[-1,-1]+water[-1,-1]),
        -.10*if(elev>y, water, x-y),
        .10*if(elev[-1,-1]>x, water[-1,-1], y-x)))"

Should I save the model code in a water.mapcalc raster map as what I did above? or should I save a water.mapcalc file in the /bin folder? I tried both but it still produced errors.

I already did the said map (water.mapcalc).
and I created the following shell controlling script that needs the file (water.mapcalc). This is the shell script (myscript5.sh) and I saved it in the \bin folder of GRASS:

#!/bin/sh
export GRASS_PNGFILE=/tmp/grass6output.png

r.mapcalc water = 120
d.rast water

r.mapcalc < water.mapcalc
d.rast water

r.mapcalc < water.mapcalc
d.rast water

g.copy rast=water,water.1

r.mapcalc < water.mapcalc
d.rast water

r.mapcalc < water.mapcalc
d.rast water

g.copy rast=water,water.2

Then, I run the sh myscript5.sh in the CMD of GRASS and it produced this message:

PNG: GRASS_TRUECOLOR status: TRUE
PNG: collecting to file: C:/Users/Rhen/AppData/Local/Temp/grass6output.png,
GRASS_WIDTH=768, GRASS_HEIGHT=508
/c/Program Files/GRASS-64/bin/myscript5.sh: line 7:
water.mapcalc: No such file or directory
PNG: GRASS_TRUECOLOR status: TRUE
PNG: collecting to file: C:/Users/Rhen/AppData/Local/Temp/grass6output.png,
GRASS_WIDTH=768, GRASS_HEIGHT=508
/c/Program Files/GRASS-64/bin/myscript5.sh: line 10:
water.mapcalc: No such file or directory
PNG: GRASS_TRUECOLOR status: TRUE
PNG: collecting to file: C:/Users/Rhen/AppData/Local/Temp/grass6output.png,
GRASS_WIDTH=768, GRASS_HEIGHT=508
Copy raster <water@PERMANENT> to current mapset as <water.1>
/c/Program Files/GRASS-64/bin/myscript5.sh: line 15:
water.mapcalc: No such file or directory
PNG: GRASS_TRUECOLOR status: TRUE
PNG: collecting to file: C:/Users/Rhen/AppData/Local/Temp/grass6output.png,
GRASS_WIDTH=768, GRASS_HEIGHT=508
/c/Program Files/GRASS-64/bin/myscript5.sh: line 18:
water.mapcalc: No such file or directory
PNG: GRASS_TRUECOLOR status: TRUE
PNG: collecting to file: C:/Users/Rhen/AppData/Local/Temp/grass6output.png,
GRASS_WIDTH=768, GRASS_HEIGHT=508
Copy raster <water@PERMANENT> to current mapset as <water.2>

the error is: water.mapcalc: No such file or directory. I think there is a problem on how I saved the water.mapcalc file. I do not know what format (raster/.txt/.sh file) and where will I save it (should I also save it in the /bin folder together with my script). I tried all of those and it still did not work. I cannot figure out where the problem is.


texts: Adding M coordinate with different measurings in PostGIS

texts:
I have a line table made of several linestring features in succession, so that all linestrings together would form one continuous line, such as:

id  |  geom
----+-----------------------
1   |  LINESTRING(0 0, 1 1)
2   |  LINESTRING(1 1, 2 2)
3   |  LINESTRING(2 2, 3 3)

I want to add measuring to this line. Normally I'd merge all line segments into one and do ST_AddMeasure(), but in this case each segment should have different measurings (e.g. the segment (0 0, 1 1) would go from 0 to 10, segment (1 1, 2 2) would go from 10 to 12, and segment (2 2, 3 3) would go from 12 to 60).

In the line table I have columns with the start and end measures for each line segment, so generating different measurings separately is easy enough. The problem comes from adding them all together in one continuous linestring. ST_Union() drops the M-index, same with ST_UnaryUnion() and ST_MemUnion().

Hence, my question: is there a way to merge line segments without droping the M-coordinate? Alternatively, is there a way to add measure to a linestring, taking into consideration difference in measurings along the line?


texts: Using WFST in Leaflet: can not move/delete features

texts:
I write a Leaflet app which makes use of two plugins Leaflet.Draw and Leaflet.WFST. The first allows to create/edit/delete features, the second saves changes to my PostGIS database hosted by GeoServer.

Here is the code:

// Basemap
var map = new L.map('map').setView([57.704, 11.965], 6);
L.esri.basemapLayer('Topographic').addTo(map);

// Layer to draw in
var wfst = new L.WFST({
    url: 'http://localhost:8080/geoserver/ows',
    typeNS: 'drawing',
    typeName: 'Stockholm_points',
    crs: L.CRS.EPSG4326,
    geometryField: 'the_geom',
    showExisting: true,
}).addTo(map);

// Drawing
var drawControl = new L.Control.Draw({
    draw:{circle:false, circlemarker:false, rectangle:false,
          },
    edit:{featureGroup: wfst } });
map.addControl(drawControl);

map.on('draw:created', function (e) {
    var layer = e.layer;
    wfst.addLayer(layer)});

map.on('draw:edited', function (e) {
    var layers = e.layers;
    layers.eachLayer( function (layer) {
        wfst.editLayer(layer);
        });
});

// Save button
L.easyButton('fa-save', function () {
         wfst.save();
     }, 'Save changes').addTo(map);

It makes good job when I create a feature.
However, when it comes to editing of existing features (for example, I want to drag a marker to another location or delete several markers), my changes don't get saved. Sometimes one of "deleted/edited" markers disappear/change location (as I would expect it to do), but never all them.

In GeoServer's log I see that my app communicates with server and sends requests to update features.. but with a little success.

What is wrong?

Here is a snippet from logfile. I tried to change location for three features and save edits. As a result only one feature appeared to be updated when I refresh my app.

2018-08-17 12:30:19,513 INFO [geoserver.wfs] -
Request: transaction
    service = WFS
    version = 1.1.0
    baseUrl = http://localhost:8080/geoserver/
    group[0] = wfs:update=net.opengis.wfs.impl.UpdateElementTypeImpl@12d70e9 (filter: [ Stockholm_points.2 ], handle: null, inputFormat: <unset>, srsName: null, typeName: {http://geoserver.org/drawing}Stockholm_points)group[0] = wfs:update=net.opengis.wfs.impl.UpdateElementTypeImpl@12a9a36 (filter: [ Stockholm_points.3 ], handle: null, inputFormat: <unset>, srsName: null, typeName: {http://geoserver.org/drawing}Stockholm_points)group[0] = wfs:update=net.opengis.wfs.impl.UpdateElementTypeImpl@1d7334b (filter: [ Stockholm_points.4 ], handle: null, inputFormat: <unset>, srsName: null, typeName: {http://geoserver.org/drawing}Stockholm_points)
    update[0]:
        property[0]:
            name = {http://geoserver.org/drawing}the_geom
            value = POINT (10.832519531250002 59.91097597079679)
        filter = [ Stockholm_points.2 ]
        inputFormat = x-application/gml:3
        typeName = {http://geoserver.org/drawing}Stockholm_pointsupdate[0]:
        property[0]:
            name = {http://geoserver.org/drawing}the_geom
            value = POINT (11.206054687500002 59.66774058164963)
        filter = [ Stockholm_points.3 ]
        inputFormat = x-application/gml:3
        typeName = {http://geoserver.org/drawing}Stockholm_pointsupdate[0]:
        property[0]:
            name = {http://geoserver.org/drawing}the_geom
            value = POINT (10.173339843750002 59.567723306212955)
        filter = [ Stockholm_points.4 ]
        inputFormat = x-application/gml:3
        typeName = {http://geoserver.org/drawing}Stockholm_points
    releaseAction = ALL
2018-08-17 12:30:19,536 INFO [geoserver.wfs] -
Request: getServiceInfo
2018-08-17 12:30:19,550 INFO [geoserver.wfs] -
Request: getFeature
    service = WFS
    version = 1.1.0
    baseUrl = http://localhost:8080/geoserver/
    query[0]:
        filter = [ none ]
        srsName = EPSG:4326
        typeName[0] = {http://geoserver.org/drawing}Stockholm_points
    outputFormat = text/xml; subtype=gml/3.1.1
    resultType = results


texts: GeoServer WPS built-in functions vs. GDAL utilities

texts:
I have been doing a series of back-end GIS operations such as:

Clipping a raster file by using user's inputs (i.e. bounding rectangle)
Returning a styled WMS based on the new raster file.

By using GDAL utilities for the clipping and Curl within PHP for the GeoServer response, and OpenLayers for the front-end.

The above works fine.

But now I would like to turn that into a simpler architecture, so is it possible to use WPS for creating the clipped raster in the server and immediately generating a WMS based on that. In other words, can I do the whole process within GeoServer and without using GDAL utlities?

I am new to the WPS aproach, so far I haven't found how to do the clipping, creating the workspace, the datastore, the layer and the styling, and publish a new WMS back to the user, all in one go.


texts: Python Script to JSON as Geo-processing service

texts:
I have been having some trouble with this script that I created. Ideally, I would like to query between to a layer file and feature class to find certain features in the class withing a certain distance from the layer file.

Then take certain attributes from those selections and store them in a python dictionary which will be converted into a JSON format and uploaded as a geoprocessing service.

I have excluded the remaining code because this bulk is what seems to be dis-functioning.

As a standalone script, it works fine and gives me the information intended. As a service -- here is the error received

import arcpy
import sys, os
import json
arcpy.env.overwriteOutput = True
segid = arcpy.GetParameter(0)                                                                                connection1 = r'Database Connections\DDOTLRS as LRSVIEWER.sde'
connection2 = r"Database Connections\Connection to DDOTGIS as DDOTEDITOR.sde"
supports = os.path.join(connection2, 'TOA.Supports')
                                                                                def getFeatureOIDsNearSeg(selectthesefeatures = None, segid = None , distance=150):
    '''  Function to spatially query for features within a certain distance of a street segment.
    user has option to provide distance '''                                                                                streetsegs = os.path.join(connection1,'Ddotlrs.RH.LRSE_StreetSegment')                         streetsegs_lyr = "streetsegs_lyr"                                                                   selectthesefeatures_lyr = "selectthesefeatures_lyr"

try:
    arcpy.MakeFeatureLayer_management(streetsegs, streetsegs_lyr, "STREETSEGID = " + str(segid)) #only selecting this block street seg ID
except ValueError:
    print "Error.  Cannot find this StreetSeg ID.  Check ID or source (currently set to {0})".format(segid)

try:
    arcpy.SelectLayerByLocation_management(in_layer=selectthesefeatures_lyr,
                                           overlap_type="WITHIN_A_DISTANCE",
                                           select_features=streetsegs_lyr,
                                           search_distance=str(distance) + " FEET",
                                           selection_type="NEW_SELECTION")
    #need to find out what was selected...
    selectthesefeatures_lyr_sel_desc = arcpy.Describe(selectthesefeatures_lyr)
    selectthesefeatures_lyr_sel_FIDS = selectthesefeatures_lyr_sel_desc.FIDset.split(";")
    selectthesefeatures_selected = ','.join(map(str, selectthesefeatures_lyr_sel_FIDS)) #capture all the IDs for the selected supports
    OIDs = "OBJECTID in ( " + selectthesefeatures_selected + ")"
    return OIDs
except ValueError:
    print "Error.  Problem selecting supports along this segment {0}.".format(segid)
return OIDs


texts: Exporting large high resolution map to PDF in ArcMap errors?

texts:
I am attempting to export a large E size map at 300 DPI in ArcMap 10.2.2 but each time I export there are errors on the map.  Once the PDF is generated there are areas where symbols are stretched out in a line across the width of the map.  I have tried it on ArcMap 10.1 and 10.2.2 and I get the same errors.  I know I can avoid these errors if I export at a lower DPI, but I want a really high resolution print of this particular map to be prominently displayed.  All other parameters for the export to PDF are the default values.

Here is an example of the errors:


texts: Find the GPS coordinates where two boundary lines intersect

texts:
I have a detailed topo map for an irregular shaped property. It has boundary lines expressed like this, for example; S76˚45'08"W and S59˚38'24"W. The map also notes how long each boundary line is. I need a beginning reference point in GPS coordinates to mark a physical point on the property (using my cell phone's compass feature), and then need to use the info on the map to stake out the property border.

How do I calculate the GPS coordinates from where these two lines intersect?


texts: Can't solve "OperationNotSupported" error on WFS GetFeature request

texts:
I'm working on a web map viewer using OpenLayers, ExtJS, GeoExt and GXP. I've added a LayerTree (gxp_layertree) that shows layers from various WMS servers from different GeoServer versions. I've also configured a FeatureManager (gxp_featuremanager) with 'autoLoadFeatures = true' and a FeatureGrid (gxp_featuregrid).

My problem appears when I click on a layer tree node and a WFS GetFeature request fires: the response is always "OperationNotSupported" "No such operation WFS null WfsDispatcher".

I've checked the XML body content with the embedded 'demo requests' in GeoServer and it's OK. That is:

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" maxFeatures="100" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <wfs:Query typeName="feature:mce_prov_escuelas" srsName="EPSG:4326" xmlns:feature="mce">
    <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
      <ogc:BBOX>
         <ogc:PropertyName>the_geom</ogc:PropertyName>
         <gml:Envelope xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
            <gml:lowerCorner>-62.406926496825 -26.770412516161</gml:lowerCorner>
            <gml:upperCorner>-57.557744198657 -22.49170214764</gml:upperCorner>
         </gml:Envelope>
      </ogc:BBOX>
     </ogc:Filter>
  </wfs:Query>
</wfs:GetFeature>

It's a POST request and here's a copy of the headers (from Firebug):

REQUEST HEADERS
Accept      text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8
Accept-Encoding gzip, deflate
Accept-Language es-ar,es;q=0.8,en-us;q=0.5,en;q=0.3
Content-Length  723
Content-Type    application/xml; charset=UTF-8
Host        localhost
Referer     http://localhost/visor/index.html
User-Agent      Mozilla/5.0 (X11; Ubuntu; Linux i686; rv:24.0) Gecko/20100101 Firefox/24.0
X-Requested-With    XMLHttpRequest

Any clue?

I've tried everything I know and Googled it without many results.


texts: ArcGIS Server - problem publishing python script that uses urllib

texts:
I have problem publishing this simple python script as geoprocessing service with ArcGIS Server.

import urllib
import arcpy

u = arcpy.GetParameter(0)
url = "http://geoegl.msp.gouv.qc.ca/cgi-wms/mapcache.fcgi/tms/1.0.0/carte_gouv_qc_public@EPSG_3857/" + u

image = urllib.urlopen(url).read()

arcpy.SetParameter(1,image)

This is a simple proxy that uses urllib to redirect images. When I try to publish this script, ArcMap always go idle for a long long time. I think the cause of the problem is urllib. If I try to publish a script that do not import urllib, I have no problem.

After a long idle time, the Publisher returns the following errors:

High    Unresolved    00178    Data: C:\Windows\bfsvc.exe, used by Script Script cannot be copied to the server

High    Unresolved    24032    Data source used by Script Script is not registered with the server and will be copied to the server: C:\temp\arcgis\arcpy\

Can someone tell me why trying to publish a python script that import urllib into ArcGIS Server could return those errors?


texts: Merging two multipolygon tables into one and afterwards dissolving boundaries using PostGIS

texts:
I am trying to merge two multipolygon tables into one and at the same time - in the same query - dissolve the boundaries between the polygon districts using postgis.
The two tables are the red polygons and the beige polygons in the map below.
As mentioned, the first thing I want to do is to merge the two tables into one and second I want to dissolve the boundaries based on values in a column name. Name are e.g. 'JS22C' and 'JS22A' in the map. I cannot figure out if you do this by doing a join operation or by using a UNION/ST_UNION

I have been trying the following:
CREATE OR REPLACE VIEW schema.tablename AS
SELECT row_number() over() as id, st_union(a.the_geom, b.the_geom), name
FROM schema.table1 a, schema.table2 b
GROUP BY name;

But it says that name is ambiguous which is quite obvious because the column name exists in both tables .
Can anybody help out?


texts: How to get distance along line at specific point in QGIS/PyQGIS

texts:
I have a line geometry and a point geometry. I need to know the distance along the line at the point. I can see the QgsGeometry.interpolate() function but that's the reverse of what I want to do. I have a point, I need the function to return the distance. Is this feature missing in the API?

For those versed in PostGIS, I need something like the ST_Line_Locate_Point() function but in PyQGIS.


texts: Efficiently working with converting vrt rasters to greyscale in QGIS?

texts:
I have large raster tile sets working well in QGIS (Ordnance Survey UK data). Sometimes I want to grey/fade these to make my data clearer to see. Attempting to access properties on QGIS for the VRT effectively crashes QGIS for the larger tile sets (it works for the smaller tile set, but takes a long time... but for the layers with more tiles I've not waited long enough to see if it works because more than a minute or ten becomes silly).

Can I access the properties in another way - for instance by taking a text editor to something (the project file presumably)? Should I do this then display the layers twice (once colour and once grey) so they can be turned on and off at will?


texts: Shape layer turns non-editable if saved filtered

texts:
I am working in QGIS 2.18.3.

I have a shape layer with point geometry, which I filter by some field, and then I work on its attribute table. Once I am done for the day, I save the project so I can continue the next day.

I store data on my computer, where nobody else connects to it, I do close the project afterwards, and finally I do switch off my computer.

Here is my shapefile with no filter (can be edited),

then my shapefile after I applied the filter (still can be edited).

But every day, when I open the project, I find the shapefile can not be edited, so I have to un-filter it and then save it with another name, delete de previous one and keep on working on the new one which allows me now to keep on editing.

I have tried duplicating the layer, but the duplicate is still non-editable.

Why does this happen? Is it a consequence of saving the project with a filtered shapefile?


texts:
I have a RasterBrick and I am trying to add new layers to it, (with addLayer) but every time I try I receive the following message:

  Error in sapply(x, fromDisk) & sapply(x, inMemory) :    operations are
  possible only for numeric, logical or complex types

The brick and the layers have the same crs, res, and extent, and I reprojected the brick to match the layers.

Is there anything else that might be keeping the two from merging?

The only thing I can think of is that the brick is stored on disk and the layer is stored on memory- but would that actually effect anything?


texts: Unit of Google Maps API geometry services

texts:
When i call the method

google.maps.geometry.spherical.computeDistanceBetween(from, to);

it returns a number,

is this distance in KiloMeters or Miles??

How can we find this and doc on this API?


texts: Using Collect Values Tool with Union in ModelBuilder?

texts:
I would like to feed a list of feature classes (only one polygon per feature class, but up to 150 feature classes) into the Union tool in ModelBuilder by using a Feature Class Iterator and the Collect Values Tool.

Unfortunately, it does not work since "Collect Values" only works with Merge, Append, etc.

Does anyone know a workaround or another way to do this?


texts: options for using "identify" command with public api

texts:
Disclaimer: I'm new to GIS

Goal: Create a report of how many customers we have in each senate district.

Progress: I found this utility from the state of Missouri which allows me to search for an address and it displays the corresponding district.  I did some debugging and found that it looks like it is doing two calls to an ArcGIS server (Using Branson, MO as example):

https://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer/findAddressCandidates?SingleLine=branson, mo&f=json&outSR={"wkid":102100,"latestWkid":3857}&outFields=*
https://ogi.oa.mo.gov/ArcGIS/rest/services/LEGIS/legisDistrict/MapServer/identify?f=json&geometry={"x":-10377036.95984281,"y":4389554.630050001,"spatialReference":{"wkid":102100,"latestWkid":3857}}&tolerance=1&returnGeometry=true&mapExtent={"xmin":-10377151.615385238,"ymin":4389439.974507573,"xmax":-10376922.304300383,"ymax":4389669.285592428,"spatialReference":{"wkid":102100}}&imageDisplay=1599,845,96&geometryType=esriGeometryPoint&sr=102100&layers=visible:1

The first call returns a list of candidate objects like this:

{
    "spatialReference": {
        "wkid": 102100,
        "latestWkid": 3857
    },
    "candidates": [
        {
            "address": "Branson, Missouri, United States",
            "location": {
                "x": -10377036.95984281,
                "y": 4389554.630050001
            },
            "score": 100,
            "attributes": {...},
            "extent": {
                "xmin": -10383716.235155277,
                "ymin": 4381233.546051718,
                "xmax": -10370357.896260085,
                "ymax": 4397882.30382274
            }
        }
    ]
}

The second call returns a list of result objects:

{
    "results": [
        {
            "layerId": 1,
            "layerName": "Senate District",
            "displayFieldName": "District",
            "value": "29",
            "attributes": {
                "OBJECTID": "29",
                "District": "29",
                "Population": "181191",
                "webLink": "http://www.senate.mo.gov/15info/members/mem29.htm",
                "photoLink": "http://www.senate.mo.gov/15info/graphics/d29-photo.gif",
                "Name": "David Sater",
                "Shape": "Polygon",
                "District_Number": "29",
                "Shape.area": "8040296221.235743",
                "Shape.len": "505379.224253"
            },
            "geometryType": "esriGeometryPolygon",
            "geometry": {}
        }
    ]
}

attributes.district being the value I want.

As you can see above, the second call sends a parameter called geometry which is a combination of location and spatialReference, as well as mapExtents which is a combination of extent and spatialReference.  It also sends some other things.

The point is, this seems to require a lot of object/string manipulation to do both calls.  The question is, are both calls necessary?  Is there another way to call /identify, perhaps directly with latlng which I already have?  I've looked at the documentation and I don't see anything, but to be honest I don't even understand what mapExtent is or why I need it.  Why would it be different for every point I query?


texts: Creating lines that connect the centroids of adjacent polygons to one another

texts: I have about 10 polygons, and I have added their centroids to the map. I want to create lines connecting the centroids of each adjacent polygon. If polygon 1 borders polygon 2 and 3, I would like a line going from 1 to 2 and from 1 to 3. Hopefully this makes sense, and someone can help me out. Also, I would like to indicate the distances between the centroids.


texts:
I need to plot the incline of roads in England. I could estimate this using elevation data, however this would be a significant amount of work.
Does this data already exists (ideally free with a permissive license)?
I've spent this morning googling and haven't been able to find anything, but thought there is a good chance a dataset like this already exists, so thought I should check before embarking on creating it myself.


texts: Randomized organic warp of adjacent polygons that preserves topology and fixed nodes

texts:
I'm looking for algoritms to modify a set of polygons to make them appear more organic. They still need to maintain their topology and furthermore some of the nodes might be fixed. Please see the example image below where the left hexagons needs be transformed to something like shown on the right.

Note the red fixed nodes, those are not to be moved.

I do know how to "randomize" a line, but not how to take topology and bounds into account as its likely the polygons will intersect if not careful. Any ideas on how to achieve this? I would love to know the names of good algorithms to use, or even better Python libraries or example code.


texts: Creating type of ranking map on ArcGIS Desktop?

texts: I wanted to create a ranking continuous legend ...but everytime symbology come in class interval with different box of colour ramp .


texts:
I am using PostgreSQL 9.6 and PostGIS extension 2.4.4. On my machine, I get an error while running ST_Transform query:

SELECT ST_Transform(ST_SetSRID(ST_MakePoint('27.9252', '-26.0987'),4326),32619)

returns

transform: couldn't project point (27.9252 -26.0987 0): latitude or longitude exceeded limits (-14)

But when I run it on our testing/prod environment hosted in AWS, I don't get the error message and it suceeds with result

POINT (9591132.411356987 -11522592.397666786)

Locally, I am using Mint 18.3 and it is configured in French.

I could not figure out configuration difference between my local environment and our hosted environment.

How can I get this query to be run without error ?

SELECT PostGIS_Full_Version() returns:

Locally

POSTGIS="2.4.4 r16526" PGSQL="96" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" LIBPROTOBUF="1.2.1" TOPOLOGY RASTER
Hosted environment

POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d6" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.1.4, released 2017/06/23" LIBXML="2.9.1" LIBJSON="0.12" RASTER


texts:
OpenLayers map is requesting tiles out of range and all the map is pink.
If I inspect single tile, the address is next:

  /Tiles/15/16384/16384.png

But the folder with generated tiles does not contain such values:
The range is only: /Tiles/15/16099 - /Tiles/15/16119

The tiles are generated by Maperitive.
Is it a Maperative problem or I should set additional options to OpenLayers.Map?

map = new OpenLayers.Map("demoMap");
var cLayer2 = new OpenLayers.Layer.OSM("New Layer",
"http://localhost:15157/Map/Tiles/${z}/${x}/${y}.png", { numZoomLevels: 19 });
map.addLayer(cLayer2);
map.zoomToMaxExtent();


texts: How to Resolve Inconsistent OpenLayers vector position Over Google map

texts:
I have an OpenLayers map with a Google base layer and an OpenLayers vector layer on which people can draw their own bounding boxes using the DrawFeature control.

I have found that, as you zoom in, the vector is not plotted covering the same co-ordinates at each zoom level. Rather it seems to zoom at a different ratio to the underlying layer.

Strangely, if I resize the browser window, the vector then stays at a consistent position through all the zoom levels (i.e. becoming larger or smaller but covering the same geographic range, as one would want it to), but slightly offset from the position it was initially drawn at. The remains the case even if I make the browser window full size again.

See example below:

I draw a vector touching the Belgian border.

I zoom out and the vector is now well within Germany.

I resize the browser window and the box anchors to a consistent set of co-ordinates but offset from its initial position, now overlapping the Belgian border

Please note that this is not the same issue as the vectors scrolling more slowly than the map when dragging the map. The vectors seem to stick to the same area when dragged, but zooming causes them to move relative to the underlying Google layer.

Is this a known issue and can anyone help resolve?

In answer to the question about my initialisation of the layers, my code is as follows:

var miny=-2000000; var minx=6000000; var maxy=600000; var maxx=10003750;

    // allow testing of specific renderers via "?renderer=Canvas", etc
var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;

//define the layers needed
    g_phys  =   new OpenLayers.Layer.Google("Google Physical",{type: google.maps.MapTypeId.TERRAIN, numZoomLevels:20}, {isBaseLayer:true});
    g_sat   =   new OpenLayers.Layer.Google("Google Satellite",{type: google.maps.MapTypeId.SATELLITE, numZoomLevels:20},{isBaseLayer:true});
    bboxes  =   new OpenLayers.Layer.Vector("Bounding Boxes", {rederers:renderer});

//delete feature on select
bboxes.events.on({"featureselected":function(e){bboxes.removeFeatures(e.feature); erasefeature(e);}});
//refresh on visibility change
bboxes.events.on({"visibilitychanged":refreshMap});
//do stuff on modification of vector
bboxes.events.on({"afterfeaturemodified":modbox});

//create new set of option for the map
var options = {
            projection: new OpenLayers.Projection("EPSG:3857"),
            units: "m",
            maxResolution: "auto",
            maxExtent: new OpenLayers.Bounds(miny, minx, maxy, maxx),
            theme: null,
            layers: [
                //define the base layers for the map
                g_phys,
                g_sat,
                bboxes

                    ]

        };

//create new map object based on "map" div
map=new OpenLayers.Map('map',options);


texts: Trouble clipping two vector layers (EPSG is the same in both)

texts:
I am trying to clip a layer but it gives empty output (No error message, just an empty layer)

The desired input layer is a hexagonal grid (148980 features) created in MMQGIS in EPSG: 32629.

The clipping layer is a layer that I created myself, by merging two line vectors and then using the tool Lines to Polygons, also in EPSG: 32629, which you can see below.

Both layers are quite large, since they represent a portion of the Exclusive Economic zone of Portugal (around 100 000km2)

I'm working in QGIS 2.18.16


texts: GeoServer - Layer Renders Differently based on WMS vs. GWC

texts:
I have a streams layer that renders with a label.

I've noticed that the label renders differently if I use a WMS GetMap request, vs. if I use a GeoWebCache (GWC) request.

Here is an example:

This is a WMS GetMap request result:

This is a GWC request result:

Notice how the Merced River label is positioned differently the two images?

These are both the same layer, using the same style.

This is very weird.  I'd think that they'd render the same, regardless of if it's already been cached.  It's a problem for my workflow, which combines WMS with GWC and relies on the labels rendering in the same place.

@Son of a Beach asked why I need to render using both GWC and WMS.  GWC is so much faster if they tiles are pre-baked.  However, we want to be able to select a stream by clicking on it, and have it be in a "selected" state, i.e. render with a different style.  So I just have a second layer on the map, e.g. streams_selected, which includes a filter with its WMS Get Map request.  The filter is for a particular feature ID.  It renders on top of the GWC (unselected) layer.  I need the labels to line up so it doesn't look funky.


texts: ArcPy Calculate Field how to convert distances from decimal degree to meters

texts:
How do I use arcpy.CalculateField_management() to convert a column with distances in degrees to meters? I don't see this specific example in the docs. This post, answer by "dlee-esristaff" implies it can be done in two steps but I am having trouble implementing it.
Specifically, I don't see how to specify my source field (DISTANCE) and target field (dist_m) in the arcpy.CalculateField_management() inputs. Do I need to write a short function to do this, and then call the function? (If so, I am not sure how to write the function, I think !shape.geodesicLength@meters! is used but I am not sure how.)
Using ArcMap 10.7 Advanced license.
Also, please note that everything is done in WGS84, instead of UTM (hence distances in degrees) because the source SHPs are from all over the world and we don't have local UTM EPSGs. The input SHP here is a point SHP which is output from the Point Distance tool.
arcpy.AddField_management(myshp, 'dist_m', "FLOAT")

# expression = ??? somehow get !shape.geodesicLength@meters! from ['DISTANCE'] into ['dist_m']
# codeblock = ??? somehow get !shape.geodesicLength@meters! from ['DISTANCE'] into ['dist_m']
# ??? arcpy.CalculateField_management(myshp, ['dist_m'], expression, "PYTHON3", codeblock)


texts: Geoserver: SQL-View for layer creation ignores PostgreSQL array_agg function

texts:
I am running Geoserver 2.11.2 running on Tomcat 7 on a Win7 64bit system, PostgreSQL 9.6.

I am trying to create a parametric view for layer creation like this

SELECT
%FOO% AS thing,
array_agg(id) AS ids,
array_agg(name) AS names
FROM
settings
WHERE
%FOO% IS NOT NULL
GROUP BY
%FOO%

However the only attribute showing up is 'thing'.
When using another aggregate function (for example string_agg) instead of array_agg like this

SELECT
    %FOO% AS thing,
    string_agg(id,',') AS ids,
    string_agg(name,',') AS names
    FROM
    settings
    WHERE
    %FOO% IS NOT NULL
    GROUP BY
    %FOO%

everything works fine with all attributes (thing,ids,names) showing up as expected. What is going on here? Does the term 'array' have to be escaped/aliased in some sort of way?


texts: How to access CRS of shapefile in standalone script?

texts:
I'm trying to access the CRS of a layer and use this to create a grid. I followed the answer in this post which works great in QGIS but I would like it to work in a standalone script. Here is a snippet:

for fname in glob.glob("C:\Users\moi\Desktop\Test\Grid\\" + "*.shp"):

    extent = QgsVectorLayer( fname, '', 'ogr' ).extent()
    centerx = (extent.xMinimum() + extent.xMaximum()) / 2
    centery = (extent.yMinimum() + extent.yMaximum()) / 2
    width = extent.xMaximum() - extent.xMinimum()
    height = extent.yMaximum() - extent.yMinimum()
    CRS = QgsMapLayer.crs(fname).authid()

    general.runalg("qgis:creategrid", 1000, 1000, width, height, centerx, centery, 1, CRS, None)

Running the above gives me the following error message:

  TypeError: QgsMapLayer.crs(): first argument of unbound method must have type 'QgsMapLayer'


texts:
I am interested in producing input files for SWMM from GIS data, hopefully from QGIS. Does any one know a free tool to do it?

I found this plug-in qgis-swmm, but I do not see it in the repository of plugins. Does it work?


texts: How to Calculate Custom Zonal Statistics?

texts: I have a vector polygon shapefile outlining counties for which I have a raster of nighttime-luminosity data. I've been able to compute Zonal Statistics in ArcGIS 10.0 using Zonal Statistics as Table and simply selecting the polygons as the first feature class option. However, I'd like to get some more stats than those that are available with the pre-packaged Zonal Stats. Specifically, I'd like to get counts and values -- for instance, the number of pixels with a given luminosity value in each polygon. Does anyone know if this is possible? If not in ArcGIS, in QGIS, perhaps? Thanks.


texts: Displaying lines in QGIS from ST_MakeLine() not working with UPDATE, but works with CREATE TABLE

texts:
When using UPDATE to create geometries based on the SCHNUM (school number) of students to their attend school, not all the lines show up in QGIS. However, when the basic query result is written to a new table, all the lines display just fine.

Without UPDATE, this query works just fine (wrapped in CREATE TABLE):

 select
      SCH.schnum
   , a.studentnumber
      , ST_MakeLine(a.geom,
                    sch.geom) as geom_new
    from
      "public"."students" as a
      join dpsdata."Schools_Current" as sch
        on  a.schoolnum = sch.schnum

Here is the result, which is what I expect:

When using the UPDATE function, which sets the GEOM_LINE column in the STUDENTS table, only some of the lines are displayed in QGIS.

Here is the SQL:

update
  public.centerprogram_mi_autism_2018 as students
set
  geom_line = data.geom_new
from
  (

    select
      SCH.schnum
      , a.studentnumber
      , ST_MakeLine(a.geom,
                    sch.geom) as geom_new
    from
      "public"."centerprogram_mi_autism_2018" as a
      join dpsdata."Schools_Current" as sch
        on
          a.schoolnum = sch.schnum
  ) as data
where
  data.schnum = students.schoolnum

And here are the results - when run the first time (dashed black lines):

And the weird part, if I run the same query again, the lines that display are different (dashed black lines):

Troubleshooting steps and results:

In QGIS, selecting the features by rectangle seems to select all the
features, including the ones that aren't being displayed
SRID of all inputs and outputs seem to be fine (4326)
School number values all return correctly in the sub-query
Coordinates of the input points are actually the same as the coordinates of the nodes of the result lines
Results are the same in QGIS 2.18 and QGIS 3.2
PostgreSQL 9.5 w/ PostGIS 2.3.7

I've used this plenty of times before but can't figure out what's going on.

Perhaps the UPDATE syntax isn't correct, but again, it seems only the display is incorrect?


texts: Add layer from a folder to Map using java script api

texts: I know you can add layers to map through map service using javascript api. Can we directly add layers (placed in a folder on ftp) to the map using javascrip api. If so, can please explain it and provide me a sample. Thanks.


texts: Zoom GeoTools JMapPane to retain data area on resize

texts:
By intention, when the JMapPane is resized, the contents of the pane do not change  [1].  So, as the application window gets narrower, content from the right-hand side is trimmed off.

An alternate GIS visualisation standard is to zoom the contents of the window to ensure the same mapping area is covered as the application size changes - see the video at [2].

Is there a simple setting that will force the second behaviour in GeoTool JMapPane ? Or, do I need to handle component resized, retrieve the previous data area, and force the map to display that area again?

[1]. http://docs.geotools.org/stable/userguide/unsupported/swing/jmappane.html#how-map-position-and-scale-are-calculated

[2]. https://youtu.be/lovWETSsXe4


texts: Calculating Solar Radiation falling on vertical surfaces using DEMs and ArcGIS for Desktop?

texts:
I am using ArcGIS10 Area Solar Radiation tool to calculate irradiation on horizontal surfaces (building roofs) but need to do the same for vertical surfaces (building facades).

Any ideas on how to do this using ArcGIS and DEMs?

I am open to broader suggestions on what tools to use, just needs to run for large portions of a city (neighbourhood scale) at once, instead of building level.


texts: Using GDAL in Command Line with QGIS 3.10.6 on Mac

texts:
I'm fairly new to GIS. I want to run GDAL in the command line on my Mac running MacOS Catalina. My understanding is that QGIS comes with GDAL. GDAL plugins seem to work in my installation of QGIS.
I have tried copying the console output from running GDAL plugins within QGIS to a terminal window and I get a 'command not found' error. Same deal with just typing the commands into terminal. I have successfully installed GDAL on its own via Homebrew and used it in the command line on another computer, so I am familiar with its usage.
Are there steps I've missed before I can use GDAL from the command line?


texts: Reusing ArcObjects code written at 9.3 in ArcGIS 10.2?

texts:
-- update:

i've narrowed down the scope of the question to something more specific in a new question:

argis 9.3 to arcgis 10.2 arcobjects connection AutomationException: 0x5 - Access is denied. in 'Invoke'

-- end of update.

I have legacy code, written in 9.3. makes extensive use of arcobjects, now the company bought 10.2, and after updating all urls, it turns out that all arcobject calls are no longer supported, and it seems we have to write a server object extension, which is discouraged even in esri's docs,  to wrap the old code.

Am I missing something here? are arcobjects calls really no longer supported in 10? is there a way to reuse the old code that I might be not aware of?


texts: Reclassifying certain values using ArcPy in ArcGIS Pro?

texts:
I need to reclassify certain values to different values in a raster using Python in ArcGIS Pro 2.2.2. I can't use reclassify because I don't want to change all values and I can't go through each one of them because there is a lot. When I use Raster calculator I get an attribute error.

Is there a way to do that?


texts: Computing an averaged latitude and longitude coordinates

texts:
How can I compute the average between several latitude and longitude spots?

Should I just compute the Arithmetic mean for both lat and lng?


texts:
To create applications that manipulate ArcObjects there is the VBA and Java runtimes.
To create application that process data with ArcToolbox tools, there is ArcPy for Python.

Today, I discovered that there is a C and Java APIs for ArcSDE! What are ArcSDE APIs for? Do they provide exactly the same functionalities as ArcSDE command lines?

Are there any other ArcGIS API?


texts: QGIS web client email link URL problems on windows

texts:
I've got QGIS server (2.6) and QGIS web client working on a windows machine. Eventually I'll get it running on linux, but for now I need to get it working on my localhost, because I need to prove this system can work for a bigger rollout on another host.

I'm trying to get the map button that allows you to email a link working. When I click the link, I get the URL:

http://localhost/web/site/qgiswebclient.html?map=/C:OSGEO4W32/apache/htdocs/web/projects/project.qgs&startExtent=18748077.826615%2C-5847068.7041556%2C18748900.718997%2C-5846618.442286&visibleBackgroundLayer=Google%20Satellite&visibleLayers=project&initialLayerOrder=Highways%2CCadastral9%2CStudy_Area_Boundary&selection=

the map gui loads, but the map just sits there "loading". I'm pretty sure it's not working because the URL is being created with this path:

?map=/C:OSGEO4W32/....

instead of

?map=C:/OSGEO4W32/....

I'm assuming this URL is created by the GetUrlParams.js in my web client files. That script is here:
GetUrlParams.js

Does anyone know  how to edit this script to correct the location of the forward slash from before the C: to after the C:?

Update: inspecting the part suggested below, I checked what urlParams.map is giving me with an alert:

urlParams = Ext.urlDecode(urlArray[1]);
    if (norewrite) {
        if (urlParams.map == undefined) {
            alert(errMessageStartupMapParamString[lang]);
            urlParamsOK = false;
        } else {
            alert(urlParams.map);
            wmsURI = serverAndCGI + "?map=" + urlParams.map + "&";
            printURI = printServer + "?map=" + urlParams.map + "&";
            wmsMapName = urlParams.map;
        }
    }

and that correctly gives me what should be there:

C:/OSGEO4W32/apache/htdocs/web/projects....

But that's not what's being written to the URL when I press the email link:

?map=/C:OSGEO4W32/....

So to me there's something else going on somewhere else when that email link is generated, but I don't understand where.
Does anyone know if the email linking feature is built into Openlayers, or is this some sort of plugin?
I've searched quite of few of the obvious places in the different javascript files for the Web Client and I can't find references to the email link - to try and narrow down where that URL is being generated.


texts: asynchronous execution javascript dojo (return a variable result of a query)

texts:
I would like to know if it's possible to return a variable in a function after finding the result in a queryTask.

Apparently this is an asynchronous execution but, to be honest, I really need to execute that ..

Here is my simple example I hope you find a solution to remedy it.

  myfunction: function (num) {

    var rft = 0;
    var gedlayer = new FeatureLayer("http:mylink/MapServer/1", {
        mode: FeatureLayer.MODE_ONDEMAND,
        outFields: ["*"]
    });

    var qr = new QueryTask();
    qr.where = "rft=" + num;
    qr.outFields = ["*"];

    gedlayer.queryFeatures(qr, function (featureSet) {
    if (featureSet.features.length > 0) {
    var content = "";
            for (var i = 0; i < featureSet.features.length; i++)
            {
            rft = rft + featureSet.features[i].attributes.rft;
                    }
    }
    },

    function (error) {
        console.log("Error dans la requette de poste source : " + error);
    }
    );

 return rft;
    }    ,


texts: Bulk convert DBF to CSV in a folder ArcGIS 10.1 using Python

texts:
I would like to convert 1000 dbf files to CSV. All dbf files are in one folder and they have sequential names (S:\output_tables\tcb1.dbf, S:\output_tables\tcb2.dbf, S:\output_tables\tcb3.dbf, up to S:\output_tables\tcb1000.dbf).

I would like S:\output_tables\1.csv, S:\output_tables\2.csv, S:\output_tables\3.csv up to S:\output_tables\1000.csv

I would like to eventually run in a ARCGIS 10.1 model so ArcPy would be the preferred solution but I'm open to other options.


texts: Failure to display points from a MySQL table in QGIS

texts:
I connected a MySQL table containing points with QGIS.
Adding new records from QGIS works well.
The problem is this: QGIS does not display the points on the map. However QGIS manages to interpret the coordinates because it zooms on the area containing the points.
The attribute table in QGIS looks correct.
I work with QGIS desktop 3.4.14.1 in French. I also tried in English I have the same problem.

Precision : when I export the layer in 'kml' the coordinates of the points are exported correctly.


texts:
I have the following info from a data provider and need to create a custom coordinate system using these parameters. The end goals is to map a file in this coordinate system against other data sources that are in a known coordinate system. Anyone have ideas?

*

  use the Geocentric Datum for Australia,1994 (GDA94)  as their mapping
  datum, with the  Lambert Conformal Conic (LCC) projection used in the
  CPR2 system System 1980 (GRS80) for which the  major semi-axis a = 6 378 137
  metres exactly and an  inverse flattening (l/f) of 298.257 222 101.

  The parameters for the LCC projection are: Standard parallels : 15
  degrees South , 39 degrees South. True origin of X,Y Cartesian
  coordinates : Intersection of Lat 27 degrees South and Long 134
  degrees East False Origin for E,N Cartesian coordinates : 5,000,000
  metres West and 5,000,000 metres South of the true origin.

  All coordinates given in LCC are E,N ie false origin and in centimetres in
  CPR.    When viewing the data in GDD, the coordinates will display in
  metres.LCC uses the reference ellipsoid  Geodetic ReferenceAll

*

I've tried to create a .prj file for use in ESRI software (AppData\Roaming\ESRI\Desktop10.5\ArcMap\Coordinate Systems\zzz.prj). Here is the result of that but defining a data doesn't correctly line up after defining a projection using this prj file:

INCORRECT:

  PROJCS["zzz",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",5000000.0],PARAMETER["False_Northing",5000000.0],PARAMETER["Central_Meridian",134.0],PARAMETER["Standard_Parallel_1",-15.0],PARAMETER["Standard_Parallel_2",39.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

CORRECT ANSWER (Standard_Parallel_2 needed to be : -39.0)

  PROJCS["GDA 1994 LCC CPR2
  T",GEOGCS["GCS_GRS_1980",DATUM["D_GRS_1980",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",5000000.0],PARAMETER["False_Northing",5000000.0],PARAMETER["Central_Meridian",134.0],PARAMETER["Standard_Parallel_1",-15.0],PARAMETER["Standard_Parallel_2",-39.0],PARAMETER["Latitude_Of_Origin",-27.0],UNIT["Meter",1.0]]


texts:
I would like to know the precision of a Geohash with a given length. If there is a 'simple' formula you can use to calculate it, that would be extra-cool.
Wikipedia lists the precision up to 8 characters:
#   km
1   ±2500
2   ±630
3   ±78
4   ±20
5   ±2.4
6   ±0.61
7   ±0.076
8   ±0.019


texts: Improve OGR SQL query for polygon intersection using WKT geometry

texts:
I want calculate the intersected area between a shapefile and WKT geometry.

My target big layer: a.shp
My query small layer: 'POLYGON((1150046 1264030, 1160100 1263793, 1160079 1260606, 1150218 1260822, 1150046 1264030))'

I have the following ogrinfo expression that is actually working:

ogrinfo a.shp -dialect sqlite  -sql "SELECT id, sum(ST_Area(ST_Intersection(
GEOMETRY, ST_GeomFromText('POLYGON((1150046 1264030, 1160100 1263793, 1160079 1260606, 1150218 1260822, 1150046 1264030))')))) as m2 FROM a WHERE ST_Intersects(GEOMETRY, ST_GeomFromText('POLYGON((1150046 1264030, 1160100 1263793, 1160079 1260606, 1150218 1260822, 1150046 1264030))'))"

The SQL query is:

SELECT col1, sum(ST_Area(ST_Intersection(
GEOMETRY, ST_GeomFromText('POLYGON(( ... ))')))) as col2
FROM a
WHERE ST_Intersects(GEOMETRY, ST_GeomFromText('POLYGON(( ... ))'))"

Could this SQL query be written better to get faster results?


texts: Merging feature classes in memory gives ERROR 001156?

texts:
arcpy.AddField_management("in_memory/originalfc","ALPHABET","TEXT",1)
arcpy.Select_analysis("in_memory/originalfc","in_memory/A","attr<=0")
arcpy.CalculateField_management("in_memory/A","ALPHABET",'"A"',"PYTHON")
arcpy.Select_analysis("in_memory/originalfc","in_memory/B","attr>0")
arcpy.CalculateField_management("in_memory/B","ALPHABET",'"B"',"PYTHON")
arcpy.Merge_management(["in_memory/A","in_memory/B"],"in_memory/merged")

When the merge executes I get an error 001156 on OID 1 that says i can't write value 'A' to output field ALPHA.

If I try to append one to the other I get the same error. They were originally the same feature class why won't they go back together?


texts: Versioning commercial datasets in ArcSDE + Oracle

texts:
We are building an application with several spatial tables, provided by a commercial vendor, at its core. The vendor provides updates quarterly. We maintain several tables of our own, derived from these commercial tables. Our application would benefit from versioning the commercial tables along with our derived tables, such that we can see a consistent view of things at any point in time.

Upon each update from the vendor, I would like to "diff" the changes and apply them to our database, then update our derived tables accordingly. Can this be accomplished with ESRI's implementation of versioning? Would this require registering our tables with the GDB? (I prefer to do the heavy lifting with PL/SQL, so I'm not sure how I feel about this.)


texts: How to view attribute information in shapefile using QGIS?

texts:
I have a polygon shapefile (*.shp) and want to view its attribute data.

How can I do this using QGIS?


texts: Model Builder prevent 'Add Field' In Multivalue Parameter

texts:
I made a Field Picker script to allow the user to select field names which are then used in another script for more processing with update and insert cursors.

Is it possible to prevent them from adding their own field with the Add Field button?

If this button cannot be removed, is there a better way to get user input where they can select Field names and I can pass those as strings to another script?


texts: Calculating parity of address range using to/from left/right ranges using ArcGIS for Desktop?

texts: I have been working on a road center line feature class with left/right, to/from address fields (in ArcGIS Desktop 10.2). Using these fields, is there a way to automate the process of deciphering if the "to" and "from" segments are even, odd, both, or neither(zero) for both "left" and "right" of each segment?


texts: r.topidx error: required parameter output not set

texts:
I'm calling GRASS r.topidx from QGIS 2.2.0. I can get other tools to work on the DTM, for example r.watershed. But r.topidx I get the Oops message:

Oooops! The following output layers could not be open
output: C:/Users/User/Documents/GIS DataBase/Gills/Terrain/TopIdx.tif
The above files could not be opened, which probably indicates that they were not correctly produced by the executed algorithm
Checking the log information might help you see why those layers were not created as expected
This algorithm requires GRASS to be run. A test to check if GRASS is correctly installed and configured in your system has been performed, with the following result:
GRASS seems to be correctly installed and configured

Looking at the log, I find:
C:\PROGRA~1\QGISVA~1\bin>r.topidx input=tmp14037940151265 Output topographic index layer=Output topographic index layer1a269dc6704c40399e5421e3bc50f16e --overwrite
Sorry  is not a valid option
Sorry  is not a valid option
Sorry  is not a valid option
Sorry,  is not a valid parameter
Sorry  is not a valid option
Sorry  is not a valid option
Sorry  is not a valid option
ERROR: Required parameter  not set:
    (Output topographic index map)

But, as you can see from the Ooops (very annoying) message, the output file is set to /TopIdx.tif. which should be OK?

Any idea what might be going wrong here?


texts: Feature Layers fail every second time the page loads, with "Error: Unable to complete operation"

texts:
If I create a very simple map based on the Feature Layer On Demand sample, and change the URL to my own feature layer, the map consistently fails every second time:

  Error: Unable to complete operation. {code: 400, message: "Unable to
  complete operation."
  message: "Unable to draw graphic (null): Unable to complete operation."

If I refresh the page, the layer loads correctly, with no error messages. Refresh it again and the error shows. This applies to multiple layers on the server, not just one, so I don't believe it's related to corrupt geometry or whatever.

I vaguely remember seeing this problem a few years ago, but I can't remember the resolution. The layer is in ArcGIS Server 10.11 and I'm running the simple map via IIS on Windows 7, and accessing it via http://localhost

The only change I've made to the sample is the layer's URL, the map's starting extent, and the infoWindow contents.


texts: ArcGIS-Server how to create a WMS with a name different from title?

texts:
Is it possible to create a WMS layer with ArcGIS Server that will have a name property différent of the title property. How can I do that with ArcGIS Server?

I want to have a layer in WMS like that, title different of the name property:

<Layer queryable="1" opaque="0" cascaded="0">
  <Name>lidar_mhc</Name>
  <Title>Hauteur de canopée LiDAR (dégradé)</Title>
  <KeywordList>
  <Keyword>produit dérivé</Keyword>
  <Keyword>lidar</Keyword>
  ...


texts: Calculate date field using date variable and time attribute in a feature class

texts:
I have a string variable (fileDateString) with the date in this format: 20150404

In my feature class there is a text field with time in twenty-four hour UTC format. There is also 2 leading spaces before the time. E.g. "  1300"

I have added a Date field to my feature class and I would like to calculate a date/time. I also want to be able to subtract 6 hours from the date/time (convert UTC to CST) for the final result.

I thought about first creating a text field with the date and time, then converting that to date but am wondering if there is an easier way?

expression = "\"" + str(fileDateString) + "\"" + "!DateUTC!"


texts: Getting rid of polygons and lines within polygon?

texts:
I wanted to create features for all the streets in my area and I did so by differentiating parcel data over my boundary. I got what I wanted but it looks like it also grabbed the small slivers between parcels. I do not want these. How do I go about cleaning this up? I can't delete the small polygons as it is all one single feature.

See images here


texts: Has anyone implemented the leaflet.shapefile plugin?

texts:
I'm trying to implement the leaflet.shapefile plugin to enable a user to upload a zipped shapefile and overlay it on my webmapping site as can be done in the demo version of the site. http://leaflet.calvinmetcalf.com/

When the page I've created loads before I even click the button to upload a file I get an error which is from the the setWorkerEvents() function (line124) in the script.js file. When I click to upload a file it goes through the steps but seems to stop when the worker.data function is called from Line 76 of the code script.js.

I created a gist https://gist.github.com/dodriscoll/4a96f825ebe8123070b4 with 3 of the files being used, I think the issue is with the script.js file.


texts: Create Archimedean spiral with geographic coordinates

texts:
I don't know if this is the right place to ask but here it his.

I want to create an Archimedean spiral based in some location. Lets say I have this set of geographic coordinates (48.860770, 2.337655). This will be the initial point of the spiral and I want to create 10 new locations.

I have this formula for calculating in x, y coordinate:

for (i = 0; i <= 10; i++) {
    angle = 0.1 * i;
    x = centerx + (a + b * angle) * Math.cos(angle);
    y = centery + (a + b * angle) * Math.sin(angle);
}

How can I change this to convert to geographic coordinates?

In my example, centerx would be 48.860770 and centery would be 2.337655


texts: QGIS -> Add VectorLayer DXF import fails: "invalid"

texts:
Trying to import a DXF file with Add VectorLayer is resulting in the error "not a valid or recognized data source."

Here is the DXF file that gives the error.  It is created by exporting from OpenSCAD.

However, I believe it's a valid file.  I am able import it in the Autodesk A360 Viewer. Here it is in the viewer.

I tried to use the DXF2SHP plugin, but subsequently read that it has been deprecated and is no longer maintained.

Any idea why this is failing?

Using QGIS 2.16.3


texts: I have created a Python script which compresses an ArcSDE geodatabase and now I would like to view the State ID after the compress.  Is there a way with Python I can get this info?


texts: My computer died. However, I had everything backed up on Dropbox. So when I got a new one, I was able to just download everything onto it, retaining the exact same file structure under C:. After installing QGIS (the same version, 3.4), everything looks good except that the Bookmarks panel is empty. My understanding is that QGIS now stores bookmarks in the project file. Either way, I'm puzzled about them having disappeared. A search of this site turned up nothing pertinent to this problem.


texts: Accurately find the bounding box of geography polygons

texts:
I'm displaying polygons on a map and have had issues with the spatial index where the index performs well only for certain data sets (large areas with low density / high density small map areas).

I've added bounding box columns (min/max, x/y) to the table with spatial data and using an index on these fields improves the speed of queries for the map by 100%.

However for large polygons the bounding box I've calculated is quite inaccurate - I'm using:

alter table contactgeography
add
CGMinLong AS (CONVERT(float, ContactPolygon.EnvelopeCenter().Long - ContactPolygon.EnvelopeAngle())) PERSISTED,
CGMinLat AS (CONVERT(float,ContactPolygon.EnvelopeCenter().Lat - ContactPolygon.EnvelopeAngle())) PERSISTED,
CGMaxLong AS (CONVERT(float,ContactPolygon.EnvelopeCenter().Long + ContactPolygon.EnvelopeAngle())) PERSISTED,
CGMaxLat AS (CONVERT(float,ContactPolygon.EnvelopeCenter().Lat + ContactPolygon.EnvelopeAngle())) PERSISTED

to add the columns. This will often return an area smaller than the original polygon (even the MS documentation says its not perfect).

I'd like to retrieve the points of the polygon in SQL and loop through them to find the min/max lat/lon but can't see a way to do this in a computed column.

Can the

ContactPolygon.EnvelopeCenter().Lat - ContactPolygon.EnvelopeAngle()

be improved to return more accurate results?

Is there an easy way to return the min/max points from the polygon itself?


texts:
Trying to create a list of well counts. In the code below, wellCntList only contains the current county count and not a list of ALL the counts. Sure I'm missing something simple.

for cname in countyList:
    whereclause = "{} = '{}'".format("COUNTY", cname)
    wellCnt = 0
    wellCntList = []
    with arcpy.da.SearchCursor(Wells_Intersect_Layer, "COUNTY", whereclause) as cursor:
        for row in cursor:
            wellCnt = wellCnt + 1
        wellCntList.append(wellCnt)
    print wellCntList


texts: Should be expecting any major issues with my custom add-ins after upgrading to 10.3 from 10.2.2?

texts:
Going from 9 to 10 was a major change as global VBA customizations had to be converted to add-ins requiring visual studio and the SDK.

There were some differences between minor "10" version as well that required some work; newer version of VS. .NEt framework, etc.

What should I expect after upgrading to version 10.3 from my current version 10.2.2.  I understand VS 2010 Express is no longer supported and I have to upgrade to 2012 when using Express? However, supported and recommended .Net framework version is still 3.5? Can someone share their experience and provide a quick rundown of what is required and what should I expect when upgrading/recompiling my add-ins.


texts: Following previous instructions to set up adding ecw files to QGIS does not work

texts:
I have followed the solutions to setting up ECW in QGIS.  However, none of them work.

In addition I get a mount - access denied error when trying to load the gdal-ecw and gdal-mrsid plug-ins using the OSGeo4W installer in advanced mode.  Any workarounds for this?


texts: Add field to layer in QGIS using Python in a plugin

texts:
I'm trying to write a plugin for QGIS that adds a field to a layer selected from a combobox and calculates area.

For that I import these modules:

from PyQt4.QtCore import QVariant, QSettings, QTranslator, qVersion, QCoreApplication
from PyQt4.QtGui import QAction, QIcon, QFileDialog

and in the run function:

def run(self):
    # here I populate a combo box in a form
    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layer_list.append(layer.name())
    self.dlg.Cmb_Capa.addItems(layer_list)

    self.dlg.show()
    result = self.dlg.exec_()
    if result:

Now here comes my problem... I find this in the documentation page of QGIS to add a field to a layer

if caps & QgsVectorDataProvider.AddAttributes:
    res = layer.dataProvider().addAttributes([QgsField("mytext", QVariant.String), QgsField("myint", QVariant.Int)])

but I don't know how to take the layer selected in the combo and apply this.

Can anyone give me a hand or a tip about how to do it?


texts: Create point on a line with given distance from end in PostGIS

texts: Does anyone know how to create points along a line (MultiLineString geometry) with a given distance in PostGIS? I'd like to create points at a regular distance over a line, which represents a shoreline.


texts:
I have a WFS layer, served by Geoserver, loaded in QGIS using the standard 'Add WFS layer' dialog. When I try to remove it using:

QgsMapLayerRegistry().instance().removeMapLayer(wfs_layer.id())

I get a full application crash, with no Python errors. This also happens when I use wfs_layer.reload() and wfs_layer.dataProvider().reloadData(). These methods previously worked on the same WFS layer. The service and content of the layer has not changed. The methods work when manipulating a (PostgreSQL) database layer.

I have made no changes to the code that I can see would affect the ability to reload data.

Any ideas or experiences?


texts: Check if listed fieldnames are in shapefiles

texts:
I have one main folder and it has multiple sub-folder. Each sub-folder has one shapefile. I want to check some attribute field names like PLACENAME, GNIS_ID, which are present in all the shapefiles or not. If not present those field name then print the shapefile name.

I have written small code and getting error. I am learning Python.

Code:

import arcpy, os, string
rootFolder = 'C:\Y4YK\Muni'
ctr = 1
for root, dirs, files in os.walk(rootFolder):
    for name in files:
        if name.endwith(".shp):
            feature_class = name
            fiedCheck = ["PLACENAME", "GNIS_ID"]
            fileds = arcpy.ListFields(feature_class)
            fieldName = [f.name for f in fields]
            for field in fieldCheck
                 if filed in filedName:
                    print "Field Exist"
                 else
                    print feature_class


texts:
In an answer to a question asked by @Canningmister, @tjmgis posted an excellent answer  Styling OS Vectormap Local in QGIS? about making OS SLD files work in QGIS.
As an amateur user who uses QGIS to produce maps for my Parish Council I have struggled with styles for OS data in QGIS.
I can see how to edit the sld files using a text editor such as TextPad which I assume will work with Mastermap area and line layers but I am confused with making the text layer SLD work and have several questions:

Are the layers imported into QGIS from ESRI shapefiles or directly from gml?
How do you post process the text layer data to add extra columns without having the data in a database and if so which one?

I am familiar with executing SQL scripts in a database. Am I missing something - I can't get my head round this.


texts:
I am trying to use gdal.FillNodata in python on a ".tif" raster. The nodata value is set to -3.402823e+38. Here is my code :

import gdal
from gdalconst import *

ET = gdal.Open("ET_sample_reclass.tif", GA_Update)
ETband = ET.GetRasterBand(1)

result = gdal.FillNodata(targetBand = ETband, maskBand = None,
                     maxSearchDist = 5, smoothingIterations = 0)

print result # return 0

This code seems to do actually nothing and I don't know how to define an output raster in python. I haven't found any example yet to enlighten me.


texts: Calculating LULC change 1973-2014 of a classified map (raster converted to shapefile) in ArcGIS 10.6

texts:
I am attempting to workout landcover change between two classified landcover maps. Ideally I want to show percentage change between the years per landcover type. I think I have done this but I am unsure whether my calculations are correct as the percentages seem to be very high...

Also, my main interest in the study is seeing how palm oil (large concessions and smallholders), has caused habitat loss (in this case loss of forest) from 1973 - 2014. Is there an elegant way of showing this in both map form and also graph/table form? This is the current landcover map I have made.

I haven't been able to work out a solution as of yet.


texts: How to Find the Average Raster Value of an area defined by a shapefile using R?

texts:
I have a set of raster images which represent a specific month over the years, and I want to make a timeline of the average values of an area using a shapefile.

How do I extract the values from the rasters and import them in R in a manner that I can use them?


texts: Automatically update field from weak entity on PostGIS/QGIS

texts:
I am having a DB centered around objects and a series of tables informing on different attributes for this object. On of these table/attribute is the diameter. It occurs the diameter is strongly linked with the material, i.e. there isn't any diameter having two different material. So the table with the diameter is linked with a weak entity (cardinality 1,1) containing the material datas.

In QGIS, I would like that when I select the diameter of an object I created (or that I am updating), it automatically fill the field relative to the material of the object. How can I do that ? Better even : I would like that the material field in the attribute table to not be editable, and so just be related to the diameter field.


texts: Raster Calculator in QGIS does not calculate the entire layer

texts:
I am working with two raster layers from WorldClim, using max temperature and minimum t. I want to obtain the range, meaning, doing tmax minus tmin. When I use raster calculator in order to do so, the result is a really small area.

The original raster

When I substract one raster layer with another (both of the same source, size, SRC, just different band values), this is the result:

How do I obtain the difference for the entire raster, and not just that small part?


texts:
I recently upgraded by QGIS from 2.18 to 3.4. Apart from a bunch of issues, I am not able to find the add feature tab. I want to create a set of polygons based on points.

Can someone help me with either letting me know what the issue is or a workaround?


texts: How to use SQL view request for WFS from GeoServer: Apply PostGIS function depending of zoom

texts:
I have a task
to allow a single WFS request to return two different responses depending on zoom:
return usual SELECT * FROM table if zoom level >n and
return SELECT ST_SimplifyPreserveTopology(geom, x) FROM table if zoom level <=n.

But I don't understand how to apply zoom level expression for parametric request to PostgreSQL. ...


texts: Any other FuseSearch.js examples in Leaflet?

texts:
Have been researching about using FuseSearch and using this example.
The only difference is that am using polygons.Its working well but not highlighting polygons on results click.
Has anyone done a similar app with polygons and has integrated the FuseSearch plugin?
Please help!!


texts: CartoDB Vis - Multiple layers mouse interaction?

texts:
My map visualization has two layers:

1) A city layer where the cities are represented as standard markers.

2) A county layer underlying the cities which has county level aggregate data.  The layer is a choropleth.

--

Can we get mouseover and / or mouseclick events working on both layers simultaneously?

I need to be able to click/hover a city or a county  and get the popup data response. (and probably highlight the county polygon)

At this point I can access one or the other depending on their order in the vis.  After looking around a bit it looks like at least at some point this was handled by a variety of workarounds.

I was hoping to get it going somewhat automagically.  If not, is this template the "right" way of doing it?  I could stuff all the data into one table but that feels like it will make more complication.  Users will be entering new marker points into the DB and it would be preferable to keep that data on its own.

https://github.com/mhkeller/cartodb-templates/blob/master/basic/point-polygon-hover.html


texts: Internal Server Error with Geocortex Essentials Manager?

texts:
I have installed Geocortex Essentials 4.1.0 and configured post-installations.
However, when I opened Geocortex Essentials Manager there occurred IIS 7.5 Detailed Error as follows:

  HTTP Error 500.19 - Internal Server Error  And Error Code 0x80070021:
  This configuration section cannot be used at this path. This happens
  when the section is locked at a parent level. Locking is either by
  default (overrideModeDefault="Deny"), or set explicitly by a location
  tag with overrideMode="Deny" or the legacy allowOverride="false".

I have attempted to fix with applicationHost.config file but confusing where to change.


texts: How to convert from lat-lon to UTM with gdaltransform?

texts:
I would like to transform the following coordinates

31.96212 -103.004715

to UTM with gdaltransform. I put coordinates into file coordinates.txt and ran

gdaltransform -s_srs epsg:4326 -t_srs epsg:25832 < coordinates.txt

and got

ERROR 1: latitude or longitude exceeded limits
transformation failed.

I also tried to reverse order to

-103.004715 31.96212

with the same result.

What I am doing wrong?


texts: Is it possible to rotate the map in layout manager?

texts: QGIS 3.12 - Often I have property boundaries that would better fit on an A3 output if rotated 45-90 degrees from that represented in the map view. I can rotate in the map view but the layout manager does not respect this rotation. Is there anyway to rotate the map in the layout manager?


texts: Splitting polygon at points using ArcPy with ArcMap

texts:
I am trying to divide irregularly shaped polygons at specific points. I don't have any fields that inform the angle at which the polygons should be divided.
I am using Create Thiessen Polygons which generally wants to build a polygon around a point, such that it is centered, rather than with the point at the edge.
Is there a way to force the polygons to be shifted such that they are along the points? Is there a different tool I should use? I have sketches below of input, current output, and desired output. (And the next step is to intersect the base polygon with the output Thiessen polygons; that is the final answer, but I have left it out here because it is a simple step.)
This is related to ArcGIS force Thiessen polygons to be larger or alternative tool?, but this is a new/separate question because of the shifting vs. centering issue.
import arcpy

#mypoly = filepath to my polygon SHP
#mypts = filepath to my points SHP
#myoutput = filepath to save Thiessen polgons output

arcpy.env.extent = mypoly
arcpy.CreateThiessenPolygons_analysis(mypts, myoutput)


texts:
I have a table in MG94 Z 54 (EPSG:28354).

My MapBasic statement is

Commit Table DCDB As "C:\DCDB_Full_Copy.TAB" TYPE NATIVE Charset "WindowsLatin1" CoordSys Earth Projection 8, 116, "m", 141, 0, 0.9996, 500000, 10000000 Bounds (-7745844.29605, 2035.05684549) (8745844.29605, 19997964.9432) Interactive
Update DCDB_Full_Copy Set X = CentroidX(obj)
Update DCDB_Full_Copy Set Y = CentroidY(obj)

It gives me

  22N14849 [LotPlan], 141.071716, -17.66865

but it should be

  22N14849 [LotPlan], 507606.510877, 8046471.56601

How can I get the coordinates in the correct format?


texts: Is there a PostGIS equivalent to ESRI Spatial Analyst COMBINE?

texts:
The ESRI Spatial Analyst COMBINE tool accepts multiple raster inputs and outputs a new raster with a unique number corresponding to unique combinations of values from the inputs.  It also provides a translation VAT.

I want to do this open source.  Does PostGIS have a way to do this?  GRASS?  Other options?  I am building my geospatial app with MapWindow and VB.net to run as a desktop Windows app.


texts:
First of all, the code which is below gives me the attribute names of a shapefile correctly if I use PyScripter IDE:

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

QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis-dev", True)
QgsApplication.initQgis()

layer = QgsVectorLayer("C:/istanbul/istanbul_geo-wgs84.shp", "sero", "ogr")

fields = layer.pendingFields()

for field in fields.itervalues():
    print "Name:%s Type:%s, Length:%s" % ( field.name(), field.typeName(), field.length())

If I want to access attribute names of a shapefile, Aptana Studio gives me an error message below:

    for field in fields.itervalues():
AttributeError: 'QgsFields' object has no attribute 'itervalues'
d:\src\qgis\src\core\qgsvectorlayer.cpp(188) : (QgsVectorLayer::~QgsVectorLayer) entered.

According to the console in Aptana Studio, my system can read vector layers:

d:\src\qgis\src\core\qgsvectorlayer.cpp(1628) : (QgsVectorLayer::setDataProvider)      Extent of layer: 28.0019662445739925,40.7878698688807120 : 29.9327258493754194,41.5270420672795026

d:\src\qgis\src\core\qgscoordinatereferencesystem.cpp(449) : (QgsCoordinateReferenceSystem::createFromWkt) authid recognized as EPSG:4326

I used OSGeo4W installer and configured my IDE based on Nathan's expressions:

How to add QGIS Python API Intellisense support in Aptana Studio3/Pydev

and finally my "aptana.bat" file to initialize the IDE:

@ECHO OFF
set OSGEO4W_ROOT=C:\OSGeo4W
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
call "%OSGEO4W_ROOT%"\apps\grass\grass-6.4.3\etc\env.bat
set PATH=%OSGEO4W_ROOT%\bin;%OSGEO4W_ROOT%\apps\qgis-dev\bin;%PATH%
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis-dev\python
set QGISHOME=%OSGEO4W_ROOT%\apps\qgis-dev
call "C:\Users\Geomatics\AppData\Local\Aptana Studio 3\studio3.bat" %*

How are we gonna tackle this problem? Thanks in advance...


texts: Get DD Lat and Lon from CRS EPSG:3111 - GDA94

texts:
I have a shapefile in CRS EPSG:3111 - GDA94  that I am trying to generate Lat and Lon, however when I use the field calculator and try to populate the lat lon I am getting values that I can't use as per the screenshot below -

I am expecting something in this format -

I have projected it in EPSG:4326 but didnt help


texts: Creating Raster ASCII file with `nlayar` and `layer_span` in ArcMap

texts: Is it possible to create this ASCII file with nlayar and layer_span in ArcMap 10.2?


texts: Continued issue with ECCODES_DEFINITION_PATH for use with cfgrib

texts:
I have found a few posts that address this general topic:

ECCODES ERROR when opening .GRIB file in Spyder/Visual Studio Code

https://github.com/ecmwf/cfgrib/issues/85

The first link is the exact error I have had.  I follow the recourse of using this command from the windows prompt:
setx ECCODES_DEFINITION_PATH "C:/Users/Scott_TSYL/anaconda3/envs/for_grib/Library/share/eccodes/definitions"

However, the same issue occurs except the new path was inserted before the old(incorrect) path, and they are separated by a semicolon (see image).
What is going wrong here?


texts: Extracting different number of points for different class in GEE

texts:
I want to write a function that takes class id and number of points for each class as argument and samples points from polygons, returning a FeatureCollection.
I could only write a function that samples the same number of points [100 points in the example] for each class with iteration as below:
GEE Link: Reproducible GEE Example
Let me know if you can't open the polygon file
// Reproducible Example for Training Point Sampling

// Create id List
var idList = ee.List.sequence(11, 13).map(function(x) {return ee.Number(x)});

// Extraction Function
var collectTrain = function(id, featureCol) {
  // cast
  id = ee.Number(id);
  featureCol = ee.FeatureCollection(featureCol);
  // Extract
  var point = ee.FeatureCollection.randomPoints(
   polygons.filter(ee.Filter.eq('id', id)), 100)
   .map(function(feat){
      feat = ee.Feature(feat);
      return feat.set('id', id);
   });

  return featureCol.merge(point);
};

// Initial FeatureCollection
var point = ee.FeatureCollection.randomPoints(
   polygons.filter(ee.Filter.eq('id', 11)), 1)
   .map(function(feat){
      feat = ee.Feature(feat);
      return feat.set('id', 11);
   });

//Iterate over idList starting from the initial point
var training = ee.FeatureCollection(idList.iterate(collectTrain, point));
print(training);

How could we write such a function that takes two lists as arguments, such as in R with purrr::map2 and works with iterating those two lists. Arguments should be like:
id = [11, 12, 13]
point = [100, 150, 120]


texts: Using Elseif Conditional Statement in QGIS Field Calculator?

texts:
I am trying to write an ELSEIF conditional statement in QGIS Field Calculator (version 1.8.0). I have used an example I found online:

CASE WHEN val < 0 THEN 'negative'
  WHEN val = 0 THEN "neutral'
  ELSE 'positive'
END

I modified the statement as follows:

CASE WHEN  "GRID_ID"  = 1 THEN 'complete'
  ELSEIF  "GRID_ID"  = 2 THEN "in progress'
  ELSE 'not started'
END

This statement would not run, the Output preview stated Expression is invalid.  The more info stated: Parser Error: syntax error, unexpected COLUMN_REF, expecting WHEN or ELSE or END

If anyone has had this error, what did you do to fix it?


texts: Updating ValueTable parameter of ArcPy Python Toolbox tool

texts:
I have a shapefile and I want to calculate multiple fields using the GP ValueTable parameter.
The script loads field names from a layer, then I should delete rows I don´t want to calculate and then I should write values to the second column to be filled in selected fields. Then it should arcpy.CalculateField.
But I am stuck at parameter settings. The script lies in the Python toolbox.
def getParameterInfo(self):
    """Define parameter definitions"""
    param0 = parameter("in_layer", "Vrstva", "GPFeatureLayer")
    param1 = parameter("field_values", "Atributova pole a hodnoty k vyplneni", "GPValueTable")

    param1.parameterDependencies = [param0.name]
    param1.columns = [["String", "Pole"], ["String", "Hodnota"]]

    params = [param0, param1]
    return params

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    if parameters[0].value:
        v_list = []
        fld_names = [f.name for f in arcpy.ListFields(parameters[0].value)]

        for fld in fld_names:
            v_list.append([fld, ""])
            del fld

        parameters[1].value = v_list

        del fld_names
        del v_list

    return

When I select a layer in the tool, it populates field names (great!) and it looks like this:

But...

If I want to delete some rows using the button on the right side, it does nothing, it does not delete the row. It just seems that it somehow "refreshes" the value table.
If I want to write a value to the second column, it lets me but after hitting enter or clicking anywhere else it erases the written value and sets nothing, empty value again. It does not hold what I am writing there.

How do I make the ValueTable work correctly?


texts: Layers not displayed in QGIS print composer legend

texts:
When creating a print composer legend some of my layers (all point shape files) are not shown in the legend on the map.

I realized that the symbol markers (3 dots) were different to the ones that were displayed (one dot) but since I am very new to GIS I have no idea what to take from this and it would be great if you could help me out.

(I tried reimporting my shape files, deleted the item, added it again but that had no effect.)


texts: How to directly copy value from attribute table associated with a selected feature?

texts: I need to copy the attribute value (only one particular attribute) of a selected point and save it in a text file. Is there any Python script or plugin which can do this for me? I need to save a lot of points so individually copying the values is going to be tedious.


texts: find the population density of some points in qgis/postgis

texts:
I have the tiger census data of population density in each block group. I also have a different kind of data: traffic accident data, where an accident has a lat/lon pair. My goal is to find the population density of each accident location.

How can I match the accident to the block groups?


texts: How to calculate pixels per meter ratio, according to Google or Bing map zoom level?

texts:
I have an image and detect two points in the image. I want to calculate the distance between these two points. I'm using Google map and Bing map.

I have

x,y pixel values for both two points
zoom level 19, or something else

What I need is the pixels per meter ratio.


texts: Polygon centroid coordinates update error

texts:
After a week of unsuccessful fighting I'm forced to ask for your help.

I have created a table 't_excavations' with fields 'id_excavation', 'geom' (MULTIPOLYGONE), 'longitude', 'latitude'

CREATE TABLE public.t_excavations ( id_excavation serial NOT NULL PRIMARY KEY, longitude numeric(9,6), latitude numeric(9,6), geom public.geometry ) WITH ( OIDS = FALSE );

Then I created a rule to auto update 'longitude' and 'latitude' columns to get XY of the polygon centroid that I add through QGIS

CREATE OR REPLACE RULE "insert_longlat" AS ON INSERT TO "public"."t_excavations" DO (UPDATE t_excavations SET longitude = st_x(st_centroid(t_excavations.geom)); UPDATE t_excavations SET latitude = st_y(st_centroid(t_excavations.geom)));

This rule works when I save changes in QGIS and coordinates are calculated fine

Now I want to update 'longitude' and 'latitude' columns if 'geom' field is changed (for example if I move the entire polygon or just one vertice)

I created a function:

CREATE OR REPLACE FUNCTION update_longlat() RETURNS trigger AS $BODY$ BEGIN NEW.longitude = ST_X(ST_Centroid(geom)); NEW.latitude = ST_Y(ST_Centroid(geom)); RETURN NEW; END; $BODY$ LANGUAGE plpgsql VOLATILE COST 100; ALTER FUNCTION update_longitude() OWNER TO postgres;

and a trigger:

 CREATE TRIGGER tr_update_longlat BEFORE UPDATE ON t_excavations FOR EACH ROW EXECUTE PROCEDURE update_longlat();

And now when I try to add new or edit existing polygon in QGIS I get this error message from postgres:

  Can't change layer t_excavations

  Errors: ERROR: not added 1 object.

Errors of the source:

  Error PostGIS while adding objects: ERROR: column "geom" doesn't exist
  LINE 1: SELECT ST_X(ST_Centroid(geom)) ^ QUERY: SELECT
  ST_X(ST_Centroid(geom)) CONTEXT: function PL/pgSQL update_longlat(),
  row 3, assignment operator

What I do wrong with the function/trigger?


texts: Clip raster by polygon, reassign values of clipped raster, merged reassigned raster with original

texts:
I'm trying to clip out part of a raster (DEM), change the values within the clipped raster and then merge/join the edited raster back into the original. The context is that I have a high resolution DEM (50 cm) representing a hill and want to show the final shape of the hill after proposed earthworks have taken place (the removal of three metres of soil in a given area).

I've followed this process:

Extract by Mask (polygon representing proposed earthworks location)
Create edited raster of earthworks area by reassigning values of clipped raster
('clipped_raster - 3')
Used raster calculator as follows:

The output I'm getting is limited to the clipped raster, I'm losing the extent of the original raster.


texts:
Processing Toolbox missing after I copied a python script from the Scripts folder to my local folder, and QGIS crashes when I try to active the "Processing" under the Plugin Manager

Re-installation also doesn't work! QGIS crashes when I try to activate it using Plugins Manager

It looks like, no one has encountered this problem!!!

There is something getting saved somewhere, or some settings getting saved somewhere, these are stopping the Processing Toolbox getting activated even after repeated installations. Cleaning of all the folders also doesn't help. Installing of stable versions also doesn't solve this problem.

If some source code is available where the Processing plugin within the Plugin Manager can be changed, perhaps might solve the problem. Developers may please help!

Otherwise, the only solution, to my mind,  appears to be to reload the OS (in my case, Windows 64-bit) and install the QGIS application if one really wants to work with the Processing Toolbox!

Interestingly, when I added another login account in this PC, the Processing Toolbox is coming in the QGIS in this new login account! However, I did not install QGIS in this login account separately. This shows some settings got permanently changed in the previous login account where the Processing Toolbox still doesn't appear!


texts: Loading CSV as a GIS layer into QGIS 3.4 not possible with certain coordinate reference systems

texts:
I tried multiple times to load a CSV (CRS: EPSG 4647) into QGIS 3.4.
While it works without problems in QGIS 2.x versions, there was no success loading it into QGIS 3.4 ('QGIS Madeira').

Do you have any suggestions of what I should pay attention to?


texts: Calculating line segment lengths within polygon using ArcGIS Desktop?

texts: I have a layer with line segments and a layer of polygons (in this case, buffers around points). I need to create a new field for each polygon that contains the total length of all line segments contained within that polygon, including portions of the line segments with the polygon that straddle the borders.


texts: Looking for Proxy Server (like MapProxy) but for Vector/Feature Data

texts:
Currently, I'm using MapProxy to pull raster tiles from web services, cache those tiles locally, then serve them back out.  So MapProxy is acting as a performance enhancing proxy.

Is there any piece of software that can function similar to MapProxy, but for feature data, like Vector Tiles or WFS?


texts: ArcGIS Pro and how to use Python to calculate an output from multiple table columns

texts:
I'm not sure how to work with tables in ArcGIS Pro using Python. I've done some coding in other areas but whenever I see the code snippets related to table operations my mind just goes cross-eyed. I'd like to rectify that.

I have a table with 5 columns; [GP_Min, GP_Max, Z_Min, Z_Max, ZUp_or_Down]
These relate to the dwelling units per acre for a cities general plan in a min-max setup such that for a general plan such as 'Agriculture' the GP_Min = 0 and GP_Max = 0.5. The conforming zoning would also be 'Agriculture' and has a Z_Min = 0 and Z_Max = 0.5.

However, not all 'Agriculture' general plan designations have a zoning of 'agriculture'.

What I would like to do, in layman's terms, is something like this:

If the general plan is not equal to the zoning then, if the general plan dwelling units per acre is greater than the zoning dwelling units per acre, define this as Up-Zoning, if it is less than the zoning dwelling units per acre, define this as Down-Zoning.

Any resources or guidance?