texts: How can I Setup Python to recognize QGIS

texts:
I want to use QGIS and work with Python out of the console, but I'm having troubles letting the system recognize the paths of some DLLs.
Since I couldn't find a good answer on the web, I wonder if I did something wrong.
I'm working with Windows 7 64 bit version, using Python 2.7 and QGIS 1.8.
In the system environment, I've created and pointed the variable PYTHONPATH to the python directory (which in my case is "C:\Quantum_GIS_Lisboa\apps\qgis\python"). Then I've modified the "Path" variable adding the address where "qgis_core" and "qgis_gui" are found (in my case "C:\Quantum_GIS_Lisboa\apps\qgis\bin"). The instructions I've followed are those of the "PyQGIS Developer Cookbook" ("http://www.qgis.org/pyqgis-cookbook").
It seems the problem relies in recognizing the DLLs, which are in the folder I've pointed to for the variable "Path".
Here is the error I have in the Shell when I try to import the module "qgis.core":

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
     import qgis.core
ImportError: DLL load failed: Impossibile trovare il modulo specificato.

Moreover, if I import the sys module and print the path to see what's there, I'm not finding what I expect:

>>> import sys
>>> print sys.path
['C:\\ProgramData\\Microsoft\\Windows\\Start Menu\\Programs\\ArcGIS\\Python 2.7', 'C:\\Python27\\ArcGIS10.1\\Lib\\idlelib', 'C:\\Quantum_GIS_Lisboa\\apps\\qgis\\python', 'C:\\Windows\\system32\\python27.zip', 'C:\\Python27\\ArcGIS10.1\\DLLs', 'C:\\Python27\\ArcGIS10.1\\lib', 'C:\\Python27\\ArcGIS10.1\\lib\\plat-win', 'C:\\Python27\\ArcGIS10.1\\lib\\lib-tk', 'C:\\Python27\\ArcGIS10.1', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\bin', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\arcpy', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\ArcToolbox\\Scripts', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages\\win32', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages\\win32\\lib', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages\\Pythonwin']

Hope someone could help me solving this, I'd be really thankful!

I've just tried to copy and paste the whole content of Python27 folder within the main "Quantum GIS Lisboa" folder in the Python27 main folder in my hard drive C, since I saw ArcGIS and Python were working fine like that. Then I tried to launch the IDLE found in "lib/idlelib" ("idle.pyw") from the new folder but this didn't solve the problem.
Everything seems to point to the right folders in the system environment.
Here is the print of the sys.path after the "reverse" trick:

['C:\\Quantum_GIS_Lisboa\\apps\\qgis\\bin', 'C:\\Quantum_GIS_Lisboa\\apps\\qgis\\python', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages\\Pythonwin', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages\\win32\\lib',
'C:\\Python27\\ArcGIS10.1\\lib\\site-packages\\win32', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\ArcToolbox\\Scripts', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\arcpy', 'C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\bin', 'C:\\Python27\\ArcGIS10.1\\lib\\site-packages', 'C:\\Python27\\ArcGIS10.1',
'C:\\Python27\\ArcGIS10.1\\lib\\lib-tk', 'C:\\Python27\\ArcGIS10.1\\lib\\plat-win', 'C:\\Python27\\ArcGIS10.1\\lib', 'C:\\Python27\\ArcGIS10.1\\DLLs', 'C:\\Windows\\system32\\python27.zip', 'C:\\Quantum_GIS_Lisboa\\apps\\qgis\\bin', 'C:\\Quantum_GIS_Lisboa\\apps\\qgis\\python', 'C:\\Python27\\ArcGIS10.1\\Lib\\idlelib', 'C:\\Users\\Umberto\\Desktop']

Traceback (most recent call last):
  File "C:\Users\Umberto\Desktop\set_qgis_01.py", line 5, in <module>
    import qgis.core
ImportError: DLL load failed: Impossibile trovare il modulo specificato.

Am I suppose to point to the Quantum GIS main folder as well as it contains a folder named "lib"?


texts:
I have script for exporting separate data driven pages (script tool in a toolbox).

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
pageNumbers = [4,7,11] #select pages before exporting
    for pageNum in pageNumbers:
    mxd.dataDrivenPages.currentPageID = pageNum
    arcpy.mapping.ExportToJPEG(mxd, r"C:\Path" + str(pageNum) + ".jpg")
del mxd

Is it possible to change the page numbers in the brackets without changing the code manually?

I mean, after clicking the link, you will be brought to a window to choose the pages you would like to export. There may be checkboxes drawn next to the page numbers or text boxes with the page numbers.


texts:
I want to perform Image Classification(JPEG Image).
I loaded the image in QGIS.
Everytime i perform an action, i get a screen pop-up of Co-ordinate Reference System(CRS).
How should that be handled?

Secondly, i bypass CRS by selecting WGS84 and later when i am about to add the ROI , i receive the following error
**

*- "AttributeError: 'NoneType' object has no attribute 'RasterCount'


texts: CARTO.js changing the styling only on one segment of a line with multiple segments

texts:
I'm attempting to change the styling of a line segment when clicked on a layer that has 4 segments. Only one segment's style would be changed at a time. Each segment has it's own column within a CARTO dataset. I was wondering if it was possible to accomplish this using only one dataset and not creating separate datasets for each segment? As of now when a line segment is clicked, all segments get the new styling applied.

Example of what I am trying to accomplish:

Below is the code I am using to change the styling:

routeLayer.on('featureClicked', featureEvent => {
  routeStyle.setContent(`
    #layer {
        line-width: 6;
        line-color: #c1d82f;
        line-opacity: 1;
     }
  `);
});

JSFiddle example.


texts: 'object' does not contain a definition for 'AddValue'

texts:
I have been able to get past this error, but I don't know why it was happening.

I'm attempting to set symbology on a layer using IUniqueValueRenderer in ArcMap 10.4. When I try and use the AddValue method I receive the 'object does not contain a definition for 'AddValue'.

This error is encountered at runtime.

var uvRenderer = new ESRI.ArcGIS.Carto.UniqueValueRendererClass() as IUniqueValueRenderer;
uvRenderer.FieldCount = 1;
uvRenderer.set_Field(0, "SampleField");
uvRenderer.DefaultSymbol = graySymbol as ISymbol; //a symbol for a gray bordered polygon with no fill
uvRenderer.UseDefaultSymbol = true;

int idIndex = SampleAreasFC.Fields.FindField("SampleField");

if (idIndex != -1)
    //the error occurs on this line. SampleAreaFeature.get_Value(idIndex) returns a non-null value.
    uvRenderer.AddValue(SampleAreaFeature.get_Value(idIndex).ToString(), "SampleField", redSymbol as ISymbol);

The way to fix this error was to do this instead:

 string value = SampleAreaFeature.get_Value(idIndex).ToString();
 uvRenderer.AddValue(value, "SampleField", redSymbol as ISymbol);

Even more strange was while that worked, the following generated the same error as before:

 var value = SampleAreaFeature.get_Value(idIndex).ToString();
 uvRenderer.AddValue(value, "SampleField", redSymbol as ISymbol);

Since var is just an implicit type, it should be functionally identical to the former, but still generated the same error as before.

Does anyone know what would cause this behavior?

UPDATE:

I also asked this question on StackOverflow because I thought there was a chance this may be a more general .NET issue rather than specifically ArcObjects. That question is here. However, that is not resolved either.

This isn't really the ArcObjects question I initially thought, but it may still be useful to leave up.

The core issue is that get_value returns a dynamic type and that shadows the ToString method making my 'var' declaration on 'value' be a dynamic type. The issue is then that the .NET RuntimeBinder isn't able to find AddValue at runtime.

In this case get_value returns an integer and ToString correctly gives me the string "347". AddValue can be found via reflection at runtime, so it seems like there may be an issue with finding the function on the associated COM object when passing a dynamically typed variable.

So the more specific question is now, why can't RuntimeBinder find AddValue when a valid string is passed in as a dynamic variable?


texts: Installed pywin 32 -- still getting error -- "No module named win32com.client."

texts:
I'm trying to install pywin32 (32 bit 2.7) so that I can run python from ArcGIS 10.2.1 (32 bit), and have access to the pywin modules. But I am getting the error  "No module named win32com.client." I installed from here.

It created a whole separate installation folder, so if I run pywin from there, and try to import arcpy, it can't find the module. And if I go into ArcGIS and try to import win32com.client, it can't find the module. So I'm not sure how to install it so that I can get all the modules loaded into the LIB of my ArcGIS installation.

I have no admin priviledges at all, I have to ask an admin, so I can't even see the env variables. But I know in ArcGIS, it's running python from that ArcGIS 10.2 folder. Was this installed wrong? Or do you have any ideas on how I can fix this? I have to edit ArcFM features, and apparently you can't do it via python without pywin.


texts: Is there an open source international GeoCoding database if you only need city level resolution?

texts: I'm writing an application that needs to plot a user's location on the map, and the user could be anywhere in the world, but I only need city/town level resolution, not street address.  Is there a single open source database that contains this information?


texts: QT Designer form - Combobox stores value but doesn't display it

texts:
I have a custom form created with QT Designer on which I have several comboboxes among other stuff. All the fields have the correct names corresponding to the database.

One of these comboboxes has a strange behavior. When I choose a value from the dropdown it is properly stored in the database. I can see it in the attributes table. But when I'm trying to edit the values on the form, every other value (either on other comboboxes, text fields or checkboxes) displays correctly, except the one combobox that always displays the first value in the dropdown menu despite the fact that some other value may be stored in the database.

Any thoughts on what's wrong?


texts:
I'm currently developing an ArcGIS based flex web viewer. ArcGIS is great for most of the things but there are some issues like performance and license costs.

That's why I'm looking into possible alternatives. At the moment I've found 3 possible alternatives, but I'm not sure what they are capable of, so I'm wondering if some people on here have more experience with them.

I have been reading the docs and will try to experiment with some of these, but it will take a lot of time to find out all these requirements. So it would be of great help if you could point me in the direction of software I should try out.

These are our requirements:

fast
easy scalability
advanced labeling (different label colors based on attributes, labels placed in a rectangle)
render cached data
network tracing possibilities (for sewers, so trace upstream and downstream)
printing (generating pdf's)

nice to have:

good flex or javascript viewer (I know this is different software, but I'm not sure if all 3 packages support those equally well)
easy to start with


texts:
How to configure a available feature types in GeoServer using python gsconfig?
I have created work space and data store using gsconfig. Now , how i can configure available layers in a exist data store (data store other than postgis table, data stores like directory of shape files)?

I have tried: as follow

from geoserver.layer import Layer
from geoserver.resource import FeatureType, Coverage
from geoserver.store import coveragestore_from_index, datastore_from_index, \
    wmsstore_from_index, UnsavedDataStore, \
    UnsavedCoverageStore, UnsavedWmsStore
from geoserver.style import Style
from geoserver.support import prepare_upload_bundle, url, _decode_list, _decode_dict
from geoserver.layergroup import LayerGroup, UnsavedLayerGroup
from geoserver.workspace import workspace_from_index, Workspace

def publish_featuretype(self, name, store, native_crs, srs=None):
        '''Publish a featuretype from data in an existing store'''
        # @todo native_srs doesn't seem to get detected, even when in the DB
        # metadata (at least for postgis in geometry_columns) and then there
        # will be a misconfigured layer
        if native_crs is None: raise ValueError("must specify native_crs")
        srs = srs or native_crs
        feature_type = FeatureType(self, store.workspace, store, name)
        # because name is the in FeatureType base class, work around that
        # and hack in these others that don't have xml properties
        feature_type.dirty['name'] = name
        feature_type.dirty['srs'] = srs
        feature_type.dirty['nativeCRS'] = native_crs
        feature_type.enabled = True
        feature_type.title = name
        headers = {
            "Content-type": "application/xml",
            "Accept": "application/xml"
        }
        headers, response = self.http.request(store.resource_url, "POST", feature_type.message(), headers)
        feature_type.fetch()
        return feature_type
publish_featuretype("self" name="Layername", store="store_name", native_crs="4326", srs="EPSG:4326")

I got the error as follow

AttributeError: 'str' objetcs has no attribute 'workspace'

code has taken from : http://pydoc.net/Python/GeobricksGeoserverManager/0.1.0/geobricks_geoserver_manager.core.geoserver_manager_core/


texts:
I've written a python function that calls a tool (Import CAD to GDB) in a loop. When I run it I have seen that it hangs sometimes, and I have to force the script to quit. I am wondering if there is a way to kill the tool execution after a certain amount of time, and move on to the next item in the loop? Here is what I have so far - I've included a timer but it still hangs:

import os, arcpy, arcgisscripting, time, sys
gp = arcgisscripting.create()

def ConvertCADtoGDB(srch_dir,gdb_dir,timeout):
    for dirpath, dirnames, files in os.walk(srch_dir):
        for file in files:
            f_name, f_ext = os.path.splitext(file)
            in_dgn = dirpath+'\\'+file
            out_gdb = gdb_dir+'\\'+f_name+'.gdb'
            if f_ext == '.dgn':
                if not os.path.isdir(out_gdb):
                    start = time.time()
                    end = start + timeout
                    while time.time() < end:
                        print 'creating',out_gdb
                        gp.ImportCAD_conversion(in_dgn,out_gdb,'','Explode_Complex')


texts: Extract raster cell values by polygon field

texts:
I have a raster as a RasterLayer (raster) and a polygon layer (bound) with 38 individual polygons as a SpatialPolygonsDataFrame. I used the polygon to mask from a larger raster. Each polygon contains either 4 or 5 individual pixels.

I want to extract the cell values by the name of each individual polygon. I know I can do this:

v <- extract(raster, bound)

v.counts <- lapply(v,table)

But this ends up looking like this:

> v.counts

[[1]]

0.0463505908846855 0.0503523759543896 0.0531094446778297 0.0665162950754166
             1                  1                  1                  1

[[2]]

0.0755022615194321 0.0766473188996315 0.0860441327095032  0.612392961978912
0.689212620258331
             1                  1                  1                  1                  1

Is there a different way to do this or a way to get the name field associated with the bound polygon to appear in the parentheses in the above output? For example [[Reservoir A]] instead of [1]?

Could this have anything to do with how the individual polygons are named in the shapefile? For example if you view the "polygons" portion in view, each polygon is named [1], [[2]] etc. :

Is there a way to just simply change the name of the polygon to get my intended output?


texts: Is it possible to open WMC(Web map Context) file in QGIS?

texts:
I am new to WMC and in our application we had a requirement to extend the WMC functions.

For that I need to know whether QGIS have the support of WMC files.

Or we need any plugins for WMC support in the QGIS.

Is WMC (wmc.xml) files able to open in QGIS as like in the below link in OpenLayers
OpenLayers WMC example


texts: How to set default style of layer using rest api in Geoserver

texts:
I want to change the default style of the layer in Geo server using the REST API.

For this i have this code:

  curl -v -u admin:geoserver -XPUT -H "Content-type: text/xml"
  -d <layer><defaultStyle><name>polygon</name></defaultStyle></layer>"
  http://mindcrewgis.com/geoserver/rest/layers/pymble:pymble_water

when  run this command  it is working fine form me.
But when i convert this command in PHP then it is not working for me
My PHP code for above curl command is

<?php

    $url="http://mindcrewgis.com/geoserver/rest/layers/pymble:pymble_water";
    $contentType = 'text/xml';
    $ch = curl_init();
    curl_setopt($ch, CURLOPT_URL, $url);
    curl_setopt($ch, CURLOPT_USERPWD, "myusername:mypassword");

     $data="<layer>
            <defaultStyle>
            <name>polygon</name>
            </defaultStyle>
            </layer>";
     curl_setopt($ch, CURLOPT_PUT, true);
     curl_setopt($ch, CURLOPT_HTTPHEADER,
                array("Content-Type: $contentType",
                'Content-Length: '.strlen($data))
            );
        curl_setopt($ch, CURLOPT_POSTFIELDS, $data);
        curl_setopt($ch, CURLOPT_HEADER, false);
        curl_setopt($ch, CURLOPT_RETURNTRANSFER, true);

        $rslt = curl_exec($ch);

       $info = curl_getinfo($ch);

?>

But this code is not working for me . There is something wrong in execute curl using PHP.


texts: Python Kivy DLL load failed on Using ArcGIS ArcPy

texts:
I am running Kivy and Python properly when it is only Kivy as

import kivy
kivy.require('1.0.6')

from kivy.app import App
from kivy.uix.label import Label

but when I add ArcGIS Python library called ArcPy like below

import arcpy
import kivy
kivy.require('1.0.6')

from kivy.app import App
from kivy.uix.label import Label

I am getting following error:

Message File Name   Line    Position

  Traceback
           5
        C:\Python27\ArcGIS10.4\lib\site-packages\kivy\app.py    327
        C:\Python27\ArcGIS10.4\lib\site-packages\kivy\uix\widget.py 219
        C:\Python27\ArcGIS10.4\lib\site-packages\kivy\graphics__init__.py  89
      init kivy.graphics.instructions (kivy\graphics\instructions.c:14640)  kivy\graphics\vbo.pxd   7
      init kivy.graphics.vbo (kivy\graphics\vbo.c:5482) kivy\graphics\compiler.pxd  1
      init kivy.graphics.compiler (kivy\graphics\compiler.c:2983)   kivy\graphics\shader.pxd    5
      init kivy.graphics.shader (kivy\graphics\shader.c:11990)  kivy\graphics\texture.pxd   3
      init kivy.graphics.texture (kivy\graphics\texture.c:31800)    kivy\graphics\fbo.pxd   5
      init kivy.graphics.fbo (kivy\graphics\fbo.c:7683) kivy\graphics\fbo.pyx   84       ImportError: DLL
  load failed: The specified procedure could not be found.

can you please let me know why this is happening and how I can fix it?


texts: How to solve error in Raster calculator in ArcGIS 10.3?

texts:
I am not sure how to solve the error message to solve my raster calculator algorithm.

The process I am trying to run is to sum pixel values of two raster files of the same area but that represent different years.

The error message is found below:

Executing: RasterCalculator ""rice_v2_2000_india_new.tif" + "rice_v2_2001_india_new.tif"" F:\Rice_Fallow_Mapping\Data_Biradarsir\India_rice\Fallows\rice_paddy\rice_xiao_v2_india_newboundary\geographic\test
Start Time: Thu Sep 15 19:00:54 2016
Raster(r"rice_v2_2000_india_new.tif") + Raster(r"rice_v2_2001_india_new.tif")
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 2, in rcexec
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 24, in <module>
    from arcpy.toolbox import *
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\toolbox.py", line 355, in <module>
    from .management import Graph, GraphTemplate
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 26, in <module>
    import _management
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\_management.py", line 14, in <module>
    import _graph
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\_graph.py", line 27, in <module>
    import numpy
  File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\__init__.py", line 180, in <module>
    from . import add_newdocs
  File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\add_newdocs.py", line 13, in <module>
    from numpy.lib import add_newdoc
  File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\__init__.py", line 4, in <module>
    from type_check import *
  File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\type_check.py", line 8, in <module>
    import numpy.core.numeric as _nx
  File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\core\__init__.py", line 20, in <module>
    import function_base
  File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\core\function_base.py", line 6, in <module>
    from .numeric import result_type, NaN, shares_memory, MAY_SHARE_BOUNDS, TooHardError
ImportError: cannot import name shares_memory

This is the output

Failed to execute (RasterCalculator).
Failed at Thu Sep 15 19:00:54 2016 (Elapsed Time: 0.03 seconds)


texts: I'm new to QGIS, and am using v2.18.3 with the cloud plugin.  I've published a map to the cloud, but am having problems with an attribute field in one of the layers.  It's a URL field, and the data contains http strings for some rows and https for others.  The http strings are properly displayed as hyperlinks and open up a new tab when clicked.  The https strings seem to be processed as normal text (no hyperlink, no action performed when clicked.)  Is there some trick to getting the software to process https strings as hyperlinks?


texts: Publishing Feature Services from spatial view on ArcGIS Server

texts:
I have a feature class layer and a handful of non spatial plain tables in a 10.1 ArcSDE geodatabase on MS SQL Server. What I want to achieve is to publish a feature service from ArcGIS Server so as for users to edit both features and attributes from the non spatial tables.

I tried to make a spatial view by joining the feature table with other non-spatial tables by referencing foreign keys. I was able to publish the spatial view as feature service but not to  make any edits from the REST service (Add, Update, Delete).

Just wondering why this was not successful and if a spatial view cannot be edited through feature service?


texts: Error when passing shapefile to PostGIS through submit button and PHP

texts:
I have the following issue when I try to pass a shapefile to PostGIS.
When I run this part of code through the terminal:

shp2pgsql -s 2100 -d /var/www/dkar/MarCO/Code/python/output/cluster/clusterNodes_.shp | psql -h localhost  -p **** -d ****** -U *******

It works perfectly and a new table in PostGIS is created, but when I try to do the same through a PHP script which I execute on a submit button then I get an error. I can not say for sure what's the error why I run the whole process through PHP exec so I just get an error code 2.

Also when I do the process through the terminal, I am asked for entering the db password. I also tried to include the password in the command like:

     `... -W password`

but it didn't make any difference.

The whole thing runs smoothly locally but it doesn't work in the remote server.

Any ideas what I might miss here?

Perhaps the path to the shapefile should be different in each case?


texts:
I recently requested a NOAA "token" to access their data. However, I keep receiving a "Token parameter is required" 404 error. Here is the URL I am using (note token is not the real one):

http://www.ncdc.noaa.gov/cdo-web/api/v2/locationcategories?Token=AbCdEfGhIjKl1234

I have been scouring google for about an hour now on just a simple example to see what I am doing wrong, but no luck. I'm hoping someone here can point out (and probably get a chuckle) what I am doing wrong.


texts: Generating a DEM and DSM from correct LiDAR point classification

texts:
I have downloaded classified LAS data from http://lidar.unh.edu/map/

I am trying to generate accurate DEM and DSM for further analysis. I created a LAS dataset with a projected CRS of  NAD_1983_UTM_Zone_19N in meters and the z CRS in meters as well. My question is on which classifications do I choose for the DEM and DSM. Here is the filter options I have in my dataset:

According to unh LiDAR data report class 2 is ground.

What I have done so far:

DEM I chose the class 2 points and used the las dataset to raster.

DSM In the predefined settings I chose first return and used the las dataset to raster tool.

Is this an accurate way of generating these two rasters? Do I not have to take into account the the unassigned class 1, noise class 7, reserved 11, reserved 17, reserved 18?

Additionally (I can ask this as a separate question if it gets requested to), when using the las dataset to raster tool the sampling value is defaulted to 10. I would like to change it to 2 or 3 (would be meters) to match the LiDAR resolution. I also plan on using the generated DEM and DSM for a least cost path analysis and I want them to be higher resolution. Will making the sampling_value 2 or 3 throw off the results or should I leave it at 10?


texts: How to use Spatial Criteria with ArcPy Data Access Search Cursor?

texts:
The ArcGIS documentation for the Data Access Search Cursor states:

  The records returned by SearchCursor can be constrained to match attribute criteria or spatial criteria.

But I can't see any example of how to specify spatial criteria when calling the search cursor.  For this example I have a geometry and wish to constrain the results of the search cursor to only features which both intersect that geometry and meet other attribute criteria.  e.g.

with arcpy.da.SearchCursor(fc,field_list,
                          "Admin_Code is null and shape intersects admin_area_geometry")

I understand that I could perform a search cursor with the attribute query and then use the geometry intersects method but this would be slow to have to loop through each feature.

How can I specify both the spatial criteria when creating a search cursor?


texts: Defining units of measure used for raster layer in ArcGIS Pro using standard ISO 19115-3 XML metadata style

texts:
I have set my metadata style to the above mentioned style in options in ArcGIS Pro before starting my project. Now I have imported multiple raster layers under my map document and each one has different units of measure and different metadata fields to be filled. When I go to the Catalog Pane for filling the metadata, i get the complete view and certain fields but its for the complete document and it does not shows any field for filling in "units of measure".
How can I assign in the metadata which units of measure have been used for a specific layer using this metadata style and how can I define metadata for individual layers?


texts: How to create a layer from existing attribute data of another layer?

texts: I am using Arc Map 10.2 and I have a shapefile consisting of parcel data.  I would like to create a shapefile/layer from the existing parcel shapefile by using a table that consists of a list of addresses that match up to the attribute associated with the parcel layer.


texts: Openlayers map WMS stretched by latitude

texts:
I am working on a small JavaScript application for displaying map with WMS layers from 2 servers. The problem is, that image is stretched by latitude, it is wider, than it should be.

This is is how it is working now:
http://postimg.org/image/pxk8lfnnd/

and this is how it should be, and how it Google maps displaying:
http://postimg.org/image/al321lg6l/

This is JavaScript code:

var map = new OpenLayers.Map("map");
var ol_wms = new OpenLayers.Layer.WMS(
"NP Podyjí WMS",
    "http://mapy.nppodyji.cz/cgi-bin/mapserv.exe?map=f:/mapserv/projects/map/nppodyji_wms.map",
{
   layers: "nppodyji_wms,hranice_npt,hranice_np,cyklo_CR",
   transparent: "true",
  format: "image/png"
},
{ isBaseLayer: false, visibility: true }
);

var ol_cz = new OpenLayers.Layer.WMS(
"Česká republika",
"http://geoportal.gov.cz/ArcGIS/services/CENIA/cenia_arccr/MapServer/WMSServer",
{layers: "0,3,4,6,7,10,12,13,14,15,16,17,19,20,21,22,23,24"}
);

map.addLayers([ol_cz, ol_wms]);
map.addControl(new OpenLayers.Control.LayerSwitcher());
map.zoomToMaxExtent();

setTimeout(function() {
var fromProjection = new OpenLayers.Projection("EPSG:102067");   // Transform from WGS  1984
var toProjection   = new OpenLayers.Projection("EPSG:900913"); // to Spherical Mercator Projection
 var position       = new OpenLayers.LonLat(15.8186 + (15.9626 - 15.8186) / 2,48.8027 + (48.8866 - 48.8027) / 2 ).transform( fromProjection, toProjection);
 map.setCenter(position, 13);
}, 100);

I tried many map.projection and map.displayProjection changes, but I'm not really sure, if it's projection problem.

It is only an example of bigger application, which take points from database and displays them dynamicaly. Displaying is correct; problem is, that client don't want it in different shape

Projection should be EPSG:102067, in HTML I have also included proj4js.js.

EDIT:
I've solved only half of the problem:
Without any display settings, it is displayed like this: http://postimg.org/image/czs81gzwj/
When I used:

var ol_cz = new OpenLayers.Layer.WMS(
    "Czech",
    "http://geoportal.gov.cz/ArcGIS/services/CENIA/cenia_t_podklad/MapServer/WMSServer",
{layers: "blablabla"},
    { projection: new OpenLayers.Projection('EPSG:3857')}
);

it looks like this: http://postimg.org/image/vm8rkxgoz/ so exactly what I want!
However problem is with another WMS server. When I was trying to change projection, or displayProjection in:

    var ol_wms = new OpenLayers.Layer.WMS(
    "NP Podyji WMS",
    "http://mapy.nppodyji.cz/cgi-bin/mapserv.exe?map=f:/mapserv/projects/map/nppodyji_wms.map",
    {
        layers: "nppodyji_wms,hranice_npt,hranice_np,cyklo_CR",
        transparent: "true",
        format: "image/png"
    },
    { isBaseLayer: false, visibility: true } ,
    { projection: new OpenLayers.Projection('EPSG:XXX'),
    displayProjection: new OpenLayers.Projection("EPSG:XXX")
}
        );

It has basically no effect :/. I was trying several projections in "projection" attribute, also in displayProjection, but only visible output was http://postimg.org/image/yv0oe22tf/ (note the pink square on the bottom). So it get some parameters, but they are wrong, or maybe he server is wrong, I'm not really sure.

The 2nd WMS did not fit in the baselayer (1st WMS)


texts: Configure a layer in GeoServer so it is displayed with a transparent background?

texts:
I use GeoServer to serve shapefile data as layers, and a program called 'Sitmun' to display them. Sitmun allows for the superposition of several layers. But, since my GeoServer seems to be delivering the layers with an opaque background, superposition is not working as desired. Layers taken from a governmental service – which I do not manage – have got a transparent background and work just fine.
How can I configure my GeoServer to serve layers with transparent backgrounds?


texts: ArcMap 10.6; ERROR 000539: Error running expression: rcexec()

texts:
While running a tool box (connectivity index) I am receiving an error at the end; ERROR 000539: Error running expression: rcexec()!.
Can anyone help me with it?


texts: Calculating Distance Between Two Give Points using ArcGIS API for JavaScript 3?

texts:
Using ArcGIS API for JavaScript 3.22 and having two given Points as

var pt1 = new Point(-106.61, 35.1107, new SpatialReference({ wkid: 4326 }));
var pt2 = new Point(-108.61, 35.1107, new SpatialReference({ wkid: 4326 }));

can you please let me know how I can calculate the distance from pt1 to pt2?

var map;
    require([
             "esri/map",
             "esri/geometry/Point",
             "esri/SpatialReference",
             "esri/geometry/webMercatorUtils",
             "dojo/domReady!"
 ], function(Map, Point, SpatialReference, webMercatorUtils) {

 map = new Map("map", {
       basemap: "topo",
       center: [-106.61, 35.1107],
       zoom: 8
 });
 var geometryService = new esri.tasks.GeometryService("https://utility.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer");

 var pt1 = new Point(-106.61, 35.1107, new SpatialReference({wkid: 4326 }));
 var pt2 = new Point(-108.61, 35.1107, new SpatialReference({wkid: 4326 }));

});


texts: Is there any way to lock layers in QGIS to prevent them from being accidentally deleted like you can in Illustrator or CAD?


texts: Raster and vector layer did not overlay with each other even in the same CRS

texts:
I am a newbie in the GIS world, I hope this community can help me to improve my skill.
However, their coordinate is different how to change their coordinate to be in the same value?
How can they have different coordinates?
raster data(image for raster data for reference)

These are the information about the layers
vector

raster


texts: Python shapely intersection with buffer in meter

texts:
Having a line and a polygon, I use shapely to determine the intersection between the two with a defined buffer. Both the line and the polygon are in WGS84 coordinates, while the buffer is in meter.

from shapely.geometry import LineString, Polygon
polygon = Polygon([
    (51.9981816, 4.3677729),
    (51.9983431, 4.3683897),
    (51.9980583, 4.3685846),
    (51.9977777, 4.3687773),
    (51.9976163, 4.3681605),
    (51.9978555, 4.3679964),
    (51.9981816, 4.3677729)
])
line = LineString([
    (51.9976836984107, 4.368050461186726),
    (51.99757198557364, 4.368489603711765),
    (51.9973971301288, 4.3684738261236475)
])

Finding the intersection is easy:

line.intersection(polygon)

Finding the intersection with a buffer in degrees (because coordinates of the objects are in degrees):

line.intersection(polygon.buffer(0.2))

But how to find the intersection with a buffer in meters? I understand that Shapely is not aware of any coordinate reference system; also no area metric CRS is known (OceanoGraphy is globally).


texts: Segmentation fault while writing an output raster by blocks with rasterio

texts:
With the answer on my previous question, I was able to set up reading a raster by blocks. I did some computation and now I try to write the result into another raster. However, the program fails with error code 139 and Segmentation fault (core dumped). My code:

with rasterio.drivers():
    with rasterio.open(asterPath,'r') as src:
        with rasterio.open(ndviPath,'w',driver='GTiff', width=src.shape[0], height=src.shape[1], count=1, dtype=np.float32) as dst:
            for index, window in src.block_windows(2):
                red = src.read(2, window=window)
                nir = src.read(3, window=window)
                #compute ndvi as NumPy 2d array, dtype float32
                dst.write_band(1,ndvi,window=window)

As a sanity check, I tried to fill ndvi array with zeros. It worked outside the loop, but not within, so I must have misunderstood something about writing by block with rasterio. All I know about rasterio is this manual (including the other doc files in the same repository). So how to save the NDVI data to a raster by block?

I'm on Ubuntu 14.04 and I use Python 2.7, Gdal 1.10.1, NumPy 1.8.2 and Rasterio 0.26.0 installed from source the way suggested for Linux at the Rasterio site.


texts: How to calculate the river slope for a lot of rivers automatically in ArcMap?

texts:
I'm working with ArcMap 10.2 and I'm pretty new in this area..
What I want to do is compute the river slope for a lot of rivers from a DTM...I know that there is a way to do it manually passing by excel..but I was wondering if there was an automatic way to do it..also maybe using other programs..
I have all the territory of Swizerland to analyse..


texts: QGIS plugin with additional Button and function

texts:
I will want to create a simple QGIS plugin where the user put a number (line_edit) and if that number is something the user take back three new codes(lineEdit_2,lineEdit_3,lineEdit_4) and if user press ok(buttonBox)  then plugin do something else in the some GUI.
For that I create a new additional pushButton in Qdesigner and I create additional def function(run1) for execute code for pushButton.
First I add new pushbutton in __init__ :
def __init__(self, iface):
    self.dlg = pluginDialog()
    self.dlg.pushButton.clicked.connect(self.run1)

My gui:
[![qdesigner][1]][1]
Plugin work without error but don't show me the codes (lineEdit_2,lineEdit_3,lineEdit_4) in the plugin window if I press pushbutton,any idea ?
My python code :
def run1(self):
    distance = str(self.dlg.lineEdit.text())
    if distance==0:
        area1='some'
        area2='hello'
        area3='world'
    self.dlg.lineEdit_2.setText(str(area1))
    self.dlg.lineEdit_3.setText(str(area2))
    self.dlg.lineEdit_4.setText(str(area3))

def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something else
        pass

And I try to replace run1 with:
def run1(self):
    self.dlg.show()
    result = self.dlg.exec_()
    if result:
        distance = str(self.dlg.lineEdit.text())
        if distance==0:
            area1='some'
            area2='hello'
            area3='world'
        self.dlg.lineEdit_2.setText(str(area1))
        self.dlg.lineEdit_3.setText(str(area2))
        self.dlg.lineEdit_4.setText(str(area3))
        pass

but no change.
[1]: https://i.stack.imgur.com/kTMMZ.jpg


texts:
Is there any place where we can get soil Shapefiles of India (if not detailed information, can we get region wise/state wise)?
I  have seen this data base for few regions like USA.

To compute "Q" in a basin or at a point, I am looking for soil shapefiles to use in Arcgis.


texts: Select by Expression and Extract by Expression generate different results

texts:
I have a polygon layer and want to extract all Polygons with the field "street" = 1.
When I use "Select by Expression" I get the following correct result:

However, when I use "Extract by Expression" I get the following wrong results:

When I look into the attribute table of the extracted layer all values under the field "street" = 1 but a couple of them shouldn't be.
What could be the problem?
I'm using QGIS 3.10


texts: How to change the "application_name" of a PostGIS-connection in QGIS?

texts:
PostGIS offers a view called pg_stat_activity, where you can see all the active connections to the database. This view includes the interesting column application_name. When a connection comes from a QGIS instance, the text "QGIS" appears there.

An application can send this parameter to PostGIS in the connection string, and I suppose that QGIS do it so. The connection string looks than more or less like that:

host=... port=... dbname=... user=... password=... application_name=QGIS

I would like to change this value in QGIS to include the username (that I will get from the operating system). Instead of "QGIS", I would like to see something like "QGIS for user1". Is there a was to do that, may be by changing the connection string using a pyQGIS script or any other solution?

Edit:

I just tried following thing: using the python console I made use of the function setEncodedUri which seems allow all possible connection parameters. But it doesn't work: The application_name appears in the layer properties in QGIS but in PostGIS I still just see "QGIS" and nothing more. This are the lines I was using for this test:

    uri = QgsDataSourceURI()
    uri.setEncodedUri("host=...&port=...&dbname=...&user=...&password=...&application_name=QGIS for user1")
    uri.setDataSource("MySchema", "MyTable", "wkb_geometry", "")
    vlayer = QgsVectorLayer(uri.uri(), "MyLayerName", "postgres")
    QgsMapLayerRegistry.instance().addMapLayer(vlayer)

Edit2: I had a look in the source code of QGIS 2.8: the trick with setEncodedUri can't work, because the code calls SET application_name='QGIS' just after opening the connection and that overwrites my setting!

I'm still looking for a working solution...

Edit3

It possible to set the environment variable PGOPTIONS to something like "-c application_name=QGISuser1" or even easier the environment varaible PGAPPNAME (see fallback_application_name) before starting QGIS. It's used to change some postgres configuration variables for the current session. I tried it with over variables and it worked. But as QGIS overwrites the value of application_name after opening the connection, it doesn't work for this variable. And I'm now thinking, that there are probably no ways to do what I want with the current version of QGIS.


texts: Locating a layer to download for my project

texts: I'm trying to find a layer that I can download that will allow me to show the administrative areas of Scotland. I've found plenty sources for downloads but can't seem to find what i'm looking for. This is purely an issue for me as I'm so new to GIS, it is absolutely mind boggling downloading government datasets.


texts: How to set 2 vector layers into map canvas, pyqgis 3.0

texts:
I have 2 shapefiles which I need to lap. The first one is city grid, and the second one representing borders of the city. Anyway, my code is:

extent = QgsRectangle()
layers = []
QgsProject.instance().addMapLayer(layerInput)
extent.combineExtentWith(layerInput)
layers.append(layerInput)
QgsProject.instance().addMapLayer(layerInputShape)
extent.combineExtentWith(layerInputShape)
layers.append(layerInputShape)
canvas.setExtent(extent)
canvas.setLayers(layers)

But for some reason it's not working, I didn't get the error message. :/


texts: Point disappears when zooming in - QGIS 3.4.1

texts:
I had a vector layer stored in GeoPackage format with a text field (lbl_text) that was being used to create labels on a map. After creating 3 points I decided it would be nice to rotate my labels, so I added a rotation field (lbl_rotation). This prompted multiple notifications stating that the rotation field already existed and that the changes could not be made. I could not save changes to the layer so I cancelled and discarded them.

Despite the rotation field not previously existing and the notifications that creating this field had failed, the rotation field was now present. So I linked it to the data defined rotation option under Properties > Labels > Rotation. Now my labels could be rotated using Label Toolbar > Rotate Label.

But when zooming in the point created after this event disappeared. I could not see any reason for this and checked for others having a similar issue. A common solution was to use the geometry checker, so I tried this and found no errors (they're all points after all, what could go wrong?).

Deleting the point and recreating it did not fix the issue. The 4th point created always disappeared when zoomed in, regardless of how many times it was deleted and recreated, or where on the map it was created.

I added a 5th point and it did not have this issue. Eventually I deleted all points and recreated them all. None of them had this issue.

Edit: now all points on this layer exhibit the same problem behaviour, though they didn't when first posting about the issue.

Does anyone have any idea what might have happened and how to fix the problem?

Edit: The following has recovered the project so has been noted to help others, but it does not identify or repair the original problem so I am not marking it as the answer. Note that simply duplicating the affected layer resulted in duplicating the problem too.

To recover from this problem:
 Right click the affected layer and save as a new GeoPackage.
 Ensure Add saved file to map is checked or re-add the saved layer.
 Right click the affected layer and select Styles > Copy Style.
 Right click the new layer and select Styles > Paste Style.
 Delete the affected layer and rename the new one.

Points when zoomed out.

The same 3 points when zoomed in. One cannot be seen.


texts: Accelerate neighbor analysis in QGIS PostgreSQL / PostGIS

texts:
For a project I have to show the relationship between 2 or more objects.
For this I use the code below in QGIS in step 1 I determine the relation in step 2 the items are counted.
In the test set I use 5128 objects here QGIS takes quite a long time, the original data set is more than 1 million objects.
I wonder would this also be possible in PostgreSQL / PostGIS to make the whole faster and reproducible?

Code for QGIS 3
Step 1

aggregate(
 layer:= ‘Test',
 aggregate:='concatenate',
 expression:=Label,
 concatenator:=', ',
 filter:=touches($geometry, geometry(@parent))
 )

Step 2

aggregate(
 layer:= ‘Test’,
 aggregate:='count',
 expression:=Label,
 filter:=touches($geometry, geometry(@parent))
 )

product

In the result you can see that an (label) object has a relation with an other object 1 to N=(Count)


texts: Script to show what transformation is being used in ArcMap

texts:
Does anyone know of a script that will search through a folder of MXDs (or an individual MXD) and report what coordinates system the data layers are in, what coordinate system the data frame/map extent is in, and if different state what transformation has been used (if any)?

Would really like to check content and whether correct transformations are being applied without having to open up every MXD for a project


texts:
I have been given a list of locations within Denmark. Most of the are in decimal degrees. However, some of the positions are in different formats and I can't figure them out. Does anyone recognize these or have a guess as to what they might be?

I have some that look like this:

lat: 55,800,643,991,678,900 lon: 12,297,413,013,875,400
lat: 55,804,257,011,041,000 lon: 1,229,334,301,315,240
lat: 55,807,527,964,934,700 lon: 12,333,790,976,554,100

and some that look like this:

lat: 583625 lon: 6236347
lat: 583722 lon: 6236625
lat: 562843 lon: 6218293


texts:
I had create a table and line feature class in ArcCatalog.
i want to deal with those tables in ArcMap.

Read
Add
Delete
Update

how could i do that?


texts: Graph of the variation in time of the nival surface

texts:
How do I calculate the area of  snow for all the images in a collection image? When I apply the script it only calculates me from an image. And after  to generate a chart
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
                .filterBounds(table)
                .filterMetadata('CLOUD_COVER', 'Less_than', 20 )
                .filter(ee.Filter.eq('WRS_PATH', 232))
                .filter(ee.Filter.eq('WRS_ROW', 77))
                .map(function(image){return image.clip(table)})
var ndsi = landsat1.normalizedDifference(['B3', 'B6'])
var mask_ndsi = ndsi.gt(0.4)
var mask = ndsi.updateMask(mask_ndsi)
var area = mask.multiply(ee.Image.pixelArea()).reduceRegion({
reducer: ee.Reducer.sum(),
geometry:table,
scale:30
})


texts: ArcGIS JS API Polygon's WKID conversion with gsvc.project

texts:
I'm trying to convert polygon from wkid 3301 to wkid 4326. I've been using this code snippet for converting points, but I'm not sure how I would accomplish this to create a polygon object. I'm not even sure if gsvc.project() is the best solution. Any suggestions would be appreciated.

dojo.require("esri.tasks.geometry");
var outSR = new esri.SpatialReference({ wkid: 4326 });
mapHandler.gsvc.project([ geometry ], outSR, function(projectedPoints) {
    pt = projectedPoints[0];
});


texts: Google Earth Engine - Dynamic Visualization Parameters

texts:
I made a script that calculates the maximum NDVI values of the last 100 days of an area. Therefore, depending on WHEN you run it, the resulting image.
Is there any straight straightforward method to achieve this?
I tried a image reducer and then parse it to the visualization parameters, but can't make it work.

Best would be to get real-time "Stretch by XXX standard deviations", or "min-max".
How do I enter these values to the visualization parameters?

var newDate = new Date();
var dateto = ee.Date(newDate);
var datefrom = dateto.advance(-120, 'day');

var geometry = ee.FeatureCollection('TIGER/2018/States').filterMetadata('NAME','equals','Nebraska');
var collectionNDVI = ee.ImageCollection('MODIS/006/MOD13Q1')
                .filterDate(datefrom ,dateto)
                .filterBounds(geometry)
                .select("NDVI");
var aoi= geometry.geometry().bounds()
var NDVI = collectionNDVI.select(['NDVI']).max().clip(aoi)

var STD= NDVI.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry: NDVI.geometry(),
  scale: 250,
  maxPixels: 1e9
});
print(STD)
var nbrParams = {min: 1, max: 5000 , palette: ['red', 'yellow', 'green']};
Map.addLayer(NDVI,nbrParams, "maxNDVI COMBINED")

Thanks,


texts: What's the difference in these two shapefiles?

texts:
I'm trying to plot points on these maps in R, with points from "Coordinates" in the same folder.

The WL.shp plots fine with the data, but the NG.shp (which should only be larger and of the same area + surroundings) won't show the dots?

I tried different UTM zones and GIS programs (MapInfo and QGIS).
Can anyone help to what the issue is, what am I missing to make this work?

## appropriate library.
library(StatDA)
library(rgdal)
## load data
setwd("C:\\...")
SS_data <- read.csv2("Coordinates.csv",header=TRUE)
#WL.shp file
shp.greenland<-readOGR("C:\\Users\\Andreas\\SkyDrive\\Speciale\\02. Data\\02. Statistik","WL")
## NG.shp file
shp.greenland<-readOGR("C:\\Users\\Andreas\\SkyDrive\\Speciale\\02. Data\\02. Statistik","NG")
## cordinates
Y=SS_data[,"Latitude"]
X=SS_data[,"Longitude"]

el=SS_data[,"Xx"]
## Create pdf
pdf("Au dot plot.pdf",width=50,height=20,onefile=F)
## plot map with x,y coordinates
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plot(shp.greenland)#, col="grey")

# take symbols continuously growing
q=c(min(el),min(el)+(max(el)-min(el))/2,max(el))

# symbols and size
symb=c(16,16,16)
ssize=c(5,10,20)

# plot symbols:
el01=(el-min(el))/(max(el)-min(el))
el01ssize=el01*(ssize[3]-ssize[1])+ssize[1]
points(X, Y, pch=symb[3], cex=el01ssize, col="red")

#Legend
legn=round(q,2)
leg=rep(NA,length(q))
for (i in 1:length(q)){
  leg[i]=roundpretty(legn[i],2)
}
legend("topright",pch=rev(symb),pt.cex=rev(ssize), legend=rev(leg), title=paste("La ppm"),
       cex=8, col="red")
dev.off()


texts: OpenLayers 2 with different projection is giving blank page

texts:
The code does not give any exception but nothing appears in the browser. This code works with EPSG:4326. Is this the correct way to define projection?

See below for the full code:

function init()
{
Proj4js.defs["EPSG:32645"] = "+proj=utm +zone=45 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
            var projection = new OpenLayers.Projection("EPSG:32645");
    var bounds = new OpenLayers.Bounds(
                          371399.417972, 2251641.911461,
                          384784.474172,  2236278.339933
                     );

        var options = {
            maxExtent: bounds,
            restrictedExtent:bounds,
            displayProjection: projection,
            allOverlays: false,
            projection: projection,
            units: 'm'
        };

var  map = new OpenLayers.Map('map', options);

var ccounties = new OpenLayers.Layer.WMS(
                        "Countries02Plugin",
                         "WmsHandler.axd",
                        {
                            srs: projection,
                            version: '1.1.1',
                            layers: 'BHUBANESWAR',
                            format: 'image/png'
                        }
        );

map.addLayer(ccounties);
map.addControl(new OpenLayers.Control.Navigation());
map.addControl(new OpenLayers.Control.Scale());
map.addControl(new OpenLayers.Control.MousePosition());
map.zoomToMaxExtent();
}


texts: Distance between a Point and Linestring: PostGIS/GEOS vs. Boost

texts:
I am currently using the GEOS library through OGR and data from a PostGIS database to work with geospatial data. I now found myself in the need to use some of the algorithms Boost.Geometry and Boost.Graph offers, so I've adopted the OGR classes to boost's iterator/range concept. However, the two libraries yield different results, and I am unsure why. From the compiler perspective everything works; however, a simple test case fails.

Step by step, here is what I did.

I define two linestrings and a point and assign WGS84 as projection, in PostGIS:

ST_GeomFromText('LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0 4.0)', 4326),
ST_GeomFromText('LINESTRING(15.0 15.0, 16.0 16.0)', 4326)
ST_GEomFromText('POINT(15.0 15.0)', 4326)

Then I calculate the distance between the first and the second linestring and the distance between the first linestring and the point:

db=# SELECT ST_Distance(
    ST_GeomFromText('LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0 4.0)', 4326),
    ST_GeomFromText('POINT(15.0 15.0)', 4326));

   st_distance
-----------------
 15.556349186104

db=# SELECT ST_Distance(
    ST_GeomFromText('LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0 4.0)', 4326),
    ST_GeomFromText('LINESTRING(15.0 15.0, 16.0 16.0)', 4326));
   st_distance
-----------------
 15.556349186104

Using OGR/GEOS to do the same calculations in the client program yields the same results.

Now I compare those with the result of the boost adaption. I registered OGRPoint with the boost library with the following code:

BOOST_GEOMETRY_REGISTER_POINT_2D_GET_SET(
    OGRPoint,
    qreal,
    boost::geometry::cs::cartesian,
    getX,
    getY,
    setX,
    setY)

The most important part is the fourth line, where the coordinate system is defined. Since all my coordinates are in WGS84 projection, I chose "cartesian" as the appropriate coordinate system.

And here is the relevant code from my unit test:

OGRLineString ls1, ls2;
OGRPoint distantPoint(15.0, 15.0);
OGRSpatialReference *sref = new OGRSpatialReference();

OGRErr rc = sref->importFromEPSG(4326);
QVERIFY(OGRERR_NONE == rc);
ls1.assignSpatialReference(sref);
ls2.assignSpatialReference(sref);
distantPoint.assignSpatialReference(sref);

for (int i = 0; i != 5; ++i) {
    ls1.addPoint(static_cast<qreal>(i), static_cast<qreal>(i));
}

ls2.addPoint(15.0, 15.0);
ls2.addPoint(16.0, 16.0);

QCOMPARE(
        bg::crosses(ls1, distantPoint),
        ls1.Crosses(&distantPoint));
QCOMPARE(
        bg::intersects(ls1, distantPoint),
        ls1.Intersects(&distantPoint));
QCOMPARE(bg::distance(ls1, ls2), ls1.Distance(&ls2));
QCOMPARE(
        bg::distance(
            ls1,
            distantPoint),
        ls1.Distance(&distantPoint));

delete sref;

However, using boost's boost::geometry::distance yields totally different results: For the distance linestring to linestring, the result is 0, for point to linestring, 16.9706. Using one of the other strategies defined in the boost library also yields completely different results, such as 0.294082 for boost::geometry::cs::geographic.

What is the reason for this? What am I doing wrong?


texts: Failed to create memory layers in QGIS application on Linux

texts:
I have written a python class which reads properties of map layers from an XML file and creates appropriate memory layers. I have developed and tested this code with QGIS 2.8.2 on a Windows 7 machine (importing the class at the python console and calling the class methods). As the module should be used in a Linux environment I did the same successfully with QGIS 2.8.2 on CentOs 7 - all map layers were successfully created and symbolized.

Then I refactored the code to fit in the framework of a QGIS standalone application. There are major problems with threads and such stuff, but the app starts and comes up with all layers created, and working methods.

Now comes the strange part. In the standalone CentOs version the creation of memory layers doesn't work - layers are not valid. Even a simplified test without field definitions, no crs doesn't work.

layer = QgsVectorLayer(type + '?crs=EPSG:' + str(epsg) + fields, layer_name, 'memory')

Hundreds of lines of code before and after this line are the same in all versions. Has anybody an idea how to handle this problem? Maybe a bug?


texts: Define precision in geometric predicate methods in PyQGIS

texts:
It is possible to define precision for geometric predicate methods, such as:

intersects (), contains (), equals (), touches (), overlaps (), within (), crosses()

For example, I'd like to set the precision in this situation:

feat_to_select.geometry().intersects(selected_feat.geometry())


texts: Functionality of Attribute Assistant and ArcGIS Online

texts:
We are an EGIS group using Arc Desktop 10.3.1 and ArcGIS Online as well as the collector app. With ArcDesktop we use the ArcGIS Local Govt Attribute Assistant

We are seeking to allow users to enter service request point data using either the collector app or ArcGIS online yet still retain the functionality of the Attribute Assistant we are using in the desktop version.

Specifically we wish to use the AutoIncrement Number option to give each new service request that has been dropped in the field using the collector or arcgis online a sequentially incremented service request number.  Additionally when a point is dropped in the field using the app or online we would like the zip, map page, map block, zones, etc  automatically pulled through and populated so the field user only needs to drop a point and enter the problem. (Basically we would like Attribute Assistant to behave in ArcGIS Online/Collector just the way we have it behaving in ArcDesktop)

I have done a decent amount of research here and located this article
http://blogs.esri.com/esri/arcgis/2015/01/29/collector-for-arcgis-v10-3-is-out-and-its-a-big-deal/

that seems to say that in 10.3 there will be this functionality, but I can not find out how any one has implemented them in the past.  Is this true that is functionality exists? and has anyone had any experience using Attribute Assistant in the Collector or ArcGIS Online environments? And if so how did you go about it?


texts: Average road distance for each county between two location/Lat-Long

texts:
Say, one travels from Times Square (1451 Broadway New York, NY, 10036) to White House (1600 Pennsylvania Ave NW, Washington, DC 20500). For this 4+ hour drive, I would like to calculate the average road distance that needs to be driven in each COUNTY along the way.

Note:

A road can cross the same state or same county couple or more times
We can consider the default recommended road provided by Google maps API

I would like to have a recommendation on the approach to calculate the average distance in each county.


texts: Selecting attributes which have comma delimiter

texts:
I'm currently working on a large dataset of roads. In which i have to check if the roads in 1:250 000 have the same road number as those in 1:50 000. Because of alot of generalization I have to make a radius serch in 1:250 000 to be sure i have gotten the road into the spatial join. For this to work, i used field mapping parameter "join" with "," as a delimiter. Since this doesen't work with numbers, I converted the field to text.

The problem I'm having now is that I cannot select specific numbers from the new field. It will only work on the first, so when i try to check if the road number from 1:250 000 is equal to 1:50 000 it will only go through the 1st number in the column and not the rest that follows after the comma delimiter. Does anyone have any idea how i might check all the numbers that has been converted to text, against one specific number?


texts: Concatenating numpy arrays based on their relative positions as defined by associated Affine transforms

texts:
I have split a raster into 4 numpy arrays to feed into machine learning algorithms.
I would now like to join to results of that process back into the full array.

For each I have a numpy array, and an affine transform specific to it.
Is there a good way to join them and keep them just as a numpy array? I dont want to go making rasterio objects unnecessarily. Currently I just put them in the right order in a list, and concatenate by row, then column (I know which elements in the list should be the top row and bottom row).
But I would like to do it in a way that uses the affine transform to reconise which array should go where.
In this way I could input any number of contiguous numpy arrays into my function and it would return the concatenated aray...

Currently I have this list called rf_classification_list, that has a tuple for each subarea with the numpy array and the Affine:

[(array([[1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       ...,
       [1, 0, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [0, 0, 0, ..., 1, 1, 1]], dtype=uint8), Affine(20.0, 0.0, 434400.0,
       0.0, -20.0, 6888740.0)), (array([[1, 0, 0, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       ...,
       [0, 0, 0, ..., 1, 1, 1],
       [0, 0, 0, ..., 1, 1, 1],
       [0, 0, 0, ..., 1, 1, 1]], dtype=uint8), Affine(20.0, 0.0, 442540.0,
       0.0, -20.0, 6888740.0)), (array([[0, 0, 1, ..., 1, 0, 0],
       [0, 0, 1, ..., 1, 0, 0],
       [1, 1, 1, ..., 1, 1, 0],
       ...,
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1]], dtype=uint8), Affine(20.0, 0.0, 434400.0,
       0.0, -20.0, 6873640.0)), (array([[0, 0, 0, ..., 1, 1, 1],
       [0, 0, 0, ..., 1, 1, 1],
       [1, 0, 0, ..., 1, 1, 1],
       ...,
       [1, 0, 0, ..., 1, 0, 0],
       [1, 0, 0, ....

I send it to this function:

def merge_arrays(dataList):
    """" uses numpy concatenate function to join the arrays.
    Would be better to find a way that takes into account the affine transformations...
    Should return a tuple with array and transformation to be consistent with inputs..."""
    myArrays = [i[0] for i in dataList]
    top_row = concatenate(myArrays[0:2], axis=1)
    bottom_row = concatenate(myArrays[2:4], axis=1)
    full = concatenate((top_row, bottom_row))
    return full

Thankyou for any suggestions.


texts: xml of a layer on geoserver is different when I am opening it in browser

texts:
xml of a layer on a geoserver is different when I am opening it in browser, for example:
http://mygeoserver/geoserver/rest/services/wfs/transforms/layer_test.xml
it look like this:
<transform>
<name>layer_test</name>
<sourceFormat>text/xml; subtype=gml/2.1.2</sourceFormat>
<outputFormat>HTML</outputFormat>
<fileExtension>html</fileExtension>
<xslt>layer_test.xslt</xslt>
</transform>

and when I open xml file in notepad it looks like:
<transform>
    <sourceFormat>text/xml; subtype=gml/2.1.2</sourceFormat>
    <outputFormat>HTML</outputFormat>
    <outputMimeType>text/html</outputMimeType>
    <fileExtension>html</fileExtension>
    <xslt>layer_test.xslt</xslt>
    <featureType>
        <id>schema:layer_test</id>
    </featureType>
</transform>

I'm missing part with the featuretype. What can be cause of it, am I searching on wrong place?


texts: OpenLayers how to style features when importing from GeoJson file?

texts:
I'm attempting to style features based on their attributes from a GeoJSON file, however I'm having difficulties

Styling code:

        //Style
    var styleFunction = function (feature) {
        var color;

        if (feature.get("Crime Type") === 'burglary') {
            color = "green";
        } else {
            color = "black";
        }
        var retStyle = new ol.style.Style({
            stroke: new ol.style.Stroke({
                image: new ol.style.Circle( /** @type {olx.style.IconOptions} */({
                    radius: 5,
                    fill: new ol.style.Fill ({
                    color: color,
                })

                }))
            })
        });
        return retStyle;
    };

Adding the markers/features to the map

       var vector = new ol.layer.Vector({
           source: new ol.source.Vector({
               url: '2017-01-west-yorkshire-street.json',
               format: new ol.format.GeoJSON()
           }),

           style: styleFunction
       })
       map.addLayer(vector);

However this just renders the map with 0 markers, if I remove the line "style: styleFunction" the map renders with default styling (all blue)

In case you couldn't tell, I'm very new to JS and Openlayers.


texts: How to style graduated using size-scale field

texts:
Is it possible to style based on size-scale field using graduated symbol but instead of using color ramp, I would like to display its size. I would also like to display this in the legend in the map composer. It would be nice if I could also apply this if I'll use svg marker (example 2nd image below).

From this:

to something like this:


texts: How to access tiles in Tilecache using QGIS?

texts:
I've Mapserver set up with a mapfile that points to a tileindex. In QGIS I can call the mapfile and view the rasters dynamically. I've got access to tilecache and I've seeded the same mapfile/rasters that I view in QGIS. My question is, what do I need to do to get QGIS to access the tilecache as opposed the dynamic rasters?

I'm guessing that I need to do either something to my mapfile to get it to point to the seeded rasters as opposed to the dynamic rasters or repoint the WMS in QGIS


texts: Recalculating Land Size from RestAPI service in AGOL

texts:
I am putting together a map in the web based AGOL map creator that would show population density using a restAPI from the US Census. The field that contains the area is called AREALAND and is calculated out in square feet and I am interested in converting this into square miles. I have tried to add a calculated field in AGO, but the table options don't allow me to create a field since it's an read only file. I could calculate out the density in the pop-up, but I can't get it to show the density on the map. Is there another way to create a virtual field?

UPDATE
So I found the query link at the bottom of the page and I am assuming that this is where I could do a simple query (AREALAND / 27878400) , but I am not sure where this calculation needs to be done or the formatting of the query (VBA JS?).


texts: PostGIS - How to do ST_UNION on each feature of a geojson FeatureCollection

texts:
I have a simple geoJSON featureCollection composed of 2 features:

featureCollection = {u'type': u'FeatureCollection', u'features': [{u'geometry': {u'type': `u'MultiPolygon', u'coordinates': [[[[-184557.61264043, 384896.54253906], [-134666.391073443, 239616.414560895], [-308616.222736376, 238788.813082666], [-184557.61264043,384896.54253906]]]]}, u'type': u'Feature', u'properties': {u'feat_id': 0}}, {u'geometry': {u'type': u'MultiPolygon', u'coordinates': [[[[-313889.06448095, 334448.10387786], [-206265.7286704, 279107.6954071], [-249376.21261724, 359825.19726501], [-31388`9.06448095, 334448.10387786]]]]}, u'type': u'Feature', u'properties': {u'feat_id': 2}}]}

I would like to do a ST_UNION on those features in order to have only one feature merged.

SQLQuery = '''SELECT ST_AsGeoJSON(ST_Union(%s))::json As f''' % (featureCollection)

First, the ST_GeomFromGeoJSON() function accepts only the geometry part of the geoJSON as input, so the featureCollection gives an error.

Do I have to loop over the featureCollection to extracts each feature's geometry or is there another solution ?


texts: enable Advanced Digitizing Panel (cadDockwidget) with PyQgis

texts:
I am working on a PyQgis Plugin where I want to enable cad in advanced digitizing panel. Following pre-conditions are met:

layer is in edit mode
project has projected crs
move feature tool is activated (iface.actionMoveFeature())

I found that the way to activate cad tool is the following way:
iface.cadDockWidget().enable()

But it does not work reliably in my case. Sometimes the cad tool is enabled and sometimes not. For the cases it is not activate with the python command, I was always able to enable it by clicking the button in the panel, which means there should not be a problem which prevents it from enabling.

Any idea what the reason could be?


texts: Embedded images in KMZ not shown on Google Earth for iOS

texts:
I've created KMZs that work well in Google Earth. Their placemarks reference embedded images files like so

<Placemark>
<name><![CDATA[IMG_1288]]></name>
<Snippet maxLines="2"><![CDATA[]]></Snippet>
<description><![CDATA[
]]></description>
<Point>
<coordinates>174.812836,-36.826668,0</coordinates>
</Point>
</Placemark>

The KMZ structure is

doc.kml
images/
       IMG_1288.jpg

The same KMZ file does not work on iOS though. The popup window opens but it displays a broken image icon. What do I need to do to embed an image in a KMZ so that Google Earth for iOS can display it?


texts: Why do I get a 503 error when pushing data to the CartoDB SQL API

texts:
I'm using ogr2ogr to push a large (200,000 feature) dataset to the CartoDB API. I'm using where_clauses to upload 500 features at a time, and sleeping every 20,000 features for a minute to give the server a rest.

Even after implementing this batched behavior, I still occasionally get a 503 error. Are there best practices for pushing data to the API, regarding sleep intervals/number of features to push at once?

I'd like to replicate the shapefile sync capability in cartoDB, but my datasets are too large for shapefiles, and so am using this method instead.


texts: How to create mask bands in GDAL and Python?

texts:
I have small program in Python that uses gdal and numpy to perform a logarithmic scale of some data into the RGB 0-255 range.

Now this input data has a nodata value set to mask the valid values, and I would like to have the output to retain that mask.

Obviously I can't use nodata for my output as I am ussing the full byte range for colour values.

From what I can tell from the GDAL API I need to create a mask band in the output file.

I added the code to do this but suddenly the program started running really slow and the output file size became huge.

What it the best way to do this?

Here is my code:

#!/usr/bin/env python

import gdal
from gdalconst import *
import numpy
import sys

infile = sys.argv[1]
outfile = sys.argv[2]
format = 'GTiff'
type = GDT_Byte

indataset = gdal.Open(infile, GA_ReadOnly)

gdal.SetConfigOption('GDAL_TIFF_INTERNAL_MASK', 'YES')
out_driver = gdal.GetDriverByName(format)
outdataset = out_driver.Create(outfile, indataset.RasterXSize, indataset.RasterYSize, indataset.RasterCount, type)
outdataset.CreateMaskBand(gdal.GMF_PER_DATASET)

gt = indataset.GetGeoTransform()
outdataset.SetGeoTransform(gt)

prj = indataset.GetProjectionRef()
outdataset.SetProjection(prj)

for iBand in range(1, indataset.RasterCount + 1):
    inband = indataset.GetRasterBand(iBand)
    inmaskband = inband.GetMaskBand()

    outband = outdataset.GetRasterBand(iBand)
    outmaskband = outband.GetMaskBand()

    for i in range(inband.YSize - 1, -1, -1):
        scanline = inband.ReadAsArray(0, i, inband.XSize, 1, inband.XSize, 1)
        scanline = numpy.multiply(numpy.divide(numpy.log1p(numpy.minimum(10000, numpy.maximum(0, scanline))), numpy.log1p(10000)), 255)
        outband.WriteArray(scanline, 0, i)

        scanline = inmaskband.ReadAsArray(0, i, inband.XSize, 1, inband.XSize, 1)
        outmaskband.WriteArray(scanline, 0, i)


texts: Scripts to automate conversion from DEM to contour shape files

texts:
I am using the DEM files in the Cornell CUGIR repository to build contour layers in QGIS.  The problem I have is that while it is easy to do, it requires many manual steps in QGIS that are both tedious and prone to error.  I would like to know if there are any tools out there to automate the process.

Starting with a directory of DEM files, on a Unix box (in my case a Mac), here is what I would like to do:

convert the DEM files from metric to english units (meters -> feet)
create a contour layer with 100 ft intervals including labels for the elevation
create a contour layer with 20 ft intervals including labels for the elevation
apply the same style to all the created contour shape files

The tool I envision would either find all the applicable DEM files in a directory or take the list of files as a command line argument.  I will write this myself if there isn't one out there somewhere already.


texts:
I have obtained some boundary maps and have joined population data to the corresponding areas using the join tool in ArcGIS. Before the joining, I first checked and edited both the boundary map and population data to synchronize the two sets of data.

For instance, an area A may be represented by a population X but may have two polygons. In this case I merged them into one. Editing produced community boundary maps with their matching population and population density. I then exported the population data which was in Excel/CVS (comma delineated) into a GDB to obtain marching OBJECT IDs. Then continued with the join which was successful.

Now my problem is when I tried to generate a population density map using quantities the entire layer will vanish from the workspace. Zooming to layer does not bring it back. I am not sure where I might have gone wrong. I have tried checking geometry for errors and had corrected the errors, I had matching coordinates. But still the same problem.


texts: QGIS 3.12 Toolbox Vector Creation - Create points layer from table corrupting field data

texts: I'm fairly new to this I have imported a table from Excel with coordinates and wanted to convert it to a points layer using the above tool. But it appears to corrupt fields that are real numbers for example 20.123 becomes 20.000, -39.12345 becomes -3.00000 and 174.123456 becomes 174.000000. Is this a known issue or am I missing something?


texts: QGIS 2.4: Basemap.at is displayed in poor quality

texts:
I have a huge problem displaying the austrian Basemap.at WMTS Service in QGIS 2.4. When I add the service, qgis logs multiple errors on failed tile requests. When I zoom into the map, it is displayed in a very poor quality (see screenshot). I really don't know how to fix this problem - can anyone tell me how to increase the quality of the tiles?WMTS URL: http://maps.wien.gv.at/basemap/1.0.0/WMTSCapabilities.xml

UPDATE 1:
I tried setting the resampling methods in QGIS to Cubic/Average, as you guys advised me to and compared the result with the same zoomlevel in ArcGIS. In this screenshot you can see the result: At the same scale level QGIS just skips an entire tile level, whereas arcgis loads this level correctly. - Any other ideas?


texts:
Following OpenLayers map printing demo here, I'm able to get my map in image format. But I'm facing two issues with the export:

Google Street base map is not included in the image output;
Whenever I activate any WMS layer, the export fails with the error TypeMismatchError: The type of an object is incompatible with the expected type of the parameter associated to the object

To get an image output of the map, I'm using Canvas.toDataURL(), which opens another window with the image of my WFS layers. According to the demo, when using Canvas.toDataURL(), images should be requested through a proxy, a point I don't understand as far as Google base layer is concerned. Any clarification will be appreciated.

I'm not sure what I may be doing wrong. Below shows how I retrieve Google base map and WMS layers:

Google base map:

var gStreet = new OpenLayers.Layer.Google(
            'Google Street',
            {numZoomLevels: 19},
            {visibility: true},
            {isBaseLayer: true},
            {useCanvas: OpenLayers.Layer.Grid.ONECANVASPERLAYER}
    );

WMS Layer:

var wmsLayer = new OpenLayers.Layer.WMS("My WMS Layer",
        GEOSERVER_HOST + "/wms",
        {layers: "wms_layer", format: 'image/png'},
        {
            useCanvas: OpenLayers.Layer.Grid.ONECANVASPERLAYER,
            isBaseLayer: false
        }
    );

GEOSERVER_HOST is proxied via Apache web server. I've configured Apache to channel all GeoServer requests as necessary. As such, I'm not making direct request to GeoServer in my client.

Any idea on what I may be doing wrong will be appreciated.


texts:
How can GeoJSON data to TopoJSON be converted?

TopoJSON files are often smaller than GeoJSON and the conversion allows you to incrementally simplify a dataset, if you wish.


texts: How can I delete polygon nodes in Carto Builder?

texts: Editing a polygon is possible by clicking on it, and the edit "pencil" icon appears, allowing you to drag or even add new coordinates. There is, however, no way to delete coordinates, as far as I can tell, except by doing an update function using SQL and manually finding and typing in the lat-long coordinates. That would be immensely tedious, since the coordinates' locations in my map are based on visual appearance rather than any quantitative relationship.


texts: Aggregate number of features by country in Overpass

texts:
Does Overpass (or Overpass Turbo) provide a mechanism for query like “show me the number of nodes tagged foo=bar by country”?

Output being something like this:

Country    Count
Estonia    5
Latvia     8
Lithuania  6
Poland     42683
Czechia    2653
Germany    1458

Doesn’t have to be nicely formatted, CSV would be fully sufficient.

Is that possible in Overpass? If so, how?


texts: Getting name of raster in point shp.file Attribute Table using ArcPy? AttributeError: 'unicode' object has no attribute 'contains'

texts:
I am trying to get name of raster for corresponding point if point belongs to raster    and fill column "name"  in   points Attribute Table:

but no names in Attribute Table column "name" of point's shp. file
What can be a problem?

I am getting   error

File "", line 22, in
    if iR.contains(pnt[0]):

AttributeError: 'unicode' object has no attribute 'contains'

env.workspace = r"C:\Users\1_fire_dNBR" # folder with set of rasters

inputPoints = r"C:\Users\select3.shp"

# List of rasters in specific folder
RExtents = [arcpy.Describe(R).extent for R in arcpy.ListRasters()]

#Names:
#RNames = [arcpy.Describe(R).name[:-9]]
RNames = [arcpy.Describe(R).name[:-10] for R in arcpy.ListRasters()]  #  extract  names of rasters
print RNames

with arcpy.da.UpdateCursor(inputPoints, ["SHAPE@","name"]) as iCur:

    for pnt in iCur:
        counter = 0
        for iR in RNames:
            if iR.contains(pnt[0]):

                print iR
                i.value.append(RNames[counter])
                iCur.updateRow["name"]
        counter =+ 1


texts: How can I make it so  that the gradient classifer (look closely, there are bits of red in some places) is more pronouced? Or that the outlines of the zip codes are less pronounced?


texts: How can I fetch Coverage Band Details from GeoServer in OpenLayers?

texts:
I want to display the gray index value in the image above in GeoServer on my map when clicked
At the bottom of the example there is (layer preview-> openlayers) I want something like this


texts: How to make a rectangle of an exact size/aspect ratio in QGIS?

texts:
Is there any way in QGIS to make a vector layer with a rectangle of an exact size or aspect ratio? Ideally, it would be nice if the rectangle from extent tool could be limited to a set aspect ratio, but barring that, is there any way to adjust the size of a vector rectangle directly?

[EDIT]
This is for QGIS 3.4


texts: Setting up if/then/else for choice of two options or boolean in ArcGIS ModelBuilder?

texts:
I'm using ArcGIS Desktop 10.6 and have a model as shown in the image below. I'm calling it option 1. The model exports GIS data to CAD DWG format and is comprised of two submodels, "Setup" and "Export".  I'm going to be setting up a second model (although I may house it within this model), option 2, also exporting GIS data to DWG but with different configuration options.

How can I configure ModelBuilder so when the model is run, it will offer a choice of option 1 or option 2 and then based on the choice, it will run the respective path?

I have no Python experience, therefore, am trying to do it in ModelBuilder. Python suggestions are welcome as I'm not sure if this is even possible to do in ModelBuilder. I was experimenting with the Boolean variable as you can see in the screenshot but am not sure how to configure it.


texts: Get layers and apply renderer to layer list from SDE ArcMap

texts:
I have the code that gets the layers from SDE, but when I add them to the map they have a different symbology. There are about 80 layers and I am accessing the SDE by getting the workspace from a selected layer in the TOC.

The symbology renderer layer files are located in a directory. So all I want to do is loop through the directory and get the renderer for each layer in my list, but my code only gets the renderer for one file not for all.  How can I compare feature class names so it will get the renderer for each layer by using the feature class name?

Public Sub AddSDECoverageLayer()

Dim pMxDoc As IMxDocument
Dim pMap As IMap

Set pMxDoc = ThisDocument

 Dim intX As Integer

'Add each FeatureLayer to the map document
'This should assign layer names to the layer it gets from SDE
'But it only gives the first layer in the list a layer name.
'There are about 80 layers so can't do if statement 80 times
'What other way can I do this?

For intX = 0 To ListBox2.ListCount - 1

       Dim xsClass As String
       xsClass = ListBox2.List(intX)
       Dim xs As String
       If xsClass = "PROPERTY_ADMIN.topo_ownership" Then
       xs = "MM Ownership"
       End If
       If xs = "PROPERTY_ADMIN.MM_Rights_of_Way" Then
       xs = "Rights of way"
       End If

      Dim pFLayerXS As IFeatureLayer

      Set pFLayerXS = importData(pFWorkspace, xsClass, xs)

      pMxDoc.FocusMap.AddLayer pFLayerXS

      pMxDoc.ActiveView.Refresh
      pMxDoc.UpdateContents

Next

End Sub

---------------------------------------------------------------------------

Public Function importData(pWorkspace As IWorkspace, FClass As String, fcName As      String) As IFeatureLayer

'Function to create Feature layer from Feature class
Dim pFeatureWorkspace As IFeatureWorkspace
Dim pFClass As IFeatureClass
Dim pUnknown As IUnknown
Dim pMxDoc As IMxDocument
Dim pFeatureLayer As IFeatureLayer
Dim pDataSet As IDataset
Dim pFLayer As IFeatureLayer

 ' Get the selected layer
Set pMxDoc = ThisDocument
Set pUnknown = pMxDoc.SelectedLayer 'pMxDoc.FocusMap.Layer(0)
 If Not TypeOf pUnknown Is IFeatureLayer Then
 Exit Function
End If
Set pFLayer = pUnknown

' Make sure the layer is from SDE and get the parent workspace
Set pDataSet = pFLayer.FeatureClass

Set pFWorkspace = pDataSet.workspace
Set pWorkspace = pFWorkspace
If Not pWorkspace.Type = esriRemoteDatabaseWorkspace Then
MsgBox "Not an SDE feature class"
Exit Function
End If
Set pDataSet = pWorkspace

Set pFeatureWorkspace = pWorkspace
Dim intY As Integer
For intY = 0 To ListBox2.ListCount - 1
Set pFClass = pFeatureWorkspace.OpenFeatureClass _
  (ListBox2.List(intY))
  Next

Set pFLayer = New FeatureLayer
Set pFLayer.FeatureClass = pFClass

Dim pGFLayer As IGeoFeatureLayer
Set pGFLayer = pFLayer

Dim intz As Integer
For intz = 0 To ListBox2.ListCount - 1

'This gets the renderer but only for one layer,
'How can i make it get the renderer for all the layers in the list.
'There is 80 layers, in a directory. i just need to assign the renderer
'from this directoryto eac item in the list box

If ListBox2.List(intz) = "PROPERTY_ADMIN.topo_ownership" Then
  Set pGFLayer.Renderer = GetRenderer("V:\Arc\Layers\MM Ownership.lyr")

  End If
  If ListBox2.List(intz) = "PROPERTY_ADMIN.MM_Rights_of_Way" Then
  Set pGFLayer.Renderer = GetRenderer("V:\Arc\Layers\MM Rights of Way.lyr")
End If
  Next

 pFLayer.Name = fcName

Set importData = pFLayer

Exit Function

End Function

I appreciate your time and any ideas.


texts: Dealing with a hole (or inner ring) at the intersection of two MultiPolygon polygons in Python Shapely


texts: How can I access the r.in.lidar algo with the pyqgis console? i tried to import processing and algolist but r.in.liad under the GRASS GIS algo does not get listed. Also Checked the options iunder toolbox processing to activate the GRASSGIS 7 commands


texts: Error using ArcPy SearchCursor to get list of unique values in multiple fields

texts:
I'm trying to extract a list of unique values in two fields in a shapefile. I think I am misunderstanding how the search cursor works, because I get the same error on which ever variable (x or y) I put second.
I'm following the answer here, and it works for just one field, but not 2+
Code:
with arcpy.da.SearchCursor(myshp, ['myid', 'myname']) as cursor:
    x = sorted({row[0] for row in cursor})
    y = sorted({row[1] for row in cursor})

Error:
y = sorted({row[1] for row in cursor})
IndexError: list index out of range


texts: saving survey123 data programmatically using python

texts: I'm recording data using the Survey123 app which logs data into a standard ESRI rest service, stored in my AGOL account. I would like to save this data an SDE dataset. Is it possible to do so using python (arcpy or other library)?


texts: I'm having a hard time printing a map within QGIS. For instance I did a 18x24 inch map printed on a 18x24 inch map. But the printed a map on a 9x12 inch printed on a 24x36 piece of paper. Why is it hard to print a big/wall map out in QGIS? Is there a step I'm missing?


texts: Error when running Split by Attributes tool from arcpy- arcgisscripting.ExecuteError: 'str' object has no attribute 'polygon' -

texts:
I am trying to use the Split by Attributes tool in arcpy to split a shapefile with 18 point features into 18 separate shapefiles.  It works as expected in ArcGIS Pro, but when I run the following script, I get the error: arcgisscripting.ExecuteError: 'str' object has no attribute 'polygon'

import arcpy
arcpy.env.overwriteOutput = True
inputfc = r'C:\Users\Ian\Documents\ArcGIS\GIS test data\SJemezNamedSites.shp'
fieldname = 'short_name'
target_workspace = r'C:\Users\Ian\Documents\ArcGIS\GIS test data\output'
arcpy.SplitByAttributes_analysis(inputfc, target_workspace, fieldname)

Additionally, if I run the tool from ArcGIS Pro, and then choose to "Copy Python Command", I get:

arcpy.analysis.SplitByAttributes(r"C:\Users\Ian\Documents\ArcGIS\GIS test data\SJemezNamedSites.shp",
                             r"C:\Users\Ian\Documents\ArcGIS\GIS test data\output", "short_name")

If I simply run that script, I get the same error as I do with my script.  What am I doing wrong?

I am using ArcGIS Pro 1.3.1, arcpy with python 3, and writing my script in pycharm IDE


texts: Use specific tile levels when cascading tile service as WMS in MapProxy

texts:
When cascading a tile service through MapProxy to provide a WMS service is it possible to restrict the zoom levels the source tiles are pulled from?
The original tile service I'm using uses different imagery sources at different zoom levels e.g. when zoomed out clouds and older imagery are shown, but as you zoom in to higher resolution imagery the clouds are lost and imagery is more recent.
Ideally I'd like to restrict the zoom levels used from the source tiles and have MapProxy stitch tiles together to provide the WMS service imagery at different resolutions, providing consistent imagery no matter the resolution viewed. Is this possible and are there any configuration examples I can use as reference?


texts: Finding repeating polygons in multiple geodatabases (.gdb)?

texts:
I'm working with a database of land parcels and their corresponding crops per year (source)

My end goal is to create a predictive model for crop type, given land parcel.
For example, if in the past 4 years a farmer grew soy, soy, corn and corn, what is the likely crop he'll be growing in 5th year.

To make the task a bit straightforward, I'm checking if the parcel polygons are the same (or similar) from year to year? If that's the case I should be able to abstract the geo info away and work directly with "parcel ID", crop type, and year data.

The gdb file for each year is quite sizeable and just loading and preprocessing in R studio takes quite a while.

What is a sane way to compare ploygon data across years and extract those exact parcels which appear in consecutive years? Would using the center of each polygon or labpt be a good way to represent the polygon?

I'm using Mac Sierra and rgdal library.
An entry from 2009 gdb after converting to GeoJson looks like this:

  { \"type\": \"FeatureCollection\", \"crs\": { \"type\": \"name\",
  \"properties\": { \"name\": \"urn:ogc:def:crs:OGC:1.3:CRS84\" } },
  \"features\": [ { \"type\": \"Feature\", \"id\": 2, \"properties\": {
  \"CAT_GEWASCATEGORIE\": \"Grasland\", \"GWS_GEWAS\": \"Grasland,
  blijvend\", \"GWS_GEWASCODE\": \"265\", \"Shape_Length\": 898.222510,
  \"Shape_Area\": 4885.337158 }, \"geometry\": { \"type\": \"Polygon\",
  \"coordinates\": [ [ [ 3.98569255494544, 51.400796253376114 ], [
  3.985792633157566, 51.400648243993814 ], [ 3.985454735105648, 51.400542942730567 ], [ 3.985419563702146, 51.400544767402884 ], [ 3.985367375532934, 51.400571777661298 ], [ 3.985239821317514, 51.400662957371289 ], [ 3.984954651582088, 51.400851050506539 ], [ 3.983832093498555, 51.401616824104352 ], [ 3.982416851757605, 51.402598958984925 ], [ 3.981780494933279, 51.403048420036697 ], [ 3.981523333282359, 51.403247946509204 ], [ 3.981396023...

This is my first geo project.

Edit

I noticed that quite a lot of polygons repeat across the years but with very slight variations in borders; so using labpt as is won't work.

2009 (red) and 2016 (green)

Only 2009 (red)


texts: Plotting polygon bounding box using ArcGIS API for Python

texts:
I have a layer of polygons and would like to be able to plot the extent boundary rectangle on a map. I have my polygons in a Spatially Enabled DataFrame (polygonSDF = poygonFeatureSet.sdf) that I can successfully plot using spatial.plot(). I've so far tried two different potential approaches to this problem and in both cases have run into errors I haven't been able to solve. I don't want to use draw() as described in the Spatially Enabled DataFrames Advanced Topics tutorial because then you lose the attribute information of the original polygon, and just have a static graphic.
Can anyone either point out a fix to one of the approaches described below, or suggest an alternate way to solve this problem? I'm working with ArcGIS API for Python 1.6.1 in the ArcGIS 10.7.1 Jupyter Notebook set-up on Windows 10.
I'm calculating and storing the extent data as follows:
extentList = []
for f in polygonFeatureSet.features:
     # Calculate the extent
     (w, s, e, n) = Geometry(f.geometry).geoextent
     extent = Polygon({
     "rings" : [[[e, s],[w, s],[w, n],[e, n], [e, s]]],
     "spatialReference" : {"wkid" : 4326} })

     # Add the extent to extentList
     extentList.append(extent)

# Add the extents to the dataframe as a new column
polygonSDF["MBRextent"] = extentList

Approach 1: Specify which column of geometry plot() uses
Issue: plot() doesn't appear to have this functionality, and just uses SHAPE automatically
Approach 2: Make a copy of my polygon SEDF, and then set the SHAPE column to store the extent data instead
Issue: Despite extentList containing valid Geometry objects, once loaded into the SEDF the MBRextent column is of type 'object', and thus is unable to be recognized and used as geometry.
# Attempt 1
copySDF['SHAPE'] = copySDF['MBRextent']

# Attempt 1 Error
AttributeError: Cannot use 'geom' accessor on objects of dtype 'object'.

# Attempt 2
copySDF.set_geometry(['MBRextent'])

# Attempt 2 Error
TypeError: Input geometry column must contain valid geometry objects.

# Attempt 3
extentDF = pd.DataFrame(polygonSDF["MBRextent"]).set_geometry(["MBRextent"])

# Attempt 3 Error
TypeError: Input geometry column must contain valid geometry objects.

Approach 2B: To try to rectify the issues caused by the extent being of type 'object' I've tried explicitly setting the column type to Geometry.
Issue: Geometry type is not being recognized by pandas methods
# Attempt 1 - replacement of Line 12 when creating the extent column
polygonSDF["MBRextent"] = pd.Series(extentList, dtype=Polygon)

# Attempt 1 error
TypeError: dtype '<class 'arcgis.geometry._types.Geometry'>' not understood

# Attempt 2 - setting type after creation of extent column
polygonSDF['MBRextent'].astype(Geometry)

# Attempt 2 error
TypeError: dtype '<class 'arcgis.geometry._types.Geometry'>' not understood

# Attempt 3 - setting type after creation of extent column
polygonSDF["MBRextent"].astype(polygonSDF['SHAPE'].dtype)

# Attempt 3 error
NotImplementedError:

Note: This same question has also been posted on GeoNet. Trying to get wider visibility here due to lack of answers so far.


texts: Is it possible to run the address locator in ArcGIS via Python?

texts:
I have a question regarding the address locator in ArcGIS.

It is possible to create an address locator by using Python, but is it also possible to run an address locator via Python?

This means, is there any Python code that does the same as the address locator toolbar in ArcGIS?

The code should get the address by an user and should then locate the address and zoom on it.


texts: Working with QGIS and Grass I would like to produce a ZTV. I have the bare land DEM happily ready to go but would like to add buildings and woodland to the calculation. Can anyone please point me in the direction of how to add buildings and woodland to the DEM upon which I will be doing the ZTV?


texts: Converting VBScript to JavaScript for labelling in ArcGIS Desktop?

texts:
I have a VB script used by ArcGIS 10.2 for Desktop to label different types of streets based on vales present in a specific field. I am a novice at Python but VB and JavaScript are beyond me. I need this to be JavaScript for older version of ArcGIS. Here is the script in VB:

Function FindLabel ( [STREETNAME] )
' Extracts number from [STREETNAME]
' ---------------------------
  Dim sIn, sNew, i
  sIn = [STREETNAME]  ' input field value
  sNew = ""  'result

  i = 1

  ' check if char is numeric
  While i <= Len(sIn)
    If IsNumeric(Mid(sIn, i, 1)) Then
      sNew = sNew & (Mid(sIn, i, 1))
    End If
    i = i + 1
  Wend

  FindLabel = sNew

End Function


texts: What does the predictor switch do while compressing geoTIFFs in GDAL

texts:
I am merging several geoTIFFs into one and am using the -co COMPRESS=DEFLATE create option for geoTIFFs. An additional switch is the predictor switch.

What does it do?
The GDAL website is not clear on this. It mentions horizontal differencing and floating point prediction. But I have no idea what is meant by these terms.
Thanks for your answer.


texts:
I am working on ArcGIS 10.2.1. I have a feature class containing more than 4000 records and a table having approximately 1000 records. Both feature class and table have some common fields however, with different values in them. Feature class is in ArcSDE geodatabase and table is in file geodatabase. We want to update feature class's fields with the values in table.

There are multiple fields (100 above). I was thinking of joining these two and then use field calculator but as records in feature class are more than the table and there are multiple fields which I want to update.

Is there any python script or tool I can use to update the feature class values from the table?


texts:
I have a whole bunch of GPX files that I need to load into QGIS and loading them one at a time using the GPS plugin is a bit painful.  It will not allow one to select multiple files -- is there another way?

I did find this:
Automating batch load of multiple GPX files into PostGIS
but I need to figure out a lot more about PostGIS and GDAL before trying that approach.

BTW I am using Mac version of QGIS -- just realised that this may be OS related


texts: Problem with same legend for two layers in Leaflet

texts:
In Leaflet, I have two layers (treated as overlays) and a single legend that I want to appear when any of them is checked in the Layer Control panel. If they are both checked, I don't want it to appear twice.
This is the code:

map.on('overlayadd', function(eventLayer) {
  if (eventLayer.name === 'Drumuri județene' || eventLayer.name === 'Drumuri comunale') {
    legend1.addTo(map);
  } else {
    if (others < 1) {
      legend1.addTo(map);
    }
    others++;
  }
});

map.on('overlayremove', function(eventLayer) {
  if (eventLayer.name === 'Drumuri județene' || eventLayer.name === 'Drumuri comunale') {
    this.removeControl(legend1);
  } else {
    if (others < 1) {
      this.removeControl(legend1);
    }
  }
});

The problem with it is that it works only if they are checked and unchecked separately. If they are both visible, two identical legends appear and If I uncheck them one after the other, one of the legends fails to dissapear. I can tell there is something wrong, but I don't have the programming knowledge to solve it. Can someone offer a solution?


texts: GRIB data format investigating octets of sections with Python

texts:
I'm trying to learn GRIB data format and I found this helpful guide to dive into
https://www.wmo.int/pages/prog/www/WDM/Guides/Guide-binary-2.html#Section3.
I had an experience with Gdal and Pygrib. Yet, I wanted to dive into for deeper knowledge on this data format.
How can I learn GRIB data format more deeply? There is not enough resources to learn it and I' intermediate with the programming


texts: How to generate routes that follow underlying road network?

texts:
I have a set of points (~2300 bus stops) and I have a network of roads that these bus stops lie on. When I have a bus route that goes from A to B, how do I generate the route such that it follows the underlying road and not just the direct crowfly distance.

The road data is from OSM and the stop data is user generated.
At present the data is hosted on :routemaster.hydlab.in