texts: How to detect mountain ridges and add them to a NO-GO area

texts:
I want to calculate a simpliefied mobility model that classifies an terrain map into movement areas. The Mountains shall be like walls and everything else shall be movable terrain. I started to classify areas with a slope of 0-30% as GO areas and areas > 31% as NO-GO areas. As input I have a slope raster file. However, when I do a raster calculation for each cell being 1 when there is a slope of 0-30% and 0 otherwise, I still have some white lines which seem to be ridges of the mountains. How can I include these lines in the NO-GO area? Or is there another way to reach my goal? I am using QGis.

To illustrate my problem:

This is my slope raster input file:

This is my result when doing the raster calculation step:

As you can see in the second picture, there are white lines indicating mountain ridges. How do I add them to the NO-GO areas?


texts: Why does Google satellite layer change the project CRS?

texts:
I have a project and a topomap in WGS84 (EPSG4326).
On the fly CRS transformation is activated.
I open Google satellite layer from the OpenLayers plugin.
The Google layer is WGS84 / Pseudo Mercator (EPSG3857).
But then the whole project has also become WGS84 / Pseudo Mercator (EPSG3857).

I thought on the fly transformation would reproject to EPSG4326. Why is that? Is it better to work directly in WGS84 / Pseudo Mercator (EPSG3857) systematically when we use Google Earth?


texts: Checking original geometry when using ST_Split function with MultiPolygon and LineString in PostGIS?


texts: Dynamic triggers with table_name in PostGIS

texts:
I want make a dynamic trigger with PostgreSQL/PostGIS like this. I try that triggers make more generic for use with many tables with the same pattern. When i use direct table's name in code, that's works but

CREATE OR REPLACE FUNCTION asst.inter_geom_insee() RETURNS TRIGGER AS
$body$
    DECLARE
        tbl text := quote_ident(TG_TABLE_SCHEMA) || '.' || quote_ident(TG_TABLE_NAME);
    BEGIN
        IF TG_OP = 'INSERT' or TG_OP = 'UPDATE' THEN
            UPDATE '|| tbl ||' SET insee = c.insee_com FROM cadastre.COMMUNE AS c WHERE st_within('|| tbl ||'.geom, c.geom);
        END IF;
    RETURN NEW;
    END;
$body$
LANGUAGE plpgsql;

I try in my code

EXECUTE 'UPDATE '|| tbl ||' SET insee = c.insee_com FROM cadastre.COMMUNE AS c WHERE st_within('|| tbl ||'.geom, c.geom)' USING NEW;

but nothing, just one answer from dbeaver :

  ERROR: asst.inter_geom_insee() is not a procedure

I see https://stackoverflow.com/questions/7519044/how-to-dynamically-use-tg-table-name-in-postgresql-8-2
or https://stackoverflow.com/questions/1997337/inserting-new-from-a-generic-trigger-using-execute-in-pl-pgsql
But I fail.


texts:
I am using qgis to create polygons in a postgis layer.  I can create polygons and save them without a problem, but I am unable to delete them or use the node tool to reposition the points of existing polygons.  Editing is enabled, but both the node tool and delete selected tool are greyed out (even when a feature in the current layer is selected):

I have seen posts elsewhere suggesting that these sorts of edits should be possible.  Is this a bug in the current release of qgis, or are these features simply not supported?

Update

It turns out I can't edit any attributes in the attribute table, not even when using the "add feature" button in the attribute table.


texts: OpenLayers 3.7 New Server Vector WFS with GeoServer

texts:
I have recently updated from OpenLayers 3.0 to 3.7 and I noticed in the 3.5 update the Vector API has changed, I tried to alter the code as per the bottom section of code but i am getting the error below:

XMLHttpRequest cannot load http://PCNAME:8080/geoserver/example/ows?service=WFS&version=1.0.0&re…-180803.92191974312,7068079.039279004,-170676.01567195728,7073596.83725834. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://localhost:9000' is therefore not allowed access.

Code as per below:

var vectorFeatures = function(response) {
    var features = vectorSource.readFeatures(response);
    vectorSource.addFeatures(features);
  };

  var geoJSONFormat = new ol.format.GeoJSON();

    var vectorSource = new ol.source.Vector({
      loader: function(extent, resolution, projection) {
        var url = 'http://PCNAME:8080/geoserver/Example/ows?service=WFS&' +
        'version=1.0.0&request=GetFeature&typeName=Example:layer_1&' +
        'outputFormat=text/javascript&format_options=callback:vectorFeatures' +
        '&srsname=EPSG:3857&bbox=' + extent.join(',');
        $.ajax({
          url: url,
          success: function(data) {
            var features = geoJSONFormat.readFeatures(data);
            vectorSource.addFeatures(features);
          }
        });
      },
      strategy: ol.loadingstrategy.bbox
    });


texts: Code not working after reprojecting raster (QGIS)

texts:
I'm running this code in two different projects (Porject1 contains a original raster with original CRS; Project2 contains the same raster reprojected in a customized CRS ):
poly = iface.activeLayer()
a = processing.run("qgis:fieldcalculator",{'INPUT': poly, 'FIELD_NAME': '1.3-field', 'FIELD_TYPE': 0, 'FIELD_LENGTH': 10, 'FIELD_PRECISION': 3,'NEW_FIELD': True,'FORMULA': '((length(shortest_line(centroid($geometry),boundary($geometry)))))','OUTPUT': 'memory:'})

When I run it in project 1  the code runs correctly.
But when I try to run it in project 2, it gives this error:
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.6\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 2, in <module>
  File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\tools\general.py", line 105, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\core\Processing.py", line 183, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.

The only thing I could think that could be causing this error is something related to the custom (reprojected) CRS.
The original CRS created when I add the raster to the project is:
"+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +a=1188300 +b=1188300 +units=m +no_defs"

The curtomized CRS that I am using to reproject the raster using warp is:
"+proj=aeqd +lat_0=7 +lon_0=109 +x_0=0 +y_0=0 +a=1188300 +b=1188300 +units=m +no_defs"

I need to run the code with the reprojected raster and custom CRS. Any ideas?


texts: I have a variety of raster files from the exact same study area including different values (hourly mean radiant temperature). Is it possible to find out the overall min/max values (not each individual one)?


texts: Installing Shapely with libgeos_c on Linux?

texts:
The installation instructions for shapely suggest:

"Windows users should use the executable installer, which contains the required GEOS DLL. Other users should acquire libgeos_c by any means, make sure that it is on the system library path, and install from the Python package index."

What exactly is libgeos_c? Is it the same as GEOS - Geometry Engine, Open Source?


texts: All coordinates with an specific distance to one point

texts:
I'm sorry for my english!
I have to get the following informations:
I have a specific point with long/lat.
Then i draw a circle with a specific radius (e.g. 5km).
The intersections of all roads on the line of the circle must be displayed as a point.

My idea with photoshop:

I'm not a real coder, but a mathematician...

So I had the idea to calculate all coordinates, which have the specific radius (e.g. 5km) on the first step. The second step is to check, if the point is on a road. If the point is on a road, I save it in an array. If the point isn't on a road, I delete the point for the session. The third step is to display all points of the array on a map.

That first step works. The distance between the points on the line of the circle are about 0,5m. That is an acceptable tolerance. I get a list of coordinates with PHP. E.g. with 1km radius:

B 32.0179859432° L 20.999880002°
B 32.0179846637° L 20.9997800137°
B 32.0179825847° L 20.9996800441°
B 32.0179797064° L 20.9995801021°
B 32.0179760289° L 20.9994801966°
B 32.0179715527° L 20.9993803363°
B 32.0179662781° L 20.9992805303°
B 32.0179602056° L 20.9991807874°
B 32.0179533358° L 20.9990811164°
...

But that are to many results depending on the distance.
1km: 6663 points
10km: 66630 points
...

That are to many points. I can't use a reverse geocoding service to check, if every point is on a road. That would be about 60000 requests in one query. I believe, that only 0,1 % of the points are on a road, so my solution is not very efficient.

Do you habe any idea to solve my problem??
Maybe with frameworks like openlayers, or ...?
Or API's, which can display all coordinates with an specific distance to one point?


texts: Why can't I add to ee lists with List.add?

texts:
I was trying to create a list with the following code to then classify:
  var randomStratClassifierIterator = function(times){

  var loops = 0;
  var classified = ee.List([]);

  while(loops < times){

    //Classify image and add to a list
    classified.add(ee.Image(0)));

    loops += 1;
  }

  return ee.ImageCollection.fromImages(classified).reduce(ee.Reducer.median(2));

};

But it won't add to the list. I then tested this outside of the function and I can't use List.add() to append anything at all. I switched to using a non-EE list and using .push(), which works fine but I'm wondering if I'm using .add() wrong? It seems like it should be straightforward and the docs show it should support ee objects, strings, etc (none of which worked for me).


texts: Postgres Performance with Radius Searches

texts: Im porting an application from mysql to postgres, and have several tables with 500k+ rows in each of geo points (lat, lon). In mysql i used a bounding box query with haversine to select rows within a given distance of a center point. i understand i can use either postgis's st_dwithin or the earthdistance module (or i could probably port my mysql query over too). Just wondering if anyone has a recommendation on which method would be best performance wise with these large data sets.


texts: Using arcpy to turn on 'label features' within model builder - why wont the script tool accept input data?

texts:
I'm relatively new to Python coding. I'd like to write a script tool for turning 'label features' on within an ArcGIS model. However, the script I've written won't recieve input data.

Other example scripts for labelling in this forum have specified layer file paths or used 'listLayers' within an mxd to source the input. However, I need this to take an input connection within model builder (so it can be re-used with different data without modifying the code). So I'm using 'GetParameterAsText'.

Code:

import arcpy

inputlayer = arcpy.GetParameterAsText(0)
label_field = arcpy.ListFields(inputlayer)[3]
lyr = arcpy.mapping.Layer(inputlayer)

lyr.supports("LABELCLASSES")
lyr.labelClasses[0].expression = "Abt  + label_field"
lyr.showLabels = True
arcpy.RefreshActiveView()

I've been attempting to use this on polygon features and also lines. Specifically, the problems are:

When I run it as a standalone script tool, I cannot drag a layer (feature classes) into the input from the TOC
When I put the script into the model builder, it refuses to connect feature classes or feature layers as inputs.

Is there a better option than GetParameterAsText?

Or is the code incorrect or incomplete?

Is it perhaps an issue with how I've set up the script tool within ArcMap?


texts: Assigning QgsCategorizedSymbolRenderer default ramp colors?

texts:
I have a plugin and I want to symbology a vector file with some default ramp colors but I see that in graduated symbology is not so easy as in graduated symbology.

Is it possible to assign a default ramp color to a vector file?

My code is very confusing because I tried so many examples but I couldn't get the result.

for code, (color, label) in attribute_sym.items():
    sym = QgsSymbolV2.defaultSymbol(layer.geometryType())
    sym.setColor(QColor(color))
    category = QgsRendererCategoryV2(code, sym, label)
    categories.append(category)

self.field = self.dockwidget.comboBox.currentText()
renderer = QgsCategorizedSymbolRendererV2(self.field, categories)
layer.setRendererV2(renderer)
QgsMapLayerRegistry.instance().addMapLayer(layer)

This code works but I don't want to allways define the list with the colors, because I am trying to symbology several fields. Is it some way of using the default ramp colors?


texts: GEE clip to the borders of European Union + UK + countries inbetween

texts:
Is there any easy way to crop the results to the European Union (+UK +countries inbetween) borders?
For now I have ended with GAUL database, but it doesn't seem to be accurate at all - the borders do not match the ones on the Google Map - like sometimes totally. Does anybody know some other more accurate way?
Here's what I got for now, and not happy with the results:
function clip_by_countries(img, country_list){
  var countries_of_interest = GAUL.filter(ee.Filter.inList('ADM0_NAME', country_list));
  return img.clip(countries_of_interest);
}

var countries = ['Norway','Sweden','Finland','Estonia',
    'Latvia','Lithuania','Denmark','Ireland','U.K. of Great Britain and Northern Ireland','Netherlands',
    'Germany','Poland','France','Belgium','Luxembourg','Liechtenstein','Czech Republic',
    'Switzerland', 'Austria','Slovakia','Hungary','Romania','Bulgaria','Serbia','Italy','Grece',
    'Croatia', 'Slovenia', 'Monaco', 'Andorra', 'Spain', 'Portugal', 'Gibraltar', 'Bosnia and Herzegovina',
    'Greece','Albania','The former Yugoslav Republic of Macedonia','Montenegro'];

var mls = clip_by_countries(s, countries);

Is there a database that matches the Google Map borders? Or maybe you got some other solution for that problem?
Example:


texts: Importing dxf as a vector layer or using dxf2shp function (Differences?)

texts:
I was trying to import a DXF file in QGIS. First tried using dxf2shp function but it didn´t work. I browse the dxf file and i put a name to the output and nothing was created in directory i chose.  I couldn´t convert a DFX file in order to import it with this method, i dont know what i did wrong.

I tried in other way with "add vector layer" selecting (Source: file,code:UTF-8; and Browse my dxf) and it works.

what is the difference between the 2 ways, is there a disadvantage with the vector layer option?

Do you know why first method didn´t work for me?

Regards


texts: How to run OSgo4W shell script on Windows?

texts:
I have a lot of DEMs that I want to use to create hillshades, slopes and color relief from and would like to automate the processes using a shell script within OSgeo4W shell. Here is the code I used within the shell script to automate the creation of hillshades.

for f in *.tif
do
  echo "Processing $f"
gdaldem hillshade $f $f-hillshade" -z 1.0 -s 1.0 -az 315.0 -alt 45.0 -compute_edges -of GTiff
done

I'm new to OSgeo4W and shell scripting. I'm working off this example as a starting point. When I call the shell script within the OSgeo4W I get an error message that Windows can't open the file and then it gives me an option to use the web to find the correct program or to manually select a program to use.

I had assumed that OSgeo4W installed everything I would need to run shell scripts in Windows. What am I missing?


texts:
I am new with dotspatial and I want to enable and use plugins on my own. I have a custom button that I want to add, it's a click handle measure function and shows the results in a text label.

Is there a way to use the measure plugin in this case? How can I use it ?

I just need some hints and if it is possible a sample.


texts: How to interpret GRASS v.kernel results?

texts:
GRASS manual reads:

  v.kernel - Generates a raster density map from vector points data using a moving 2D isotropic Gaussian kernel ...

Ok, but how do I interpret the results? I understand that v.kernel is the more advanced than v.neighbor function but I'm unsure as to which advantages it has.


texts: I am using QGIS. I want to create perpendicular Lines on another Line (River). The Lines should be created on points on the Line and should have a length of 40 Meters. I want to create the red lines in the picture below (for the complete River).


texts: Hibernate Spatial - Query within an X kilometer radius?

texts:
I am new to Hibernate Spatial, and am trying to perform a simple query of objects within a given radius.  I've created a number of entries in my database with properties corresponding to a latitude and longitude, using data from Google Maps and other sources.  This property is defined like this in my Entity class:

@Column
@Type(type = "org.hibernate.spatial.GeometryType")
private Point coordinates = null;

I'm now trying to figure out how to do a search of all entity objects that have coordinates that fall within a radius of x kilometers from a given point.  For example, I'd like to find objects that fall within a 50 kilometer radius of the point (12.34567, -76.54321).  However, I can't find any examples or tutorials that would explain how to do this in Hibernate Spatial.

Can anyone give me any information on how a query like this can be constructed?


texts: Cause of ImportError: No module named mapscript in QGIS mapserver export?

texts:
I'm trying to expor project to mapfile in new QGIS Chugiak but a get this error below, someone have the same error? How can I fix this?thnks.

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/Artur/.qgis2/python/plugins\rt_mapserver_exporter\plugin.py", line 65, in run
    from .mapfileexportdlg import MapfileExportDlg
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/Artur/.qgis2/python/plugins\rt_mapserver_exporter\mapfileexportdlg.py", line 30, in
    import mapscript
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python\qgis\utils.py", line 454, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named mapscript

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.4.0-Chugiak Chugiak, 8fdd08a

Python path: ['C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\\processing', 'C:/PROGRA~1/QGISCH~1/apps/qgis/./python', u'C:/Users/Artur/.qgis2/python', u'C:/Users/Artur/.qgis2/python/plugins', 'C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins', 'C:\\INTERG~1\\ERDASI~1\\usr\\lib\\Win32Release\\python', 'C:\\PROGRA~1\\QGISCH~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISCH~1\\bin', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\qgis\\python\\plugins\\fTools\\tools', 'C:/Users/Artur/Documents/GIS DataBase']


texts: How to treat PosgreSQL array fields in QGIS rule-based styles filters

texts:
In a PostgreSQL/PostGIS DB I have a table with an attribute declared as a 1-dimensional integer array INTEGER[] holding a value like {1,2,3,4,5}. Now I would like to build a rule-based style and need to test if a certain value is in the array. I'm not quite sure if I have to treat it as text and test for substring like (let's say test for 3 is in the array)

`strpos("myattribute", '{3,')) > 0 or strpos("myattribute", ',3,')) > 0 or strpos("myattribute", ',3}')) > 0`

or if it is treated as some sort of list and I may use

`3 in "myattribute"`

giving an invalid expression. Can it be converted to a list in some way or is my text approach recommended?


texts: I am trying to compile GDAL with support to convert to PDF format. When I try gdalinfo --formats the option for PDF is not present. I am using 2.1.2. The command I am using to build GDAL is nmake /f makefile.vc WIN64=YES. Is it possible to convert raster data to PDF format using GDAL?


texts: Convert Postscript (*.ps), Encapsulated Postscript (*.eps), or Adobe Illustrator (*.ai) files to shapefile?

texts:
Are there any programs which can convert from Postscript, Encapsulated Postscript, or Adobe Illustrator file formats to a shapefile? What about adding geo-referencing?

I'm also open to recipes using a number programs if there isn't a single one to fit the bill.


texts: At what stage is it best to project your data? Or does it even matter?

texts:
I am doing a series of geoprocessing steps in ArcGIS with a couple different datasets, for a global analysis of habitat loss and impacts on biodiversity. I have raster forest cover tiles and polygon range data (unprojected and/or in the GCS). I have my series of geoprocessing steps figured out but I am wondering if it is better to project the data early on in the process or later on. For the rasters I figured out that I couldn't project the tiles prior to mosaicing them - otherwise there were gaps in the mosaiced output.

Is there a standard or best practice or does it even matter?


texts: Merging polygons in shapefile based on field values?

texts:
I have a shapefile with a lot of polygons in it.

These polygons are grouped into a category: district.

How can I merge the polygons based on the same name of the category within this shapefile?


texts: QGIS just installed on MAC OSX - vector layer styles do not work properly

texts:
I am a newbie with QGIS - I have just installed QGIS 2.10.1-Pisa. I am following a youtube tutorial by Steven Bernard and I just found that the styles do not work properly (I run Mac OS X Lion 7.5.1. on a 1.7 GHz Intel Core i5 MacBook Air).

When I want to change the sybology of my map, this does not work properly. If you have a look at step 3.2.3. Changing Symbol Structure in the manual on 3.2.1. basic Follow Along: Changing Colors , QGIS does not provide me with any choices. It tells me that 'This layer doesn't have any editable properties'. I am attaching a screenshot just to be clear.

Can anybody tell me what is wrong?


texts: Export google satallite raster in QGIS 2.12- problem

texts:
When I import aerial raster from Google satellite, using the OpenLayers plugin

the raster is located in the right place.
But when I export it the raster isn't located in the right place:

I read Problems exporting aerial imagery from QGIS and didn't find any suitable answer for my issue.


texts: How can I aggregate a raster in R using a quantile function?

texts:
I'm trying to create a new, smaller NDVI raster with the upper quantile of raster values for the new aggregated cells.

library(raster)
rast <- raster("ndvi.tif")

Now I'd like to create the new raster downsampled by a factor of 10, using a quantile function for the top 25% of the data:

rast2 <- aggregate(x = rast, fact = 10,
    fun = function(i) quantile(i, probs = 0.75, na.rm = T))

This always gives an error:

Error in FUN(newX[, i], ...) : unused argument (na.rm = TRUE)
Called from: FUN(newX[, i], ...)

If I use mean or max it works just fine, only quantile gives an error. I've also tried defining a custom function but with the same result.


texts: Getting NameError with PythEveryone in ArcGIS virtual campus module?

texts:
I'm working on the python for everyone lab in Esri's Virtual campus. I'm running into an issue that I can't seem to solve, I've already done some outside research hence posting here!

The issue:

4th Lab "Use python in ArcGIS" Step 5, when I run the code I'm getting this error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 3, in <module>
NameError: name 'arcpy' is not defined

My code seems to looks fine, I've doubled checked multiple times this is what it looks like:

>>> # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
... # The following inputs are layers or table views: "RI_Schools", "Providence"
... arcpy.Clip_analysis("RI_Sewers","Providence","C:/Users/Laura.Ye/Documents/GIS & RS Labs/GIS II Project Folder/Assignment 1/PythEveryone10_1/PythonInArcGIS/sewers.shp","#")
...
Runtime error
Traceback (most recent call last):
  File "<string>", line 3, in <module>
NameError: name 'arcpy' is not defined
>>>

I added an import arcpy statement as my first line of script but didn't get to far because I received the following:

      import arcpy
      Runtime error
      Traceback (most recent call last):
        File "", line 1, in
        File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy__init__.py", line 24, in
          from arcpy.toolbox import *
        File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\toolbox.py", line 351, in
          from management import Graph, GraphTemplate
        File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 22, in
          import _management
        File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy_management.py", line 14, in
          import _graph
        File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy_graph.py", line 27, in
          import numpy
      ImportError: No module named numpy

** This import error occurred earlier except that time it said "No module named arcpy"

Overall, I've been having a number of errors and issues, and it seems like I should remove/reinstall ArcGIS and Python. Both seem to be located in two different locations on my computer. One location is the C:\ and the other is C:\users. I'm looking in the property files for both, the first has 6,145 files and 33 folders while the second has 63 files and 3 folders.

A little more background:

Within the first folder location C:(labeled ArcGIS10.2) there's a Python folder and a bunch of other folders. Also, there is a separate folder not nested with the ArcGIS10.2 folder, but alongside labeled Python27, which contains 5,477files and 249folders.

Within the second folder location, C:\users\ (labeled ArcGIS) there's a pywin32-218.win-amd64-py2.7 application file.

My conclusion:

This all seems like an unorganized mess, so I'm tempted to delete everything and start over. I'm new to using Python but have used ArcGIS over the last year and a half so when I first installed everything I didn't pay much attention to Python. About 2 months ago I slowly started to begin learning Python, outside of the ArcGIS environment, and downloaded a version of Python2.7. Now that I have a little more experience with GIS etc. I'm much better at trouble shooting and figuring things out, but this seems like a crazy mess of unorganized and broken files/programs etc.!

I am due for an upgrade to ArcGIS 10.3 (with my student license) so perhaps starting fresh would be easier than trying to fix everything. At least then I would have a better idea if what I encounter is user error or something else.

Any thoughts?


texts: Exporting full Landsat image to GeoTIFF in Google Earth Engine?

texts:
I am trying to get a complete Landsat image in Google Earth Engine, but the export shows some mistakes like this:

I have calculated the total pixel numbers and it is impossible to have such a huge amount of pixels: 7621X7761X4=236586324

The link of code is here:https://code.earthengine.google.com/38c19f74ee68ff44a0c2c22c0231cb87

I am wondering is there any wrong with this export code? And is there any way to output GeoTIFF image in patches like we export images in tfrecord?

Update:

Although I tried to set the maxPiexels to successfully export a Landsat image according to the suggestion from Sean, it still has some problems:
First, the exported data are automatically separated into some files and I do not know how:

Second, when I downloaded all these images and open them in ENVI, the total size of these patches is not correct, and all values are 0.


texts: Why is Union in ArcMap much faster than other approaches?

texts:
Taking the union* of these two large shapefiles file1 file2 locally using ArcMap takes appr. 10 min.

Over the last few weeks, I've tried to obtain the same results using other, open-source approaches without success. I am curious why only ArcMap seems to generate results in a reasonable amount of time without errors.

The approaches I've tried so far:

QGIS

The most straight-forward alternative. There are two "union" functions available in QGIS. One default (part of vector overlay), one based on SAGA. Both functions create the desired results on smaller samples but fail to run (or take forever) on my shapefiles. I tried running union for 30 minutes after which the progress bar was at appr. 4%

PostGIS

I've ingested the shapefiles on an AWS RDS instance running postGreSQL with PostGIS enabled. The approach is not so straight forward and requires a hideous long SQL query to obtain the same result. Most importantly, the query did not create a meaningful result on the full datasets. It did run properly on a sample though.

See separate question here

CREATE TABLE test.sqlislelijk AS
-- input data
with polys1 AS (
  SELECT
    pfaf_id as df1,
    geom as g
  FROM hybas06_v04
),
polys2 AS (
  SELECT
    aqid as df2,
    geom as g
  FROM y2018m11d14_rh_whymap_to_rds_v01_v01
),
-- intersections
intersections AS (
  SELECT df1, df2, ST_INTERSECTION(a.g, b.g) i, a.g AS g1, b.g AS g2
  FROM polys1 a, polys2 b WHERE ST_INTERSECTS(a.g, b.g)
),
-- per-row union of intersections with this row
diff1 AS (
  SELECT df1, ST_UNION(i) i FROM intersections GROUP BY df1
),
diff2 AS (
  SELECT df2, ST_UNION(i) i FROM intersections GROUP BY df2
),
-- various combinations of intersections
pairs AS (
  SELECT df1, df2, i AS g FROM intersections
  UNION ALL
  SELECT
    p.df1,
    NULL,
    CASE
      WHEN i IS NULL THEN g
      ELSE ST_DIFFERENCE(g, i)
    END
  FROM polys1 p LEFT JOIN diff1 d ON p.df1 = d.df1
  UNION ALL
  SELECT
    NULL,
    p.df2,
    CASE
      WHEN i IS NULL THEN g
      ELSE ST_DIFFERENCE(g, i)
    END
  FROM polys2 p LEFT JOIN diff2 d ON p.df2 = d.df2
)
SELECT
  df1 as pfaf_id,
  df2 as aqid,
  g as geom
FROM pairs WHERE NOT ST_IsEmpty(g);

Google BigQuery

Google's GIS extension for BQ looks promising but I experienced weird behavior and internal errors, probably due to mixing of GEOMETRY and GEOGRAPHY objects. Also, the syntax is similar to PostGIS and a command of 1 line in geopandas takes >20 lines in SQL.

GeoPandas

Probably my favorite approach. The out-of the box approach does't give any results after running it for several hours. My latest approach includes tiling the shapefiles, use multiprocessing on each tile and merge them. Not very straight forward and one 10x10 degree tile till takes 40+ minutes to calculate.

See Overlay Union Geopandas improve performance

It is the first time in my career that I am impressed by the performance of ArcMap. It seems that ArcMap handles this use-case much, much better than anything else. The issue seems persistent over all solutions I've tried so far. Is there something wrong in the underlying libraries of the open-source/alternative solutions?

Note that Union in ArcMap means something different than in the Shapely / PostGIS world.


texts: How to convert debug.js to .js file after editing debug.js file

texts: How to convert a library file debug.js to .js file after editing the debug.js file. I edited opensource library file olgm.debug.js . Now i want to convert it into normal olgm.js file to use in application. Can any one help me to solve this problem. olgm is not supporting setStyle and setRotation etc.. For that i created new function edited in olgm.debug.js. Now how to convert to normal olgm.js.


texts: Calculating percentage (or length) of stream network that is multiple channel streams using ArcGIS Desktop?

texts:
I have a polyline stream network of all channels within a drainage basin (which is 1.1 million square km)
A number of the rivers in the basin alternate between single thread and multichannel (i.e. anabranching or braiding).

Is there a way using argis for me to calculate the length or percentage of total which are either multi or single thread, without simply measuring the along line distances?

The ideal end product would be two values:

percentage (or total length) multichannel
percentage (or total length) single channel


texts: How to create and use a list of unique dates of images in an ImageCollection GEE

texts:
I am looking for a way to get a list of all dates where SENTINEL2 created an image in a specific ROI.
Due to the fact that my ROI may contain more than 1 image per date (several tiles) - I get a list with duplicated dates.
I have reduced the list by frequencyHistogram, but I need to use the list for a function that iterates on that list and extracts all the images from these dates. (Later on, I would like to extract SEN2 image from every date in the list and 5 SEN3 images from that specific date and 5 days ahead).

I have used the code described here: Reference topic but I cannot use the unique list because it contains the repetition of each date and my ee.Filter.inList is unable to cope with that.

So, for now I use the un-reduced list, but I get duplicated images.

This is my code Link:

//Dates of Interest
var start = ee.Date("2014-10-01");
var finish = ee.Date("2020-01-01");

///--------------------- SENTINEL2 Collection ---------------------------------------///
var sen2 = ee.ImageCollection('COPERNICUS/S2_SR')
.filterDate(start, finish)
.filterBounds(roi)
.map(function(image){
  return image
  .clipToCollection(roi)
  .updateMask(image.select('QA60').not())
})

// Year-Month-Day Extract function
function ymdList(imgcol){
    var iter_func = function(image, newlist){
        var date = ee.Number.parse(image.date().millis());
        newlist = ee.List(newlist);
        return ee.List(newlist.add(date).sort())
    };
    return imgcol.iterate(iter_func, ee.List([]));
}
// MilliSec Extract function
function millisec(imgcol){
    var iter_func = function(image, newlist){
        var date = ee.Number.parse(image.date().millis());
        newlist = ee.List(newlist);
        return ee.List(newlist.add(date).sort())
    };
    return imgcol.iterate(iter_func, ee.List([]));
}

var ymd = millisec(sen2)
var ymd_Hist = ee.List(ymd).reduce(ee.Reducer.frequencyHistogram())

print("With histogram: ", ymd_Hist) //Need to figure out how to use this unique dates list
print("No histogram: ", ymd)

var allImagesInRange = sen2.filter(ee.Filter.inList("system:time_start", ymd))

print(allImagesInRange)

My questions:

Is there a way to reduce the list by key? So I will get only the unique dates?
How can I build another list that contains all the dates from the unique dates list and add (in addition to each one of the dates written there) 5 more dates? Let's say in time T I have a SEN2 image, I would like to add to the list T+1, T+2, T+3, T+4.


texts: QGIS: Splitting layer into individual shapefiles with the attribute as the shapefile name

texts:
In QGIS 3.14 I have created a layer containing multiple polygons from x/y coordinate data, but I now need to split those ploygons into individual layers/shapefiles by their different attributes - However, I need these new layers/ shapefiles to be named after the different attributes that they are split, but I only seem to be able to get layers numbered arbitarily 1 to x
Eg:
Attrbiute Field 4 has different feature numbers - eg. F021, F032, F063, which represent the individual polygons that I need seperated into different layers.
When I split by attribute to create the shapefiles I get field_4_1, field_4_2 etc. (where the _1 and _2 are arbitary numbers and not realated to the feature), instead of field_4_F021, field_4_F032 - or ideally just F021, F032 etc.


texts: Transfer data from SQL database to postgre database using GDAL

texts:
I am new to GDAL and have been uploading CSV crime data to a postgis database using ogr2ogr.

ogr2ogr -f "PostgreSQL" -lco GEOMETRY_NAME=geometry -lco FID=gid PG:"host=server user=user dbname=dbname password=****" C:\Users\TommyH\Desktop\csv_upload_postgres.vrt -nln th_crime_2017 -t_srs EPSG:27700

I am now wondering if there is a way to transfer data from postgis to SQL or Oracle using GDAL ?


texts: Removing "Bow-tie deletion line " from my raster dataset

texts: I am working on VIIRS SDR dataset where I have created geolocated and projected raster of "BrightnessTemperature" using "Latitude" and "Longitude" dataset given in HDF5 format using ENVI. Now my study area comes under bow-tie deletion line as uploaded in image below.


texts: Connecting Street lines to Speed Profiles for Network Analyst

texts:
I am creating a network dataset with speed profiles. I have a unique connection between the two where Link_PVID and OLD_FEAT_ID have a unique connection between the two. But it seems that the Patterns needs to have the OBJECTID from the streets as per https://desktop.arcgis.com/en/arcmap/10.3/guide-books/extensions/network-analyst/copying-source-feature-classes.htm

I get the following error for all streets

SourceName: Speed_Patterns, ObjectID: 102184814, Network edge not found with positions 0.000000 to 1.000000
SourceName: Speed_Patterns, ObjectID: 102184814, Network edge not found with positions 0.000000 to 1.000000

Unfortunately both datasets didn't have a common OBJECTID at any point.

I have joined the speed patterns to the streets using the linked field and then got the objectid of the street and pushed that into the EdgeID of the Patterns 3 but this doesn't build the traffic into network.

I have also tried going back to using the tables from the ESRI San Diego Exercise 10 and replacing some sample values but still have the same issues with the error SourceName: Streets_Patterns, ObjectID: 2, Network edge not found with positions 0.000000 to 1.000000

I have also tried naming every table to be the same as the exercise.

References

https://community.esri.com/thread/106547
https://desktop.arcgis.com/en/arcmap/10.3/guide-books/extensions/network-analyst/common-build-errors.htm
https://desktop.arcgis.com/en/arcmap/latest/extensions/network-analyst/traffic-historical-10-1-and-later.htm
Error when building network dataset in ArcGIS Pro?


texts: How to import raster from PostGIS into Python as dataframe

texts:
I am trying to import the raster from PostGIS into a raster for further analysis. I know that I should use psycopg2 to connect to DB, but I do not know how to query the raster from DB, and then convert it into dataframe. Each column should show only the pixel value of the one raster

The aim to do classification, and I need to import single and multiband rasters.


texts: How to connect stops to road for network dataset?

texts:
I have three files that cover a bus route network, road network and bus stops.

I want the bus stops to be the connections between the road network and the bus route.

When I try and build the network dataset I am getting a standalone user-defined junction error. I tried snapping the stops to the bus route and building a network dataset with just these 2 files but still got the same error.

I followed the instructions from here Snapping points to lines in ArcGIS Desktop and automate using VBA? to do that.

Does anyone have any suggestions on what I can do to fix this?


texts: What parameter's data type is Short integer in arcgis 10.3 model builder?

texts:
I am defining a parameter for a python script in model builder. The script calculates a single value from a featureclass using a formular referencing values of a long integer field in the same table. This value will be used as a variable for another tool in the model which is defined as a short integer. I therefore need to set the parameter in the script setting to "output" so that I can use the variable elsewhere and set the data type.
Which data type should I use for small integers?


texts: Only Merge Line segments that don't overlap

texts:
I have 2 road/ line datasets. A & B. A is less detailed than B. I would like any line segment in B that is not in A, to be added to A. But only the lines not in A I do not want to add in overlapping lines.

I'm using QGIS


texts: I have the problem where my two map layers overlap each other, meaning that part of the map is hidden by the white area of the other layer. Thus I need to cut the excess white area to let the bottom layer show.


texts: Python: Cut out the square from the raster. The square has a 30 kilometres (18 miles) side length and has the center in the given coordinates

texts:
The following code displays the coordinates for a city in Poland given by the user, for example Kock.
import subprocess
import sys
from pyproj import Proj, transform
import warnings

warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=DeprecationWarning)

name = input("Enter the name of the object (the town, e.g. Lublin or Kock): ") #Lublin, Kock

query = "SELECT ST_Centroid(ogr_geocode('{}, województwo lubelskie'))".format(name)
command = ["ogrinfo", ":memory:", "-q", "-sql", query]

proc = subprocess.run(command,
                      stdout=subprocess.PIPE,
                      encoding="utf-8")

if proc.returncode:
    sys.exit("Failed to execute ogrinfo command.")

lines = proc.stdout.split('\n')

found = False
for line in lines:
    line = line.strip()
    if line.startswith("POINT"):
        found = True

        _, x, y = line.split()
        x = x.lstrip("(")
        y = y.rstrip(")")
        print("")
        print("{}:\n".format(name))
        print("WGS84 (EPSG: 4326) coordinates")
        print("longitude:  ", x)
        print("latitude:  ", y)
        print("")
        break
if not found:
    sys.exit("City not found.")

inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:2180')
x1,y1 = x,y
x2,y2 = transform(inProj,outProj,x1,y1)
print("PUWG 1992 (EPSG: 2180) coordinates")
print ("X: ",x2)
print ("Y: ",y2)

Output
Enter the name of the object (the town, e.g. Lublin or Kock): Kock

Kock:

WGS84 (EPSG: 4326) coordinates
longitude:   22.4458062648351
latitude:   51.654282346758

PUWG 1992 (EPSG: 2180) coordinates
X:  738269.852848907
Y:  426492.45745886303

In Output, the script displays the coordinates of the city in EPSG:4326, and then in EPSG: 2180.
I have a DEM for all of Poland in the Poland.tif file, in the coordinate system EPSG: 2180.
How to extend the existing script so that from the Poland.tif raster cut a square fragment (with side length of 30 kilometers [18 miles]) with the center where the EPSG:2180 coordinates?


texts: How to turn sqlite lat/lon columns to spatialite geometry?

texts: I have a sqlite database with a table that has a column for longitude and one for latitude. Is there a way to make this data base 'spatial'?


texts: Summarizing Attribute Table Based on Condition

texts:
I have a one-to-many joined table (i.e. duplicated AssetIDs with other different attributes), and I would like to consolidate the table (20 records) into one single record.

However, when using Summarize it takes the max compare score (3), but it always takes the lowest near_dist (0.015821 instead of 0.215665)

I have tried the summary statistics tool and dissolve, but cannot get the 2 conditions correct.


texts: Why don't GRASS tools work in SEXTANTE? (QGIS, MacOS)

texts:
I've got both SEXTANTE and GRASS installed on my platform. GRASS-SEXTANTE report an error when I'm trying to create the Modeller, however, other tools (R) work fine in SEXTANTE. Do these tools have separate configurations?
Please let me know about your experiences.
Best, N


texts: How to fix distortions in rendered OpenStreetMap maps?

texts:
I use OpenStreetMap in order to generate maps and add some cool stuff in it. Actually I have good results but I seem to have a problem with reporting a good image ratio when I create my generated file.

I put OpenStreetMap data in a PostgreSQL and I use PostGis to work with it in the EPSG:4326 projection. Everything is fine except this ratio problem.

Actually I've tried almost every things (except the good one indeed)... that's why I'm looking for your help. I think I forgot something.

In my program I use the min lon & min lat (as origin of the selection) ans max lon & max lat (as selection max coordinates)... basically I though a correct ratio would have been given by something like that (using N as size factor):

N = 10000
x_size = maxlon - minlon
y_size = maxlat - minlat
file_width = x_size * N
file_height = y_size * N

The most intriguing part is that this computation works for some parts of the world. In my examples below, the three cities use this ratio computation, but even if Montreal and Lyon don't work, Stockholm seems to work pretty well... The two others seem stretched.

Lyon ( http://www.openstreetmap.org/?lat=45.7731&lon=4.8544&zoom=14&layers=M ):

Stockholm ( http://www.openstreetmap.org/?lat=59.3278&lon=18.0616&zoom=14&layers=M ):

Montreal ( http://www.openstreetmap.org/?lat=45.5088&lon=-73.5878&zoom=13&layers=M ):

I think it's a noob problem... something silly enough to be barly invisible after hours of work on it. So I think you will find my problem easily.

Thank you.


texts:
ro<-readShapeSpatial("roads")
ra<-readShapeSpatial("railways")
co<-readShapeSpatial("LIE_adm0")
u<-union(union(extent(ro),extent(ra)),extent(co))
rast<-raster(u,res=1)
r1<-rasterize(ro,rast)
r2<-rasterize(ra,rast)
r3<-rasterize(co,rast)
rMerge <- merge(r3, r2, r1, fun=sum, na.rm=T)
plot(rMerge)

question 1

I've heard readOGR is better than readShapeSpatial, can u tell me when and why is better?

question 2

Why the script doesn't return the map but it returns this thing below? I have already change the resolution from 1 to 10 and 100 and still the same.
.


texts:
I am trying to add a NetCDF file of sea ice age into QGIS (have tried by raster, drop and drop, mesh), the .nc file contains the sea ice age, latitude and longitude. When I by adding by raster it asks which layer I want to add out the the three so I select sea ice age however it does not project it to my current map scale - it only appears very very small and therefore unusable. Why is this happening?
This is my .nc data source: https://daacdata.apps.nsidc.org/pub/DATASETS/nsidc0611_seaice_age_v4/data/


texts: Zonal Statistics doesn't work after updating processing (2.10.2)

texts:
After updating the processing tool zonal statistics doesn't work.

  in method 'SpatialReference_ImportFromProj4', argument 2 of type 'char *

Could someone else try if it's a bug or just an issue on my PC.


texts: Get max value from multiple points surrounding buildings

texts: I have a series of building represented as polygons. Surrounding these buildings I have multiple points with heights associated with them. What I am trying to do is identify all the points that are within 2 meters of each building and assign the maximum value from those points to the building polygon. The image below show what I am trying to achieve. I think I should be using spatial join, but not sure how to get the maximum values from a group of points.


texts:
I'm using google earth engine to do some simple habitat classification and now want to export the result as a KML or KMZ. When I run this code it is layer 2 that I want to export - the habitat classification map.

  //sat image
    var S1 = ee.ImageCollection("COPERNICUS/S1_GRD"),
    S2 = ee.ImageCollection("COPERNICUS/S2"),
    fc = ee.FeatureCollection("ft:1jvw7bQoChrt-3dk4bQC96TMjZnnyrw6DvaLpfEdk"),//polygon
    fc2 = ee.FeatureCollection("ft:1_uzNGOXOEdfnRmkIYIkVCiMHGx8gCr-1ywGyhjLF"),//points
    fc3 = ee.FeatureCollection("ft:1WXJEHBhLLe1ztmsjzt7WyMe1Go04Uuj96XEQ-C0-");//study area
    //Set a point of interest
    var aoi = fc;
    //Set bands of interest
    var bands = ['B2', 'B3', 'B4','B8', 'B11', 'B12'];

//Import and filter Sentinel 2 collection
var S2 = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2017-06-20', '2017-07-20')
.filterBounds(aoi)
 print(S2);
var image = S2.mosaic();
// Clip the image to the polygon geometry
var landcover_roi = image.clip(fc3);

Map.addLayer(landcover_roi, {'bands': 'B12,B8,B4','min': 0,'max': 4000},        'Sentinel 2A SWIR,NIT,R', 0);
var training = landcover_roi.sampleRegions(fc2, ["name"], 59);
var classifier = ee.Classifier.cart()
.train(training, "name", landcover_roi.bandNames());
var out = landcover_roi.classify(classifier);
Map.addLayer(out, {palette: [
'ff0000',// not forest
'387242',// forest
'e0e0d1',// volcano
], min:1, max:3});

// Create a geometry representing an export region.
var geometry = geometry;

Map.addLayer(ee.Image().paint(fc, 1, 3));
var fc = fc2;
Map.addLayer(fc2);

 //Create a band for each veg type
 var class1 = out.eq(1).select([0],['not forest']);
 var class2 = out.eq(2).select([0],['forest']);
 var class3 = out.eq(3).select([0],['volcano']);
 // Concatanate the bands into a single image
 var total = ee.Image.cat([class1,class2,class3]);

 // Export the FeatureCollection to a KML file.
Export.image.toDrive({
  image: total,
  description: 'imageToDriveExample',
  scale: 30,
  region: geometry
});


texts: Convert Google Street View coverage layer in QGIS to a Maphub.net friendly format

texts:
I have added an XYZ tile layer to QGIS with this address: https://mts1.googleapis.com/vt?hl=en-US&lyrs=svv|cb_client:apiv3&style=40,18&x={x}&y={y}&z={z}
I want to export the coverage layer in order to use it on a map that I have in maphub.net, this website only accept importing data from GeoJSON, KML/KMZ, GPX, Shapefile ZIP and IGC files or JPG and PNG images.
What format I need to use? and how can I do this?


texts: Select and count street pattern features (morphology) using ArcGIS for Desktop?

texts: I'd like to be able to select and count the number of cul-de-sacs, lollipops, dead ends, or other types of street structure in a hillside community.


texts: Looking for the simplest way to get elevation for a give wgs84 coordinate

texts:
I started a project and write an application in java using geotools. I have an simple geoobject with the attributes latitude, longitude and altitude. Now i have to determine, wheather the geoobject is an aircraft or a car on a road.

I modeled my roads as LineStrings an my geoobjects a circle polygones. So i run a simple test, if the geoobject touchs or intersects the linestring, at the given coordinate. So far it works.

But some roads have an elevation. I found sample digital elevation model files, but don't know to handle them. Is there an easy way to get elevation data respectively is there another simple way? Need i software like postgis or somethin else?


texts: Why does AddField give ERROR: 000732 Input Table: does not exist or is not supported?

texts:
I was handed off a a bunch of state shapefiles and wrote a script to add a few new fields and calculate some attributes. The first line in my script is -

arcpy.AddField_management('C:\WB_prj\cityPts\ghosts.gdb\vietnam_adm3_citypt_albers_ghost', 'City_cnt', 'SHORT', '#', '#', '#', '#', 'NULLABLE', 'NON_REQUIRED', '#')

I have all of the state shapefiles in a geodatabase. The script runs perfect for some of the files but for others I get an error that is shown below. My pathways are correct and I can manually add the fields through ArcMap but when I try and run the script for some of them I get the error. Any ideas on how to solve this? I'm running ArcGIS 10.0

  File "C:\WB_prj\scripts\ghostRegions_final.py", line 79, in
      add_fields()   File "C:\WB_prj\scripts\ghostRegions_final.py", line 14, in add_fields
      arcpy.AddField_management('C:\WB_prj\cityPts\ghost.gdb\vietnam_adm3_citypt_albers_ghost',
  'City_cnt', 'SHORT', '#', '#', '#', '#', 'NULLABLE', 'NON_REQUIRED',
  '#')   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\management.py", line 2625, in
  AddField
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Table: Dataset C:\WB_prj\cityPts\ghost.gdb
  ietnam_adm3_citypt_albers_ghost does not exist or is not supported
  Failed to execute (AddField).


texts: Two inner polygons touching and covered by two different outer polygons

texts: I need to create a topology where 1 and 2 touch but don't overlap, a and b touch but don't overlap, and 1 covers a and 2 covers b.  Is it possible to achieve this?  I've tried various snapping options, but they all seem to result in one of the inner polygons overlapping with a wrong outer polygon.  I'm using QGIS.


texts: How to find a datasets spatial reference and apply it to another dataset?

texts:
I'm struggling to use the describe function to find the spatial reference of a given shapefile and then apply it to another shapefile which has been copied to a projected folder.

I'm relatively new to using python so I apologize if my code attempt is dismal.

This is where I'm at:

import arcpy

data = arcpy.GetParameterAsText(0)
targetProjection = arcpy.GetParameterAsText(1)
projected = arcpy.GetParameterAsText(2)

# Create an error class
class NoExist(Exception):
    pass

try:

 # copy dataset to projected dataset folder
    arcpy.copyFeatures_management(data, projected)

 # describe the spatial reference of the dataset in the target projection folder
    desc = arcpy.Describe(targetProjection)

    # obtain the spatial reference
    spatialRef = desc.spatialReference

I think I've worked out how to find the target projection but am unsure of how I apply that to my other dataset.


texts: Filter for first geometry that is crossed by a line

texts:
I have a query like this

SELECT p.id
FROM polygons p
INNER JOIN lines l
ON ST_Crosses(p.the_geom), l.the_geom)
WHERE p.id = l.id;

lines can cross more than one polygon. However, I want the inner join only to be performed on the very first polygon that is crossed. How would I do that?


texts: How to add Geo PDF header to a PDF rendered with mapnik?

texts:
I rendered a PDF with mapnik and I am now trying to add geospatial header information to the PDF.

I tried the method add_geospatial_pdf_header() from mapnik.printing.PDFPrinter. I do get some Geo PDF information in my PDF but when I check it with gdalinfo or with qgis, I get some invalid angles in my corner coordinates (and the ones that are not marked as invalid don't look correct either).

What I do in my code is something like that (where my pdf extent is calculated based on my layer envelope):

m = mapnik.Map(100, 100)
m.zoom_to_box(mapnik.Box2d(my_mapnik_pdf_extent))
page = mapnik.printing.PDFPrinter()
page.render_map(m,"my_output_file.pdf")
page.finish()
page.add_geospatial_pdf_header(m, "my_output_file.pdf", epsg=2193)

Any idea on what I might be doing wrong or any experience with add_geospatial_pdf_header to share? (I checked my projection multiple times and that side seems ok)

Source code from mapnik: https://github.com/mapnik/mapnik/blob/2.3.x/bindings/python/mapnik/printing.py


texts: Optimizing PostGIS k-nearest neighbour function for large amount of polygons?

texts:
I created a function for identifying the K nearest neighbours for every municipality, which works fine, but when I use this function for a large amount of municipalities (over 2900) it takes over 2 minutes. The function is defined as:

CREATE FUNCTION KNN(k int,m_id int) RETURNS SETOF integer AS $$
SELECT  mun2.id AS m2_id
  FROM municipality AS mun1, municipality AS mun2
  WHERE mun1.id = m_id AND mun1.id != mun2.id
  GROUP BY m2_id, mun1.geom
  ORDER BY mun1.geom <#> mun2.geom
  LIMIT k
$$ LANGUAGE sql

Then I run query

SELECT m.id, KNN(30,m.id)
FROM municipality AS m
ORDER BY id

Is there any way to optimize my function?


texts: Finding distance between neighboring points in ArcGIS Desktop?

texts:
I have two shape files: one has point features indicating urban centers, the other is a polygon features file that has counties. Each county has one urban center. There are around 2800 counties. I want to find the distance between one urban center in a county to all of its neighboring urban centers (by this I mean centers located in counties that border the selected county)

For example,
if county 1 borders county 3 and county 2, I would like to find the distance between urban center 1 (located in center 1) and center 3, and between center 1 and center 2. If county 7 borders county 8, 9 and 10, I want to find the distance between urban center 7 and centers 8,9 and 10.

So far I have tried the following:

I have tried using the Near Generate Table and the Near tool. The problem with this is that I do not only want to find the distance between points that are closest to one another, but rather find the distance between a point and all of its neighboring points.
I used the point distance tool, which gives me the point distance between every single point for all 2800 points, and this is very redundant and inefficient, since I only want to find the distance between neighboring point features.

Is there any way to do this in ArcGIS Desktop 10.3?


texts: Creating Layer for Selected Features, then Moving it using ArcGIS for Desktop?

texts:
I have a problem with ArcGis 10.1 which I rarely use.

In an edit session I select a polygon. I would like to 'Create Layer for Selected Features' and can do this; the new layer appears in the table of contents. When I try and move the new layer, the old layer moves with it.

When I click on the polygon the option of selecting the two coincident polygons appears but I can't get these to separate. It seems to be impossible to do anything to one feature and not the other.

The strange thing is that yesterday the program wasn't behaving like this.

Can anyone can explain what is going on?


texts: Sampling stratified random points using R for accuracy assessment

texts:
I want to sample stratified random points using R. I found the sampleStratified function from the R raster package.

This functions samples the same amount of points for each class of my classification, but I want to have points that are randomly distributed within each class, where each class has a number of points proportional to its relative area.

Can anybody suggest another way of doing that?


texts: QGIS: Load GPX polygon (as polyline) creates additional lines (GPS Tools)

texts:
this is slightly confusing: I want to upload a polygon to my Garmin Etrex Vista HCx in QGIS 3.10.3 using GPS Tools.

I have now managed
- to convert my polygon to GPX (convert to polyline, save as GPX, specs: CRS = EPGS:4326 – WGS84, FORCE GPX TRACK = YES, unselect all attributes)
- to upload the GPX track to my Garmin device using GPS tools.

The GPX track that is exported from my polyline looks just fine and matches my polygon boundaries. However, when I load the GPS using GPS tools so that I can upload it to my Garmin (without that step it doesn´t show up in the list for upload), I get a CONTINUOUS polyline, all points are somehow connected to each other, resulting in this picture (yellow being my polygon and red the resulting track):

How can I avoid this behaviour?


texts:
Some months ago, I had an old version of QGIS (I don't remember exactly the version, probably < 3), and when I opened an mbtiles file, it was easy to find all the layers, including transportation, place, water, etc. I was also able to export a specific layer to csv. (See my previous post)
Now with QGIS 3.14, it is no longer able to see the layers, and of course I'm not able to export data to csv files (see figure below).
How can I get back the layers as before?


texts: How to encode shapefiles from LATIN1 to UTF-8?

texts: Is there a way to encode shapefiles from LATIN1 to UTF-8?


texts:
I'm building an application that can accept a shapefile/mapinfo/feature class/spatialite table/etc.  So - various formats from various vendors/organizations.

So let's say the user is going through a wizard to import data (or visualizing that data in plan view)...is there generic term for spatial data?  Something that a typical user would understand that would differentiate it from a text file or standard table?  This isn't a GIS product - it just makes use of GIS data.  It is not limited to use in one area of the world.

Anyone have any experience with lowest common denominator terminology for beginner/intermediate users (advanced users would just figure it out I suppose)?  What is most obvious?

Spatial data?

Spatial file?

Geospatial data/file?

Geographic data/file/table?

Data with one weird column/file that I can't understand?


texts: Why are labels visible through a basemap?

texts: I'm working with a custom basemap I created as a vector-tile and published it to AGOL, but I noticed labels are visible through it in both AGOL and ArcGIS Pro. Can anyone explain why this behavior possible? In the image, you can see that the basemap is covering my other layers, yet I can still see the labels for the building and safe zone layer.


texts:
I'm new on QGIS, I have shapefile for the whole country and also have for Province down to City/Municipalities, is there any way on attribute table where I can filter using an expression for example:

"PROVNAME" = 'CEBU', but this only for one area.

I tried adding multiple areas, like  "PROVNAME"  =  'CEBU' 'others' 'others'.

So basically tried that multiples areas, but getting an error:

Parser Error:
syntax error, unexpected STRING, expecting $end

Eval Error:
No root node! Parsing failed?

Is there anyway I can select multiple areas by expression, so I can extract from a large area shapefile, I'm trying to filter it from country down to province and cities, so I can extract it easily.


texts: ArcGIS ArcPy: Iterate over shapefiles and report field "Desc" Length

texts:
I've just started learning to script using python, i.e. no programming experience. I need to iterate over a folder that contains multiple shapefiles that all start with nlc... I need to iterate over each shapefile and then report the length of the field "Desc". The reason that I need to do this is that the field lengths are not the same and I've generated a model in modelbuilder that iterates over the folder, generates a feature layer (in_memory) filtering the output layer by using a sql expression to return a specified set of records and then append the output records into a single feature class. I only realised half way through running my model that the values in certain shapefiles are larger than the initial field length that I set for the new feature class.

An example to get me started, would really be appreciated, thanks in advance.

Peter

The following is the initial python script that I wrote, not realising that I didn't need to loop through the fields list:

    import arcpy
... arcpy.env.workspace = r'E:\Projects\H109009\data_received\DOA\NLC2000\tiles'
... fcs = arcpy.ListFeatureClasses("nlc*","polygon")
... for fc in fcs:
...     fields = arcpy.ListFields(fc,"Desc*")
...     for field in fields:
...         print("{0} is a type of {1} with a length of {2}"
...           .format(field.name, field.type, field.length))


texts: Defining and projecting batch of raster images using ArcPy?

texts:
Edit: the error was caused by spelling .tif as .tiff`.

The script ran fine before I added InSR, the input coordinate system.

I'm getting the following error:

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000445: Extension is invalid for the output raster format.
Failed to execute (ProjectRaster).

Here is my code:

import arcpy, sys, os

InFolder = r"E:\ReprojectionScript"
OutFolder = r"E:\ReprojectionScript\test"
OutSR = arcpy.SpatialReference(26918) # NAD83 / UTM zone 18N
InSR = arcpy.SpatialReference(26996) # Missouri East
arcpy.env.workspace = InFolder

rasterList = arcpy.ListRasters()
for Ras in rasterList :
     inRaster = os.path.join(InFolder, Ras)
     outFileName = "{}.{}".format(os.path.splitext(Ras)[0], "tif") # Remove .sid and add .tif to Ras filename
     outRaster = os.path.join(OutFolder, outFileName) # Path and Filename of output raster (.tif)
     arcpy.ProjectRaster_management (inRaster, outRaster, OutSR, "", "", "", "", InSR)


texts: Save a selection in QGIS and keep the alias of field names?

texts:
I am trying to save selected features into a new shapefile using QGIS 10.2.  The new file does not contain the aliases for the field names.

Is there a way to maintain the aliases in new files?


texts: Assign points with varying radii to central meeting locations

texts:
I have a set of points in a city, each with a maximum distance radius of 0.5 - 1.5 miles.
From this set of points, I'd like to create a small number of 'meeting points', i.e. locations that are central to a multiple points while not violating the constraint of <= maximum distance raidus). Those points far away from the others will have a meeting point of 1.

My first attempt was to buffer each point and find the intersection with the highest number of overlaps, but it proved very computationally expensive because I compared each new intersection with all other nearby intersections to determine the points reachable within the area.

My second attempt was to fill each buffer with possible meeting points every 100 meters & run a distance calculation on each one to determine which points they could reach, but again, computationally expensive and prone to miss edge cases.

I'm now thinking of creating a centroid of all points where at least 1 point is within radius(p1)+radius(p2) of all the other points, which would work as is, if they all had the same radius. Since the radii vary, I thought I could adjust the centroid by moving it closer to the mean of the points it doesn't reach until it either fits or not.

Am I overcomplicating this? It seems like there should be a more elegant solution, but I can't think of it. Any pointers in the right direction would be very welcomed.


texts: Copying/replacing first letter in column of attribute table in ArcMap with field calculator and Python Parser?

texts: How to copy / replace the first letter in a column if that letter is "N" and replace with "C" in the field name "DISTRICT" inside the attribute table in ArcMap using Python inside the field calculator?


texts: Checking two different releases of POI Features for not only spatial differences using ArcPy?

texts:
I have Point Of Interest (POI) data for two different years. I wrote code to get two feature classes where each contains the POIs that differ (spatially) from the other release. I make mostly use of the arcpy.SelectLayerByLocation_management.

# Get intersection from the "new side"
arcpy.SelectLayerByLocation_management(new_layer, "ARE_IDENTICAL_TO", old_layer, "", "NEW_SELECTION")
# This selection contains NEW WITHOUT OLD
arcpy.SelectLayerByLocation_management(new_layer, "", "", "", "SWITCH_SELECTION")
# NEW WITHOUT OLD feature layer output
arcpy.CopyFeatures_management(new_layer, new_diffs)
# Get intersection
# First switch the selection from the previous layer again
arcpy.SelectLayerByLocation_management(new_layer, "", "", "", "SWITCH_SELECTION")
# Then proceed as above
# Get intersection from the "old side"
arcpy.SelectLayerByLocation_management(old_layer, "ARE_IDENTICAL_TO", new_layer, "", "NEW_SELECTION")
# This selection contains OLD WITHOUT NEW
arcpy.SelectLayerByLocation_management(old_layer, "", "", "", "SWITCH_SELECTION")
# OLD WITHOUT NEW feature layer output
arcpy.CopyFeatures_management(old_layer, old_diffs)

It works ok, but I don't got the cases where the POI location may have changed a tiny bit (maybe to a more accurate position) but the POI itself is the same. I would NOT want to have these included in the output, or in other words, I would not classify them as a real difference. So I additionally might want to compare with something like a arcpy.Buffer_analysis and check for additional POI's names if they are equal to the one which differs spatially.

Is there an easy way to integrate this functionality into my code?

The unique IDs for the POIs change from release to release, so they are not 100% identifiable throughout different releases (yes, really, I also couldn't believe it). But it should be enough for me to just check for equal names.

EDIT:
I drew a picture to make it a Little bit more clear what I want to get.


texts: Creat point count trigger in PostgreSQL?

texts:
I wish to create a trigger in PostgreSQL for counting points in a polygon.

The tweak is that the count should automatically get updated whenever the polygon is edited.

Can anyone please help me on the same?


texts: GeoJSON shape distorted on Leaflet.js spherical mercator WMS map

texts:
I have a map with a simple WMS layer in spherical mercator projection (rendered by a OSM-Mapnik-TileCache server stack, everything is in EPSG:900913 / EPSG:3857).

First: I can't understand why, if Leaflet uses CRS EPSG:3857, whose unit is meters, I still have to provide coordinates, like maxBounds and center, in WGS84 format, which is degrees. Is it an API inconsistency?

Second: I need to render a GeoJSON feature expressed in WGS84 coordinates (EPSG:4326), but the shape is visibly distorted (see pics below). Should I convert the shape from EPSG:4326 to EPSG:3857? How? If I did so, I would have coordinates expressed in meters, while Leaflet still expects degrees!

I made a double check with a simple wms layer in OpenLayers and a GeoJSON overlay and it works fine, so I proved that the distortion is not in the data.

Any clue?

Here is the test GeoJSON, hope someone can understand the reason of the issue:

{
    "type": "Polygon",
    "coordinates": [
        [
            [14.740017498458682, 40.673078870109705],
            [14.740132563378529, 40.673283531348574],
            [14.741625561383819, 40.67246759585111],
            [14.741671820229074, 40.672308542443076],
            [14.74151049646397, 40.672262934612235],
            [14.740017498458682, 40.673078870109705]
        ]
    ]
}

GeoJSON feature on Leaflet map:

GeoJSON feature on OpenLayers map:


texts:
I have written a very simple C++/Qt application using QGIS comprised of a single "main.cpp" file as below:

#include <qgsapplication.h>

int main(int argc, char ** argv)
{
    QgsApplication app(argc, argv, true);
    return 0;
}

and my Qt Creator project file is defined as:

QT       += core gui

greaterThan(QT_MAJOR_VERSION, 4): QT += widgets

TARGET = qgis_test
TEMPLATE = app

CONFIG += c++11

SOURCES += main.cpp

HEADERS  +=

LIBS += -L/usr/local/lib -lqgis_core

QGIS_DIR = /home/rouhollah/Programs/qgis-2.8.1

INCLUDEPATH += $$QGIS_DIR
INCLUDEPATH += $$QGIS_DIR/src/core

DEFINES += CORE_EXPORT= GUI_EXPORT=

The program is compiled and linked, but exits unexpectedly when run with the following error message:

realloc(): invalid pointer: 0x00007fc3d802e0e0

Even running the program in the Debug mode, the debugger stops before hitting the breakpoint at the first line of the main function.

Note that I have installed the QGIS libraries by compiling the source codes obtained from GitHub.


texts: Looking for interpolation tool to produce river bathymetry (topography) from sections

texts:
As per the title, I am looking for a tool that can produce bathymetry of a river from cross sections. Although people who may know such tools probably will not need any explanation why this kind of interpolation is a bit complicated, here are a few words of explanation for those which may know and use such things for other purposes.

In order to receive a DEM of river valley with river bed bathymetry, I need to fence/mark with barriers the area on my original DEM where bathymetry is to be created.  Basically it is the whole area of the wet river bed where TLS did not have a chance to measure the bottom surface (bathymetry). Rivers are in various shapes, sometimes more meandering, sometimes more straight, and also the depth of the bottom and its alignment are varying a lot, which makes interpolation more complex as it has to follow the centre line of the river and also conduct the interpolation vertically and horizontally.

Dr Venkatesh Marwade has prepared such extension to ArcGIS 9.2, though I did experience some problems and was just wondering if anyone would know a piece of software or code that could generate something alike.

(Edit: see graphics in this tutorial.)


texts: Error in the positioning of a point with L.GeoJson

texts:
I am having a problem and I cannot understand why.

This is the code :

const map = L.map('map');

const tileURL = 'http://{s}.tile.osm.org/{z}/{x}/{y}.png';
L.tileLayer(tileURL).addTo(map);
map.setView([-34.603722, -58.381592], 8);

var points = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "properties": {
                "idSitio": "46",
                "emplazamiento": "ARCF0628",
                "codigo": "CAGVR     ",
                "enplanmtto": "No",
                "nombresitio": "AGUAS ARG. - VARELA 1835",
                "categoria": "C",
                "area": "AMBA",
                "propiedad": "TERCEROS",
                "direccion": "Varela 1835 Capital Federal",
                "contrato": "TELEFONICA"
            },
            "geometry": {
                "type": "Point",
                "coordinates": [
                    -34.6475166,
                    -58.4511706
                ]
            }
        }
    ]
};

var layer = L.geoJson(points).addTo(map);

var marker = L.marker([-34.6475166,-58.4511706]).addTo(map);

Points are at the same Latitude and Longitude but on the map they are positioned in different places

Any suggestions?


texts: Which .las classification codes and returns to use when creating DTM and DSM?

texts:
I am trying to figure out what classification codes and returns that I should use to create DTM and DSM?

In my layer --> "BrooksCamp.lasd" I have these classification codes and returns:

(and then 26 to 31 reserved)

I am using the LAS Dataset to Raster tool.


texts: How to access clicked coordinates in a Custom Python Expression Function?

texts:
My final goal is to be able to digitize a point layer (create new point features) and to have as a default for one of the attributes an attribute from another layer that is at the same location.

Up until now I figured how to access the layer/features/geometry/coordinates/attributes from the source layer. Then I thought that I can just use the feature from the function definition to access its coordinates to compare them with the coordinates in my source layer (they should be the same, as I'm clicking on the same points).

@qgsfunction(args=0, group='Custom')
def Extract_Attribute(value1, feature, parent):

But, this feature actually refers to any commited feature, at the moment of the creation of the feature, the feature only exist somewhere in the memory, but even in memory I think that the feature is created only after I fill in some attributes and click ok. I think that what I need to access are the coordinates that the mouse clicked (or snapped on) last. From the documentation it is clear that this is what I have to use QgsPointLocator::Match QgsMapMouseEvent::mapPointMatch()

The problem is that I have only seen this function being used when the event is already passed to the function, which is not the case for a custom express function. How can I do it anyway?

This is what I have until now:

"""
Define new functions using @qgsfunction. feature and parent must always be the
last args. Use args=-1 to pass a list of values as arguments
"""

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

@qgsfunction(args='auto', group='Custom')
def extract_elevation(source_layer_name, feature, parent):
    source_layer=None
    for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
        if lyr.name() == source_layer_name:
            source_layer = lyr
            break
    source_features=source_layer.getFeatures()
    extracted_attribute = None
    clicked_coordinates = 'something like QgsMapMouseEvent().mapPointMatch()'
    for source_feature in source_features:
        source_geometry=source_feature.geometry()
        #clicked_coordinates = source_geometry.asPoint()#this line simulates a click on the first feature of the source_layer
        if clicked_coordinates == source_geometry.asPoint():
            extracted_attribute = source_feature['some_attribute']
            break
    #QgsMessageLog.logMessage(str(clicked_coordinates))
    return extracted_attribute


texts: Making Parameter List Update/Refresh in Python Toolbox?

texts:
I am using a python toolbox in ArcGIS 10.2.2 and want a drop down list of all the versions for the person reconciling to choose from. (And I do not want to include dbo.DEFAULT and DS.Draft). But this list is populated right after a script is run to create the editor's version and if the Python Toolbox isn't refreshed manually in ArcMap, then this newest version doesn't show up in the list.

I looked at updateParameters, but either I was writing it incorrectly or it doesn't work because the user isn't actually changing a parameter that triggers it.

Regardless, I want to refresh the list at the time the script is run.

Can I use updateParameters or do a reload or something else snazzy that will make the list repopulate?

The other thing is I could switch the order of the parameters and having the user type their initials first which would then be a user action that could trigger the updateParmaters, but not sure this is the issue and seems convoluted.

import sys, arcpy, os

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "EGDB Data Steward Version"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Prep]

class Prep(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Prepare for Reconcile and Post"
        self.description = "Create grandfather version of editor for DS"
        self.canRunInBackground = False

    def getParameterInfo(self):
    #Define parameter definitions
        params = []

    # First parameter
        param0 = arcpy.Parameter(
            displayName="Editor's Version Name",
            name="editor_full_versionname",
            datatype="GPString",
            parameterType="Required",
            direction="Input")

        param0.filter.type = "ValueList"
        param0.filter.list = [v.name for v in arcpy.da.ListVersions(r"\\inpyosegis05\DatabaseConnectionFiles\DataStewardsOnly\YOSEGIS_VectorYOSE_DS.sde")
                              if v.name not in ["dbo.DEFAULT", "DS.Draft"]]  #Take DEFAULT and Draft off the drop-down list as choices.

    # Second parameter
        param1 = arcpy.Parameter(
            displayName="Data Steward's Initials",
            name="DS_initials",
            datatype="GPString",
            parameterType="Required",
            direction="Input")

        params = [param0, param1]
        return params

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    param0 = arcpy.Parameter(
        displayName="Editor's Version Name",
        name="editor_full_versionname",
        datatype="GPString",
        parameterType="Required",
        direction="Input")

    param0.filter.type = "ValueList"
    param0.filter.list = [v.name for v in arcpy.da.ListVersions(r"\\inpyosegis05\DatabaseConnectionFiles\DataStewardsOnly\YOSEGIS_VectorYOSE_DS.sde")
                          if v.name not in ["dbo.DEFAULT", "DS.Draft"]]  #Take DEFAULT and Draft off the drop-down list as choices.
    return


texts: QGIS2Leaf Plugin installation error "cannot import name AlgorithmDecorator"

texts:
I have been using the qgis2leaf plugin and succeeded in creating interactive maps, adding legend, photos, etc. But suddenly I noticed that the OSM base layer has disappeared and the polygons are not getting published the way they are styled, so i tried reinstalling the plugin. And this is the message I am getting :-

I am unable to install the plugin. I am using QGIS 2.6.1 Brighton.


texts: ArcGis online app - editor stopped working

texts:
I am trying to use ArcGIS Online.
I have a web map application but all of a sudden I am unable to create points when editing.

I've tested the map service and editing via arcmap without any issues.

When I do try and create a point I get the following error in the google chrome console:

Error: Unable to complete operation.

at Error (native)

at Object.g.load (http://js.arcgis.com/3.16/init.js:1042:343)

at http://js.arcgis.com/3.16/init.js:181:69

at c (http://js.arcgis.com/3.16/init.js:66:221)

at d (http://js.arcgis.com/3.16/init.js:66:10)

at resolve.callback (http://js.arcgis.com/3.16/init.js:67:350)

at c (http://js.arcgis.com/3.16/init.js:66:436)

at d (http://js.arcgis.com/3.16/init.js:66:10)

at resolve.callback (http://js.arcgis.com/3.16/init.js:67:350)

at http://js.arcgis.com/3.16/init.js:185:300

the first time I examined the console error log I noticed dojo errors (referenced back
to the arcgis.com site)

I've eliminated proxy and firewall issues and I have tried to edit now on three different
browsers.

Any ideas on how to resolve this issue?


texts:
This opens just fine in the browser and the server seems to respond positively to the WFS request.. but nothing loads. Can someone help me out here?

<title>LeafletTest</title>
<link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-
0.7.3/leaflet.css"/>

<meta http-equiv="Access-Control-Allow-Origin" content="*">


texts: Attribute data showing (Pop-up) in GeoSever

texts:
How can I show attribute value of shapefile in geodatabase?

<title>OpenStreetMaps and GeoServer</title>
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
 <link rel="stylesheet" href="http://localhost:8080/geoserver/openlayers3/ol.css" type="text/css">

Text


texts: How to create temporary tables in spatialite?

texts:
How do you create a temporary table with geometry column in spatialite?

I have tried:

CREATE TEMP TABLE rivers (code INTEGER, id INTEGER, name TEXT);

Followed by:

SELECT AddGeometryColumn('temp.rivers','geom', 27700, LINESTRING');

The first query succesfully creates the table, but the second query does not create the geometry column. I have tried running the queries in spatialite-gui, and in python with the sqlite3 module, with the same results.


texts: how to save a raster after converted from numpy array?

texts:
I'm trying to save a raster after calculation in numpy array. Here is what I wrote:

rasterSM = arcpy.NumPyArrayToRaster(SM, lowerLeft, cellSize, value_to_nodata=nodata)

rasterSM2 = arcpy.DefineProjection_management(rasterSM, spatial_ref)

rasterSM2.save(SMOutputName)

But the error message says that 'Result' object has no attribute 'save', then out to save a raster from Numpy array?


texts: Opposite of r.carve in QGIS - r.wall or r.dam?

texts:
I wonder if there is a way to perform an operation similar to r.carve, but instead of carving in the stream based on the low elevation on the extremes of a line, it fills a barrier based on the high elevation of the extremes of a line.
The purpose is hydrology, to play with the streams to change the stream network.


texts: Any idea how to export multiple KML files from one SHP layer?

texts: When exporting one KML, it's exceeding the Google Fusion Table's max limit of 100mb per upload. Fusion Tables however lets users merge more than one KML for a merged max size of 250mb. Any idea how to export multiple KMLs from QGIS?


texts: Point Coordinates Conversion with unknown system

texts:
I've been using QGis for few weeks now. Although this site always came very helpful, I can't find solution to my recent problem:

I have few points, and know their coordinates in "WGS-84" and "ISN-93" systems. Those points also have specified coordinates in "Local" system.

The problem is that there are also few points that only have "Local" coordinates, and I would like to find their "ISN-93" coordinates.

How can I do it?

*for example:

point1 ISN93 (685320.253,420614.733) Local (-3537.239,2284.487)

point2 ISN93 (685841.253,420586.536) Local (-3016.528,2248.842)

point3 ISN93 (?,?) Local (-5187.996,2584.508)

I was trying to use simple geometry, but when I measure distance between points, ISN93 system gives slightly different lengths than Local (difference vary between lengths to up to 0,01% of its value), so it's no use.


texts: Extract shapely polygon from matplotlib hexbin polycollection

texts:
I'm trying to extract all the vertices of hexagons that are returned from plt.hexbin. I then want to eventually convert them to shapely polygons, so I can use them with geopandas. At the moment I can only seem to extract one polygon from the PolyCollection returned by plt.hexbin
Does anybody have any advice as to how to return the full grid of hexbins?
This example will reproduce my problem,
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon
from matplotlib.patches import RegularPolygon

geoms = [Point(919.000, 533.500), Point(940.500, 528.500),
 Point(937.000, 555.000), Point(889.000, 526.500), Point(898.000, 520.500)]

l = gpd.GeoDataFrame({'geometry':geoms})

y= np.array([ 533.500,  528.500,  555.000,  526.500,  520.500])
x= np.array([ 919.000,  940.500,  937.000,  889.000,  898.000])

f,ax=plt.subplots(1)
im =ax.hexbin(x,y, gridsize=(3,3), edgecolor='black')
l.plot(ax=ax, color='black', markersize=55)

Which produces this...

But when trying to extract the hexagons it produces a hexagon with the same position..
array_of_hexes=[]
for x,y in im.get_offsets():
    hexes = RegularPolygon((x, y), numVertices=6, radius= 5 )
    verts = hexes.get_path().vertices
    trans = im.get_transform()
    points = trans.transform(verts)
    array_of_hexes.append(Polygon(points))

gpd.GeoDataFrame({'geometry':array_of_hexes}).plot(edgecolor='black')


texts: Having more than one Openlayers layer switcher associated to one map

texts:
I am developing an OpenLayers 3 web application were I would like to incorporate more than one Layer Tree alongside my map.

I have this JSFiddle set up that it has a map with multiple layers and one Layer Tree. I would like to duplicate the layer tree multiple (4) times and have a parameter (as per the one in each layer of the example) that then defines which Layer Tree the layer goes into.

I imagine it involves an extra if statement when building the layer tree at the bottom of the code, similar to the slider option already in there, but I don't seem to be able to make it work.

The solution should also ideally remove the "all layers" group at the top of each Layer Tree.


texts: What are the "must read" academic articles every GIS professional should read?

texts:
There was a related question a while back that dealt with the top GIS blogs, but I'm also very interested in "academic" sources of GIS information. I'm currently updating my repository of information and am also considering some new journal subscriptions.

Every professional field has some big names and big articles. Who are the big names and what articles should everyone be aware of in GIS/cartography/remote sensing?

I'm looking for 3 things in particular:

1) Recommendations on the top sources for GIS articles - important/reliable journals.

2) Major articles from the past - those that paved the way towards the current state of our field.

3) Recent articles showcasing major findings/methodologies.

Bonus: I've recently been searching for a good source of information on cartographic design - any articles related to this would be a great help.


texts: Using Georeference information stored in tfwx (from ArcGIS) in QGIS?

texts:
I have georeferenced some images (ArcGIS -> Update Georeferencing) which creates a *.tfwx file with the new parameters for georeferencing.

I renamed it to *.tfw but QGIS is not able to read the georeferencing information.

Does anybody has any useful tips how to fix this?


texts: How to get intersected points between polygon and line in desired order

texts:
I have polygon for example "POLYGON((5 1, 6 1, 6 5, 5 5, 5 1))" that I want to intersect with line "LINESTRING(0 2, 10 2, 10 3, 0 3, 0 4, 10 4)" and get intersection points in the order that they intersect between both geometries as if you move sequentially from the start of the line to the end. Not sure this makes sense, so here is an example:

WITH impact_zone AS (
    SELECT ST_ExteriorRing(ST_GeomFromEWKT('SRID=4326;POLYGON((5 1, 6 1, 6 5, 5 5, 5 1))')) AS geom
), line AS (
    SELECT ST_GeomFromEWKT('SRID=4326;LINESTRING(0 2, 10 2, 10 3, 0 3, 0 4, 10 4)') AS geom
)
SELECT
    ST_AsText(ST_Intersection(
        (SELECT geom FROM impact_zone),
        (SELECT geom FROM line)
    )) AS iz_path_intersection_points

The result is "MULTIPOINT(5 2,5 3,5 4,6 2,6 3,6 4)"

but I want to get "MULTIPOINT(5 2,6 2,6 3,5 3,5 4,6 4)"

How can I achieve this?

I have two ideas for now that am not sure how to implement with PostGIS:

Use intersection function that will preserve the order from a geometry
Use intersection function that will preserve M attribute of the point (I have saved there the time)

Not sure if both are achievable. Any other ideas are welcome.

Thanks in advance!


texts:
I have hydrographic data in the form of an xyz file. I would like to display it in ArcMap. I need to symbolize the points based on their depth, shallow points one colour and deep points another. With a small subset of data I was able to import the symbols, symbolize based on elevation and turn on labels at an appropriate scale.

I then tried creating a multi-point file with the full dataset and set the symbology the same however I end up with unnatural bands.

Is there a way to get the multi-point file to display the same symbology as the sample set? I want the user to be able to quickly zoom in and out of the data to be able to check the depth of any given point.


texts: Adding file geodatabases as group layers?

texts:
I have a couple hundred file geodatabases (GDBs) each containing several feature classes.

Is there a way for me to move the GDBs into an mxd such that each GDB is imported as a group layer containing the feature classes as sub-layers?

I'm specifically looking for an arcpy script to do this since there are so many of the GDBs.


texts: Calculating Area Percentage of multiple overlapping polygons using ArcGIS Desktop?

texts:
I am having some trouble with a particular problem regarding area % in Arc.

I have 5 layers, each representing a different data quality, A1, A2, B, C, Overall extent. Each layer has dozens of shapefiles.

I need to calculate the area of each data quality, as a percentage of the overall area (UK territorial sea limits).

The problem is, I need to discount any lower quality layer shapefiles (e.g. B) that sits beneath an area of higher quality layer shapefiles (e.g. A1).

The problem might be better described by asking: in the top image (all layers ticked on) what is the percentage coverage of each VISIBLE layer?


texts: Python for loop with Boolean input reclassification has a save data bug

texts:
So I am creating a script that takes an input soil raster from gSSURGO (Not shown in script as this is part of a larger model). The script is below.
import sys, arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

Erosionclss = sys.argv[1]
Drainageclss = sys.argv[2]
Floodfrq = sys.argv[3]
Pondfrq = sys.argv[4]
ErodeHaz = sys.argv[5]
Data_Source = sys.argv[6]
JoinedRaster = sys.argv[7]

env.workspace = Data_Source
arcpy.env.overwriteOutput = True

input_list = [Erosionclss, Drainageclss, Floodfrq, Pondfrq, ErodeHaz]

for input_ in input_list:
    if input_ == input_list[0] and str(input_) == "true":
        ErosionclssRcll = arcpy.sa.Reclassify(JoinedRaster, "erocl", RemapValue([["'Class 1'", 2], ["'Class 2'", 3], ["'Class 3'", 4], ["'Class 4'", 5], ["'None - deposition'", 1]]))
        ErosionclssRcll.save(Data_Source + "/Erodeclss")

    elif input_ == input_list[1] and str(input_) == "true":
        DrainageRcll = arcpy.sa.Reclassify(JoinedRaster, "drainagecl", RemapValue([["'Excessively drained'", 4], ["'Moderately well drained'", 2], ["'Poorly drained'", 6], ["'Somewhat excessively drained'", 3], ["'Somewhat poorly drained'", 5], ["'Very poorly drained'", 7], ["'Well drained'", 1]]))
        DrainageRcll.save(Data_Source + "/Drainage")

    elif input_ == input_list[2] and str(input_) == "true":
        FloodRcll = arcpy.sa.Reclassify(JoinedRaster, "flodfreqcl", RemapValue([["None", 1], ["Rare", 2], ["Occasional", 3], ["Frequent", 4]]))
        FloodRcll.save(Data_Source + "/Floodfrq")

    elif input_ == input_list[3] and str(input_) == "true":
        PondRcll = arcpy.sa.Reclassify(JoinedRaster, "pondfreqcl", RemapValue([["None", 1], ["Rare", 2], ["Occasional", 3], ["Frequent", 4]]))
        PondRcll.save(Data_Source + "/Pondfrq")

    elif input_ == input_list[4] and str(input_) == "true":
        ErodehazRcll = arcpy.sa.Reclassify(JoinedRaster, "forpehrtdcp", RemapValue([["Slight", 1], ["'Not rated'", 1], ["Moderate", 2], ["Severe", 3]]))
        ErodehazRcll.save(Data_Source + "/Erodehaz")

As you can see, I take 5 different inputs (in Boolean format) and put those in a list. That list is then iterated over by a For loop to reclassify the input raster. The for loop first checks if the current iteration matches a specified location in the input_list and whether or not that value is true. If both are true, then the associated reclassification occurs. I can run each of the reclassifications individually and it will run and save in the right location with the right name, but if I run more than one, only the first input to return as true will be saved, and I can't seem to figure out where this bug is occurring.
If you try to recreate this, the data can be downloaded from the USDA Geospatial Data Gateway. you'll need to run Join Field on the MapUnitRaster_10m, comonth, component, and muaggatt. Then pass that joined raster into the script.


texts: Print composer problems Qgis 2.8 Wien with Google Satellite

texts:
I am having serious problems in exporting a map (Google satellite map) linked to QGIS.

Using the print editor I clearly see the part of the map I want to print and also in the right scale, but when I save it as jpeg, the result is:
the correct map in a circle in the middle of the page, and all around it seems that Qgis is zooming at different scale on another part of the map [see attached file].

Any suggestion about it ? I don't want to save it as simple screenshot  because I need an HD file.

]1


texts: I'm really stuck today!!  I need my "units" set in CAD to reflect a "compass".  I need the bearings I enter to reflect a complete circle with North being 0/360d, East being 90d, South being 180d and west being 270d (as shown in the attachment).  I have never had to set my units up this way in CAD.  I've always used the survey plan type bearings (ie - N 90d 10' 2" E).  I just can't seem to figure out how to set my units properly.  I haven't phrased this question very well (sorry), but can anyone help?


texts: Counting number of rows in Attribute Table using expression?

texts:
I am trying to figure out how to streamline a formula I am using in ArcMap.

I am using the "select by attributes" tool in ArcMap to select rows based on the total number of rows in the target attribute table. Right now I am using a formula that requires that I look at the number of rows in the table, and then replace the words "number of records in table" with the actual number. Since I do this quite a few times every day, I am looking for a way to avoid having to look at the attribute table and type in the number of rows in my formula.

Here is the formula I am using. How could this be modified to count the number of rows?

MOD(OBJECTID+(2*(number of records in table/20))-0,(number of records in table/20))=0


texts: Repeating points by field using ArcGIS for Desktop?

texts:
I have a point layer in ArcGIS that has a field named "Population" that stores an integer value for every point in the layer.

I want to create N points in the exact place of the feature that has N as value in the population "field". I mean, if a point has in the population field 16, I need to create 15 more points in that exact lat long (1+15 so that there are 16 in total) of that point.


texts: Returning (x,y) from axMapControl to use in zooming by ArcObjects?

texts:
i'm new to ArcObjects and i'm trying to make a program that uses an axmapcontrol and i wanted to zoom in and zoom out using the pointer , i've searched here and i've found this link which shows how to return x and y , i wanted to know what x and y do this lines of code return

How can one get the value of the point where the AxMapControl is centered at?


texts: overzoom on OpenLayers3 or a possible alternative

texts:
I'm currently trying to view a local mbtiles file of the Netherlands using OpenLayers3. Now, in itself, it's working great. I'm using a modified version of the tileserver-php framework to server the tiles into a .pbf format to then display them in Vector tiles. The result looks like this:

Which is great! However; the problem is that I need more detail than the maxzoom allows. The Maximum zoom within the metadata of the mbtiles itself is 14. Now, I know of 'overzooming', however I have not found a way to properly do this on OpenLayers. Is there something I'm overlooking that I need to set up? So far I've tried the following:

Copied the example in https://openlayers.org/en/latest/examples/mapbox-vector-tiles-advanced.html , see below for the result of that
Set up my own resolutions
Tried to have maxZoom off or on a value of greater than 14 (which results, unsurprisingly, in empty tiles)

Additionally: I return 404's on tiles that can't be loaded. Now I'm getting deeper zooms but with awful pixelation.