texts: QGIS 2.18: Configure edit widget of a editable layer with Python with "value map" from CSV

texts:
I want to configure the edit widget (in qgis: "Properties" of a layer, "Fields", column "edit widget") for a field of a editable layer (shapefile) using a Python plugin. I want set the edit widget with a "Value Map" using a CSV file for the values.

In the documentation http://qgis.org/api/2.8/classQgsVectorLayer.html#a82f6764bb1189e9d5192cc4b4285fdba i find the method setEditorWidgetV2Config, but when i try a code like this:

layToEdit.setEditorWidgetV2Config(indField, {...} )

QGIS returns this error for any configuration value in the second parameter:

(qgis.bin:6585): Gtk-CRITICAL **: IA__gtk_widget_get_direction: assertion 'GTK_IS_WIDGET (widget)' failed

what parameters should I use to configure the edit widget as a "value map" that uses the values of a CSV?


texts:
In QGIS 3.4.4, using the Field Calculator, the expression

round(-2.040,1)

returns

-1.9

where I expected it to return

-2

Converting the field to real round(to_real(-2.040),1) did not change anything.

In QGIS 2.18.19 and got the expected result.

Is the syntax changed in 3.x or is this a (known) bug? Or am I expecting the wrong output?


texts: How to select multiple values and attributes ie block and lots

texts: Is there a way to select multiple polygons that have a value for the Block attribute and a value for the lot attribute. An example is BLOCK = '100' AND LOT = '1' plus BLOCK = '103' AND LOT = '9' plus BLOCK = '200' AND LOT = '12', etc, etc...


texts: File Doesn't Exist/Isn't Supported Error - Running Python Script through C# ArcObjects

texts:
I have a script that is run using a combination of ArcObjects and a python script. The process starts in C# Arcobjects code, I create a geoprocessor, and then I run a toolbox tool (python script) with uploaded data from the user. The data is uploading successfully, I can see it in file explorer, however my python script says it doesn't exist. Now if I go in the python IDE, the script runs just fine there using the same file in the same location. I'm not sure what's going on.

This sets up the toolbox and passes the file geodatabase to the toolbox tool and runs it:

 //Set up HabitatClassification script
 IGeoProcessorResult result;
 IGeoProcessor2 gp = new GeoProcessor() as IGeoProcessor2;
 string tlbxPath = Server.MapPath("~/ArcToolbox/HabitatMon.tbx");
 gp.AddToolbox(@tlbxPath);
 IVariantArray parameters = new VarArray();
 parameters.Add(geodatabase);

 //Run Habitat Classification script
 try
 {
   result = gp.Execute("HabitatClassification", parameters, null);
 }

Now that I have passed the file geodatabase into the tool, in the python script, I set up the workspace and get the list of feature classes in the FGDB. It perform this step just fine, so clearly it can find the feature class I am looking for and it does, indeed, exist.

#Get user uploaded geodatabase and set the workspace to it
uploadedGDB = sys.argv[1]
arcpy.env.workspace = uploadedGDB

#Get the patch feature class out of the geodatabase.
#is only one feature class within it.
patchfc = arcpy.ListFeatureClasses("*Patch*")
svyPtFC = uploadedGDB + r"\{0}".format(patchfc[0])

At this point in time, svyPtFC looks like this:
C:\Users\xxx\Desktop\GIS_NET_Testing\ArcObjectsTesting\ArcObjects\UploadedFiles\NETTesting.gdb_02_17_2018_11_29_23\NETTesting.gdb\Patches

It errors out on the next line:
result = arcpy.GetCount_management(svyPtFC)

As a result, I get this displayed on my webpage (and you can clearly see it has the correct path the feature class displayed):

Traceback (most recent call last):

File "C:\Users\xxx\xxx\PythonProjects\HabitatClassification.py", line 300, in <module>
result = arcpy.GetCount_management(svyPtFC)

File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 17419, in GetCount
raise e

ExecuteError: Failed to execute. Parameters are not valid.

ERROR 000732: Input Rows: Dataset C:\Users\xxx\Desktop\GIS_NET_Testing\ArcObjectsTesting\ArcObjects\UploadedFiles\NETTesting.gdb_02_17_2018_11_29_23\NETTesting.gdb\Patches does not exist or is not supported

Failed to execute (GetCount).

Failed to execute (HabitatClassification).


texts: How to read the count of the cells with data with non "no-data" values ,in QGis?


texts: Issue with using minimum bounding geometry tool with points in QGIS

texts:
I am attempting to create a 100% MCP in QGIS using the minimum bounding geometry tool in QGIS. My data is sensitive and I therefore cannot share any of it (including a subset). I have used a dataset of animal locations as the input layer, left the field open as I need an MCP for all individuals together, and have chosen the "convex hull" geometry type. However, the output clearly ignores some of my location points (see left-hand side of the photo).
I have tried redoing this multiple times, even re-importing the point dataset, yet I still get the same result. I also tried creating a 100% MCP in R and then importing a shapefile of it into QGIS, and it still looks the same. Has anyone ever had this issue come up? I cannot seem to find anything wrong with the few points it's ignoring, so I am not sure what else it could be.


texts: Getting underlying layer data in entry form

texts:
We have a requirement to place a point feature on a line feature (a road) in QGIS.

When placing the point, we would like the road number from the line feature to be automatically populated in the data entry form of the point feature attributes.

We are using out of the box features of QGIS to connect to WFS/WMS and MSSQL spatial. We have tried the spatial analysis features of qgis, but they always expect the results the be saved as an external .shp file, and not back to the database. We have also tried varous field calculators to make a spatial query on the fly from the roads layer but not successful. We would like to know if making spatial query on adding new point is the generally accepted practice, or like DPSSpatial suggested that we do it after the fact using triggers.


texts: Why do collinear ground control points create larger residuals?

texts:
When georeferencing a raster image of a long corridor, why do collinear ground control points cause larger residuals?

Say I am using the georeference function in ArcMap, and I have 6 ground control points on a site that is is shaped like an L. Three ground control points are close to collinear on the long end. I find that if I use just 4 of the points (exclude the 2 of the 3 collinear points) I get a lower root mean square.

Why are collinear points giving me a worse root mean square?


texts: Checking/unchecking rule-visibility of layers at once in QGIS

texts:
I've got several rule-based layers in my project (qgis 2.18.13).
Is it possible to check/uncheck the visibility of all rules (not the visibility of the layer itself) at once?
If yes - is it possible to do this on all layers of the project at once too?


texts: Use folder with multiple rasters as an input for the r.series tool

texts:
I need to use the r.series tool (that makes average raster out of multiple rasters) in a model. the idea is that each time the user will specify a folder with rasters, and the r.series tool will output a single raster that is the average of all the rasters in the folder (the number of rasters is changing for each folder).

When I use the modeler I can insert raster as parameter but only one-by-one. I know that it is possible to use "folder" as parameter but I cannot use it as an input for the r.series tool.

So.. is there a way to use a folder as an input for the r.series tool? or would you recommend a different way to accomplish that?


texts:
I am using QGIS 2.18.10. I am working in UTM Zone 30 WGS 84.

I have created a vector grid 1000km X 1000 km using the Vector>Research Tools> Create Grid

I have another shapefile which has point data in it.

I am trying to use: Vector>Counts point in polygons, and selecting the correct dataset.

However it is giving me an error: "unknown see log for more details". But that is in the log. So I can't figure out why.

What might I be doing wrong?


texts:
I am trying to use the Road Graph plugin within QGIS 1.8 (running under Windows 7 64-bit). When I click Vector -> Road Graph -> Settings I get the settings-panel:

But, I don't know how to enable the "working-panel", which should appear on the left side of the main window, according to the producers of the plugin (screenshot of pic #2 is taken from their site):

Any hints? What am I doing wrong / what's missing?


texts: Getting the results of a PostGIS query to a file with ogr2ogr with low memory usage

texts:
I have a PostGIS query returning geographies and attributes. The result set is deterministic and ordered (so OFFSET and LIMIT work for getting subsets).
The result set of my query would exceed my RAM by many times.
ogr2ogr can read directly from PostGIS and load SQL from a file so I built this command line:
ogr2ogr -f "ESRI Shapefile" output.shp PG:"host=host user=user dbname=gis" -sql \@q.sql

When I run this, my memory is exhausted and the process is killed by the system.
How can I make PostGIS+ogr2ogr "stream" the results and flush them in reasonable intervals instead of trying to get the full result set at once?
I am not looking for scripting, instead I wonder if there is some magic of GDAL I just have to enable to get paging.
A different output format would be fine if that makes it work. But only if I can convert that to Shapefile or GeoPackage without exhausting RAM again.


texts: Find polygons that overlaps in the same layer with SQL

texts:
I have one layer with 40k polygons and I would like to create new layer only with the polygons that overlap each other.
I have tried to use the SQL in the db manager in order to create some separation but seems like I cannot run the query on one layer.
this is how I tried to do that:
SELECT "45k_less".*,
ST_Overlaps("45k_less"."geometry","45k_less"."geometry") "overlaps"
FROM "45k_less"

I have also tried the same with contains but seems like I miss something here -- maybe because it's the same layer?
I saw there is way to do it with the topology checker but is there any way to do this query inside the db manager?


texts:
I have a polygon shapefile (created by someone else) that is a bunch of little watersheds. Unfortunately, it doesn't include any coordinates in the attribute table for the lat/long of the pour point of said watersheds.
Is there a way to batch process roughly 1300 polygons to find the coordinates of the pour point of each polygon or is this just not possible?


texts: MRR -setting transparent color (MapInfo)

texts:
In Mapinfo Advanced 16

I tried to merge all the 2016/10cm, 2014/15cm and 2011/50cm separately into mrr (works fine) and then to merge these three together. It started to process but (3 days) later MI crashed. Even if it did work I am afraid I would get the no imagery areas on the higher rez file (higher resolution doesn't cover the whole area) as white over the lower rez (covering the whole area).

What can I do to at least view these three datasets in mapinfo without the white -the usual way of setting the transparent color doesn't work on the mrr's.

Properties in MRR file

10cm imagery over 50cm

50cm imagery without 10cm imagery

disregard the white no-data area, that tile was missing from our project.

What can I do -the end result should be a highest rez possible mrr or a mapinfo workspace with each layer seperately not showing the areas with no data. This is for use in SSA so I am not (to) worried by the file size as it will get tiled.


texts:
I know that a lot of people asked already that question, I read all of the answers and have also seen all the videos, but I still don’t understand how to link a pdf file to an object.

I have a huge database with the biotopes of a region.  Now I want to add to a few biotopes a pdf file (Similar to How to link a pdf file to a specific feature with QGIS?). For that, I created a new column named “Action” and in that column I wrote the exact path where the pdf file is stored on my computer (but just for biotopes I need with a pdf file).
Now if I click on one biotope it just opens the attribute list of that shapefile and I can see the path to the pdf file, but it doesn’t open the pdf

In the properties I also wrote down the path, where the pdf files are stored, but it still doesn't work and I don’t really understand what is wrong?


texts: How to access the metadata link URL field of GeoServer using the OpenLayers Library?

texts:
GeoServer v: 2.4.4

OpenLayers v: 2.13

All maps are WMS.

In my JavaScript I use "layer.name" to access the Title of the specific map as specified in my GeoServer layer. Using jQuery I was able to inject the "layer.name" of the selected map in the footer of my HTML. Now I would like to add a hyper-reference of the metadata.

In GeoServer, I have specified the URL of my metadata under the header: "Metadata Links" in the field URL. (This all seems very logical and simple).

The issue I am having, is that I don't know how to access this URL... I've searched the documentation, but can't find any syntax that would allow me to easily access this field.

I would like to avoid typing in the URL manually, because I have a whole list of maps, and I have added the metadata for each map in GeoNetwork.

Does anyone know how to get the "Metadata link URL" field with OpenLayers? Or have a nice workaround?


texts:
I would like to reproject images in ENVI classic. The purpose of this is so that I can view multiple scenes in a uniform direction so that I can create an ROI and observe the differences in the ROI in each scene. Is there a way that I can reproject my images to make sure all of my scenes are true north?
The Harris website gives a tutorial on how to reproject in Envi, but not Envi Classic.


texts: Avoid MapBasic crash when unable to open file

texts:
I'm working on some MapBasic code and in the process I want it to open the .tab files in a given folder. But when it reaches a .tab file that it is not able to open (e.g. attached to a raster that does not exist, ID/DAT/MAP files missing, etc.), I get the following warning message and it won't open the files that come after it:

Snppet of code:

If FileNames(i) like "%TAB" then
    Open Table FileNames(i)
End If

Is there away to have the code continue after this type of error is encountered so that it will load all of the non-corrupted files? Or better yet have the code load all the files and give a list/report at the end of which files it was unable to open?


texts:
I am beginner of QGIS.

The purpose of using QGIS is making contour line and 3d visualization.
I practice to make contour with DEM and it is not hard, but in my work, I have to treat contour line on paper maps.

I want to find line tracing tool or line detection tool to trace contour on scanned paper map. it is binary image and, unfortunately, contour line is not closed. Until now I make shapefile layer (line) and draw contour manually, of course, it is pretty time consuming.

I hope to find some method like, when I click some point, trace the line until meeting intersection or disconnection. then draw line manually, and restart line trace.

Do you know good method to do that?

I need your help.


texts:
Is it possible and/or feasible for an Oracle trigger to periodically update ArcSDE via a REST geoprocessing service?

Does anyone have any experience with such a thing?

The reason for this is that I have an Oracle database that is separate from the SDE database, and there are number of business rules that we want to run first, before inserting the geometry into SDE.

So the basic workflow would be:
  1. Point features via ArcGIS Server via ArcSDE loaded onto web map
  2. User moves point or edits attributes
  3. Edits to points first get written to Oracle, then, via trigger to ArcSDE


texts: How can I have e.g. all highway speed cameras displayed on the map, i.e. what is the URL syntax to search OpenStreetMap for certain tags? And where is this information (tagged locations, POIs etc.) coming from, is it collected/contributed by OSM community members or does the OSM project use third-party sources (private/public/proprietary)?


texts: How to import a Composer Template into QGIS Print Composer?

texts: How can I import a .qpt file into Print Composer (Composer Manager) in QGIS. I've tried saving it in the composer_templates folder but to no avail.


texts: ArcGIS pythonaddins.ProgressDialog inconsistent progress bar drawing

texts:
I am trying to use pythonaddins.ProgressDialog in a Python add-in. I am on ArcGIS 10.4.

I will need to perform several steps that will take some time. To show how much work is done, I plan to use a progress bar which is available with the  pythonaddins.ProgressDialog.

Trying to run this sample code in the Python window in ArcMap (taken from the Esri help page):

import pythonaddins
with pythonaddins.ProgressDialog as dialog:
    dialog.title = "Progress Dialog"
    dialog.description = "Copying a large feature class."
    dialog.animation = "File"
    for i in xrange(100):
        dialog.progress = i
        time.sleep(0.125)
        if dialog.cancelled:
            raise Exception("Ooops")

The code works fine unless you increase the step in the xrange() from 1 to something bigger. For instance, if you have 10 feature class to copy, it would make sense to make a step of 10 units (100 / 10) so every increase in the bar would correspond to a copied item.

However, I am observing inconsistent drawing of the progress bar. Try running this code in Python window in ArcMap, in a Python add-in or in a script tool.

with pythonaddins.ProgressDialog as dialog:
    dialog.title = "Progress Dialog"
    dialog.description = "Copying a large feature class."
    dialog.animation = "File"
    for i in xrange(0,100,10):
        dialog.progress = i
        time.sleep(0.5)
        if dialog.cancelled:
            raise Exception("Ooops")

This is how it looks on my machine. As you can see, the progress bar is jumping back and forth. I've tried multiple steps (2,5,10,20) and observed the same behavior: as long as the step is not 1, the bar jumps back and forth in progress.

Has anyone seen this and know if there is any workaround for this?


texts:
I'm working with shapefiles of water protection areas (of Germany) to find out weather or not somme geocoded adresses are inside of these areas.

Every federal state delivers shapefiles in a different CRS and I usually solve this by choosing "Save As" to project them in my Project CRS (EPSG 4326) before further processing them.

The file for "Baden Würtemberg" however aborts with an Error message when I try to do this.

Link to the source of the shapefile:
http://udo.lubw.baden-wuerttemberg.de/public/pages/url/show.xhtml?url=http%3A%2F%2Frips-dienste.lubw.baden-wuerttemberg.de%2Frips%2Fripsservices%2Fapps%2Fgeodatenexport%2Fudo%2Fdownload.aspx%3Fid%3D32

Error message when "Save As" is used:

  Export in Vektordatei schlug fehl. Fehler: Objektschreibfehler:
  Objektgeometrie nicht importiert (OGR-Fehler: ) Objektgeometrie nicht
  ...
  importiert (OGR-Fehler: )Abbruch nach 1001 Fehlern

translated:

  Export to shapefile failed. Error: object writing error:
  Object geometry not imported (OGR-Error: )
  abort on error number 1001

What I tried so far:

I can "Save As" other shapefiles from and to the chosen directory. I also tried creating a new QGIS project with the shapefile directly on the desktop (so its not a permissen problem).
"Save as" doesn't work when not changing the CRS as well
I checked with "validate geometry" and eliminated all the errors by deleting duplicate nodes and overlapping areas of the shapefile. Different shapefiles with non-valid geometry can be saved too.
Save to a file without critical signs like spaces, umlaute (ä,ö,ü,ß).

Any ideas where this error comes from?

Am I missing out some problem with the original shapefile?

Maybe you could try to reproduce this to find out if the problem is in the original shapefile itself or caused by my GIS-configuration.

Using QGIS-2.10.1 on Win7. I'am stuck with these specs for now, will try with newer QGIS version and report back later.

Back home it works with QGiS2.8.6. Maybe there is an bug only in 2.10.1 or, more likely, I missed something.


texts: Google Geocoding API to geometry (Multi-line string)

texts:
I am using the Google Geocoding API. I need to use that response and convert that to geometry. How is this possible? What postgis function should I use and what parameters to assign?
Raw JSON data:
"geometry":{
   "bounds":{
      "northeast":{
         "lat":35.1614839,
         "lng":33.3903832
      },
      "southwest":{
         "lat":35.1595564,
         "lng":33.3889372
      }
   },
   "location":{
      "lat":35.1602836,
      "lng":33.389555
   },
   "location_type":"GEOMETRIC_CENTER",
   "viewport":{
      "northeast":{
         "lat":35.1618691302915,
         "lng":33.3910091802915
      },
      "southwest":{
         "lat":35.1591711697085,
         "lng":33.3883112197085
      }
   }
}


texts: Why TreeView Kendo is not executed properly?

texts:
I made a tree view kendo to enable the user to turn off or turn on the two layers "osmlayer" and "bing_road". but the treeview is just shown and does not work. first of all chrome debugger showed that "length" property is not defined, so I changed map.layers.length to map.length. now it says that getLayersByName is not a function. how can I make this work?

here is the code:

  Create_Legend();

var selected_layers = [];
var selected_layers_name = [];
function CheckedNodeIds(nodes) {
    for (var i = 0; i < nodes.length; i++) {
        layer_id = nodes[i].id;
        if (nodes[i].hasChildren) {
            CheckedNodeIds(nodes[i].children.view());
        }
        else {
            if (nodes[i].checked) {
                map.getLayersByName(layer_id)[0].setVisibility(true);
            }
        }
    }
}
function onLegendCheck() {
    var treeView = $("#panelLayers").data("kendoTreeView");
    for (i = 0; i < map.length; i++) {
        map[i].setVisibility(false);
    }
    selected_layers = [];
    selected_layers_name = [];
    CheckedNodeIds(treeView.dataSource.view());
}

function Create_Legend() {
    $("#panelLayers").css("background-image", "linear-gradient(rgba(43, 42, 42, 0.5),rgba(97, 176, 190, 0.5)");
    //$("#panelLayers").addClass("card");
    $("#panelLayers").kendoTreeView({
        checkboxes: {
            checkChildren: true
        },
        check: onLegendCheck,
        dataSource: [
            {
                id: 1, text: "نقشه های پایه", expanded: true, items: [
                    { 'id': 'osmlayer', 'text': 'OpenStreetMap', "checked": true },
                    { 'id': 'bing_road', 'text': 'Bing Map' },
                ]
            }
        ]
    });

    onLegendCheck();


texts: Automatically select a feature from another layer based on shared attribute in QGIS

texts:
I have two shapefiles, A (polygons) and B (points), and I need to select all features from shapefile B that have the same attribute (for example, "Name") as shapefile A, and I need to do this iteratively for a large number of features.

Normally I would browse the table of shapefile A using the Go2NextFeature plugin, select a record, copy its name and make a "select by attribute" selection on layer B.

It would be best if whenever I select a new feature in shapefile A, all features of shapefile B with the same name would automatically be selected.


texts: I am creating a tool that opens an IGxDialog which it's objects are filtered to show only workspaces, my question is how can I convert IGxObject to IWorkSpace or get the selected workspace without taking the selected object name and opening it using an IWorkSpaceFactory


texts: List all available extensions using Python to get extension code

texts:
I am building an extension which depends on the availability of another extension (not from Esri) so I'd like to find out if another extension is available using Python.

Based on the availability of this other extension, I want to set the 'Enable' property of toolbar icons.

I tried to use the arcpy.CheckExtension('') however I am getting this message
*Unrecognized extension code'. This is because I don't know the extension code.

Is there a way to list all the available or enabled extensions with their extension code?

Edited:
From the Customize > Extensions.. I can see the extension is enabled (checked) but from python I like to find if it is enabled/checked. I tried arcpy.CheckExtension('') with displayed named but no luck.


texts: Rasterize points in a map : result is not leasible. How to make the raster (points) bigger?

texts:
I have a vector layer (points) that I rasterized based on one attribute of the vector layer. I chose the GRASS GIS 7 command : v.to.rast.attribute. The raster layer appears in the map, but is not very readable.

Here is the original vector layer (points) :

Here is the rasterized result (based on one attribute of the layer) :

I would like to make the points bigger in the raster.


texts: Get Extent of Multiple Shapefiles (Batch Process)

texts:
I have 1000+ Shapefiles and want to get extent of these polygon shapefiles.
I know the manual way to do so.. Is there any way to do it through Script .


texts:
I am trying to convert an existing tile cache (created by custom process) into a usable format/structure for a GeoWebCache WMTS layer.

I can't seem to figure out the directory/naming structure required - using a layer created and seeded by GeoWebCache as a reference.  The first level is easy to understand EPSG_4326_00 - that would be spatial reference and tile matrix level.

The next level of directories, and the names of the tiles themselves is totally alien to me.  If anyone could explain - or point me at some documentation - that would be greatly appreciated.


texts: Using WCS DEM data for make slope and aspect maps in QGIS

texts:
I need to do several slope and aspect analyses (hundreds).
Instead of downloading an ASCII for each location, I want to extract/clip the elevation data from wcs and do the analysis.
This is the WCS service (DEM of Spain):
https://servicios.idee.es/wcs-inspire/mdt?version=2.0.1&request=GetCapabilities&service=WCS
Is there any way to do this on QGIS?


texts: Using map algebra get minumum and maximum value in ArcMap

texts:
I am looking for an easy solutions with Map Algebra for this one without Saving:

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/temp"

# Set local variables
inRaster = "DOM.tif"
inMask = "Mask.shp"
maskBuffer = "C:/temp/maskBuffer.shp"
temp = "C:/temp/temp.tif"

arcpy.Buffer_analysis(inMask, maskBuffer, "5 Meters")

outExtractByMask = ExtractByMask(inRaster, maskBuffer)
#is it possible to get the Minumum and Maximum directly from this Map Algebra?
outExtractByMask.save(temp)

tempRaster = Raster(outExtractByMask)
print(tempRaster.minumum)

This is a try to get the height of buildings in a loop, so I will ExtractByMask for each Feature in maskBuffer.


texts: Extra lines in polygon in ArcScene vs ArcMap

texts:
I have a polygon of a park boundary in ArcMap. It is a single line around the perimeter of the park. When I display that same polygon in ArcScene, it looks the same. However, when I drape that layer on the DEM (using "float layer on a custom surface"), all sorts of extra lines appear in the polygon. I have clipped the screen below to show a) the image straight out of ArcMap with the polygon on a DEM layer,

a)

b) the polygon in ArcScene without adding Base Height information

and c) the polygon in ArcScene once I have added "Base Height" information

What is going on with this?


texts: How to create random points in a polygon in PostGIS

texts:
I have a polygon layer, and would like to generate random points in each polygon. An initial searching in Google resulted this website which presents a function for this purpose:

http://trac.osgeo.org/postgis/wiki/UserWikiRandomPoint

The function I'm interested to use is called RandomPointsInPolygon()

When I use this function and execute my query, it doesn't give an error but also never finishes the process. I tried giving a WHERE clause and limiting my query, this time process executed without any error but no points were stored in the relevant field in my table.

Does anybody know what the problem is? and/or what could be an easy way to solve my problem? using this function or any other options that you may have know?

The code for function:

CREATE OR REPLACE FUNCTION RandomPointsInPolygon(geom geometry, num_points integer)
  RETURNS SETOF geometry AS
$BODY$DECLARE
  target_proportion numeric;
  n_ret integer := 0;
  loops integer := 0;
  x_min float8;
  y_min float8;
  x_max float8;
  y_max float8;
  srid integer;
  rpoint geometry;
BEGIN
  -- Get envelope and SRID of source polygon
  SELECT ST_XMin(geom), ST_YMin(geom), ST_XMax(geom), ST_YMax(geom), ST_SRID(geom)
    INTO x_min, y_min, x_max, y_max, srid;
  -- Get the area proportion of envelope size to determine if a
  -- result can be returned in a reasonable amount of time
  SELECT ST_Area(geom)/ST_Area(ST_Envelope(geom)) INTO target_proportion;
  RAISE DEBUG 'geom: SRID %, NumGeometries %, NPoints %, area proportion within envelope %',
                srid, ST_NumGeometries(geom), ST_NPoints(geom),
                round(100.0*target_proportion, 2) || '%';
  IF target_proportion < 0.0001 THEN
    RAISE EXCEPTION 'Target area proportion of geometry is too low (%)',
                    100.0*target_proportion || '%';
  END IF;
  RAISE DEBUG 'bounds: % % % %', x_min, y_min, x_max, y_max;

  WHILE n_ret < num_points LOOP
    loops := loops + 1;
    SELECT ST_SetSRID(ST_MakePoint(random()*(x_max - x_min) + x_min,
                                   random()*(y_max - y_min) + y_min),
                      srid) INTO rpoint;
    IF ST_Contains(geom, rpoint) THEN
      n_ret := n_ret + 1;
      RETURN NEXT rpoint;
    END IF;
  END LOOP;
  RAISE DEBUG 'determined in % loops (% efficiency)', loops, round(100.0*num_points/loops, 2) || '%';
END$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 1000;
ALTER FUNCTION RandomPointsInPolygon(geometry, integer) OWNER TO postgres;

My query:

  INSERT INTO grid(point_geom) SELECT RandomPointsInPolygon(h.geom, 1)
  FROM grid AS h;

Here I have a point geometry column column named point_geom to store the result of query, and in the same table I have polygons stored in the geom field. I would like to generate points inside each polygon of this grid table. In the query above I tried to create 1 random point for each polgyon and store it in the relevant field.


texts: Label disappear after setting data defined position

texts:
Every time I make a map I like to have more control over labels position so I manually move them. But today something weird occurred. I defined X and Y position of label to manually move them. Every time it was working but this time labels just disappear after saving the layer settings. What's wrong here ? The "Move Label" option is active but there is no feature to move...

I Use Qgis 2.14.15 on Windows 10 HOME.


texts: Getting attributes of OpenLayers map with JavaScript

texts:
I have created an OpenLayers map via the following code:

var map = new ol.Map({
target: 'mapdiv',
layers: [
new ol.layer.Tile({
  source: new ol.source.MapQuest({layer: 'sat'})
  }),
new ol.layer.Tile({
  source: new ol.source.MapQuest({layer: 'hyb'}),
  opacity:01
  })
],
view: new ol.View({
  center: ol.proj.transform([7.72,47.94], 'EPSG:4326', 'EPSG:3857'),
    zoom: 10
    })
});
console.log('map created');

Now, of course, 'mapdiv' is the ID of the div I want the map to show up. This part is working according to the plan.

I nonetheless am implementing a separate JavaScript library to manipulate maps. As such, I want to access the created map (post creation) and its information. I look at the built in DOM inspector in Firefox, to see what exactly is created inside my div called 'mapdiv'.

The direct child of my 'mapdiv' is : <div style="position: relative; overflow: hidden; width: 100%; height: 100%;" class="ol-viewport">, which in turn has the child : <canvas height="696" width="1392" class="ol-unselectable" style="width: 100%; height: 100%;"></canvas>

This canvas seems to be the relevant one here. In the Firefox element inspector console, hovering my mouse over this generates a small preview of the map.

See image:

So I would like to grab the openlayer.Map object responsible for this, and call the calculateBounds method, to know the latitude and longitude boundaries of the displayed map.

Noting that the map is displayed in an element given by class, I am aware of this. But note that I am not aware of the class name beforehand. If for example, there are multiple maps, then using the linked method won't work, because they all have the same class, despite residing inside different divs.

I thus need a method to access the OpenLayers maps in a div element via JavaScript, if the div ID is known.

I am using OpenLayers 3.


texts:
Is there a managed library for reading the coordinates from a geotiff file in .NET ?

I know it can be done thru GDAL, but I am looking for a managed alternative.


texts: QGIS and satellite images band combination

texts:
When I add a satellite images as a raster in QGIS:

Why the image is changed, it looks like it has only one band?
Where is a tool for band combination?


texts: Copying features from one layer into different layer using ArcGIS for Desktop?

texts:
I would like to have two different kinds of feature (lets say coniferous and deciduous woodland).

To represent this I have two different layers that are different colours.

In an edit session when I copy features from one layer using the standard toolbar, I have a choice of pasting the features into either layer. The problem is that the feature always pastes into the same layer, no matter which one I choose. I can only seem to create features in one layer as well, even though there are two different templates available in 'create features'.

Can someone explain what I am doing wrong? I was managing to produce features with different appearance ok two days ago, I can't work out what I am doing differently now.


texts: Fitting rectangle with maximum area inside polygon in PostGIS

texts:
I am looking to fit a rectangle inside a polygon which would have the maximum area for the rectangle.
I tried using ST_OrientedEnvelope(geom) but this would give the rectangle fitted outside the polygon as shown in the image below:

Ideally what I am looking to fit something as below:

coordinates of the sample polygon
POLYGON((529004.994 159096.691,529005.519 159096.992,529011.013 159100.145,529021.878 159106.109,529022.1 159105.95,529033.4 159097.85,529066.85 159073.7,529059.1 159060.25,529005.5 159096.35,529004.994 159096.691))
What's the best way to achieve this?


texts: Batch Share as Web Layer in Arcgis Pro 2.3

texts:
I have a federated Arcgis Portal 10.6 and I use Arcgis Pro 2.3 for publishing.

What I'm looking for

My final aim is to publish all the layers in a given map individually (not as a single Map Service with all layer at the same time). I want to be able to use individually each layer and it's cache.

What I do now

When I wish to publish a layer from Arcgis Pro, I go in the map where all my layers are and right click on the layer I wish to publish.

Then I select Sharing | Share as Web Layer :

I then have the following panel :

My Question

How do I batch this single operation ? I couldn't find a toobox matching this process. I did find a page mentionning arcpy.sharing but I would like to avoid coding if possible.


texts:
Iam struggling with one letter below the other vertical text orientation for DivIcon. Found some hints/hacks here https://davidwalsh.name/css-vertical-text

But i couldnt manage to get it working for me. My result was that the letters just lay above each other. Is that even supported for DivIcon? Another possibilty would be to create a DivIcon for each letter. But thats very unusual i think.


texts: How to get parcel boundaries with tax roll, land use, and zoning info?

texts:
How do I get parcel boundary shapefiles with tax roll, land use and zoning info included?

Is it possible to just create the parcel boundaries data from scratch, i.e. could I create parcel polygons in QGIS?

I would imagine there is a somewhat efficient way to do it considering companies like Corelogic and Digital Map Products have nationwide coverage.


texts:
I am working with QGIS on a Mac, and understand the process for importing a CSV file. However, I do not know how to create the CSVT file to define the CSV file's data types.
Can someone instruct me on how to do this?


texts: QGIS 3x Update existing field autoselection in field calculator

texts:
I want to have the "Distance [meters]" column autoselected when launching a data field calculator.
It's driven by a multitude of records, which I want to put by using the expression (individual for every ID). As I have several records, the job is quite tedious. I have to select the field for updates every time.

Is there any way to set the target field for updates and make it autoloaded when running a field calculator?
I think, some hint can be here:
Is it possible to have auto-updating area fields?
but it looks like there is no Autofield plugin for QGIS 3x
Functionalities of expression builder in Autofield plugin available in QGIS 3.4?
How can I have the field auto populated in the data attribute table calculator?


texts: How you can repair invalid topology errors within postgis?

texts:
I've imported the shapefile of corine land cover project to my postgis database and the

select gid, ST_IsValidReason(geom) from gis.clc2000 where not ST_IsValid(geom);

returns in every occasion

gid    |                                st_isvalidreason
-------+--------------------------------------------------------------------------------
     1 | IllegalArgumentException: Points of LinearRing do not form a closed linestring
.....

I know that my polygons are closed by visually instecting some of the errors.

How can I fix the topologies errors from withing the server?

edit: I'm using postgresql/postgis2.0.1, Ubuntu12.04


texts: High resolution and reasonable cost RGB satellite images

texts: I need to obtain high-resolution RGB satellite images with reasonable cost for Pakistan. The objective is to guide tourist on high-resolution GIS maps. Google Earth does not have very high resolution images for certain areas in Pakistan. Which satellites can provide very high-resolution images for this area?


texts: Managing large number of catchments/watersheds?

texts:
Does anyone know of an easy way to manage a large number (10,000's) of catchments/watersheds?

I have defined a stream network and for each segment, I have a corresponding catchment (in a unique shapefile).  The name of the catchment is the x and y co-ordinate of the catchment outlet which is also assigned to the attribute table of relevant stream segment.

The end user will only want to view the stream network and maybe a few of the catchments (certainly not all!) at a time.  I need a way to make it easy for the end user to view the catchments of interest (I guess they could search for them).

My only idea (which I am a little unsure of how to implement) is to build a tool (in ArcGIS) that loads the catchment when they click on a stream segment (with the tool selected)

For my idea, I think I would have a button similar to the identify tool the end user would use this tool to select a stream segment, the tool would then look in the attribute table to find the file location of the catchment and add this as a shape file.

Is such a thing possible to do?

I have coded up some python scripts before but have never attempted to integrate such functionality into ArcGIS.

We only have the basic version of ArcGIS, although if this sort of functionality is available in something like Arc Hydro I'd be interested to know!


texts: Avoid Crash Minidump while using Semi-Automatic Classification Plugin

texts:
I'm using Semi-Automatic Classification Plugin v4.9.5 on QGIS 2.14.3.
OS is Windows7 64 bit.

After I conduct some analysis by Semi-Automatic Classification Plugin, it is always shown "Crash minidump".

Would you tell me how to avoid this?

I've checked "QGIS 2.14. minidump crash when closing program" and so far.
Uninstall the QGIS and remove .qgis2 folder.
And reinstall QGIS and Semi-Automatic Classification Plugin.

But the message is still shown.

Regards,
Akihiko


texts: Real time positioning of vector features using database

texts:
I am a newbie in GIS and want to implement something like a vehicle tracker wherein I will get a series of coordinates in a database.

I have no idea how to go about it, I mean what should my process consist of? I have found something about GeoJSON, but since I'm using GWT should I be really using it? I guess I can directly interact with the Java data structure.

Apart from that, how to initiate the loop for animation? I have seen rotate and move methods in openlayers, but would they suffice in case of a complex motion like vehicular movement.


texts:
I have several hundred points that represent samples in a landscape:

As you can see, the points are unevenly spread and form large, dense clusters. To account for spatial autocorrelation, I want to subset or thin these points so that the points that remain are a minimum distance from each other. At the moment I've set up a pairwise distance matrix for them and am randomly removing points that are below the distance threshold until the only ones remaining meet my minimum distance criteria. However, this does not seem optimal from two different perspectives: first, different potential random subsets using this method can return different numbers of samples. Ideally, I'd retain as many points as possible. Second, different subsets may be distributed differently across the landscape. Ideally, I want a point subset with a distribution that more uniformly represents the landscape. Are there any recommended or general algorithms for optimally subsetting points (at least more optimal than my approach)?


texts: WMS GetFeatureInfo not working at larger scales

texts:
I use OpenLayers/GeoServer on an Ubuntu server. GeoServer runs on Tomcat6. I have a grouplayer of several layers all shown at certain scales.  When I click at feature at small scale, the GetFeatureInfo request works. When I zoom in and then click at the same feature, the response from server is empty.

There are no errors in the GeoServer logs or Apache logs. Also when I send that request manually via web browser, the response is empty too. At one layer, there are even two zoom levels working, only the last one is empty. BBox and pixels in request look OK.

 >>> zoom >>>

Most probably, the problem isn't in code, but here it is:

info = new OpenLayers.Control.WMSGetFeatureInfo({
        url: 'http://158.196.109.46:8080/geoserver/diplomka/wms',
        title: 'Identify features by clicking',
        queryVisible: true,
        maxFeatures: 1,
        eventListeners: {
            getfeatureinfo: function(event) {
                map.addPopup(new OpenLayers.Popup.FramedCloud(
                    "chicken",
                    map.getLonLatFromPixel(event.xy),
                    null,
                    event.text,
                    null,
                    true
                ));
            }
        }
    });
    map.addControl(info);
    info.activate();

Proxy.cgi is also set properly I think. I also tried raising Tomcat heap memory, it didn't help. So I really don't know, what should be causing this.

Anybody experienced this, or knows what to do?


texts:
What is the best way to change the cell size of a netcdf file, keeping it as netcdf?

I have looked into ArcGIS and GDAL and have not found anything yet. I prefer to use python if a solution is available.

The netcdf file is CRU monthly rainfall data, which I am running zonal statistics on to get mean value for districts in Malawi. Because the cell size is large, districts are showing up as having no data. This makes sense because the center of the cell is not located within the polygon districts with no data. My goal is to change the resolution from ~60km2 to ~10km2 which should pick up the cell values in the districts when I do zonal stats.


texts: OpenLayers2 and GeoJSON - possible to create, modify and delete features?

texts:
I wondered if anyone has some OpenLayers GeoJSON code to read, create, modify and delete features (i.e. instead of using WFS-T)?

This is how far I've come - I've no problem changing the pop-up to a HTML form to modify attributes, it's more the linking the create button to a pop-up to enter data and then catch the move feature event, and last a 'delete feature' button

  <title>OpenLayers Example</title>


texts: Make square polygon from one point and area

texts:
I have a shapefile of points with a variety of attributes, the relevant ones being "ID", "area", "lat", and "long". Such as...

The lat/long coordinates are always the SW corner of a plot and the area varies from 10 to 50 km-sq. Each plot is always a square. I would like to make a polygon shapefile from the point shapefile that displays the entire plot to be able to clip other features. To look similar to...

This is a one-time operation in this case and I would like to do it in Arc (I'm using ArcMap 10.2.2) but would be willing to use R if it is easier. Maybe it would be easiest to calculate the coordinates of each corner first? I have found answers similar to my question in:

Createing polygon layer from point data using QGIS?,
Creating polygons (squares) from point using PostGIS?, and
Calculating coordinates of square x miles from center point?,

but not quite what I want and using different programs.


texts: Specific QGIS file will not open and "QGIS quit unexpectedly" screen pops up

texts:
This problem has occurred for me several times for different files and I can't figure out what the cause is. I will be working on a QGIS file for several weeks and then all of a sudden one day I can no longer open the file and the attached screen appears. I don't think the cause is missing shapefiles with broken links.
I've done some research and think I maybe need to open the file in a text editor, but I have no experience with this and don't know what part of the text I'm supposed to alter. I have been using QGIS 3.0 on a mac.


texts:
I created a trigger that checks if geometry is valid on insert and update.

        create or replace function check_geom() returns trigger as
$$
begin
--insert a polygon
IF st_isvalid(new.geom) = 'f'
THEN RAISE EXCEPTION 'erreur de topologie' ;
END IF ;
RETURN NULL ;
end;
$$ LANGUAGE plpgsql;
DROP TRIGGER IF EXISTS trg_check_geom ON public.uniteobservation ;
CREATE TRIGGER trg_check_geom AFTER INSERT OR UPDATE on public.uniteobservation
FOR EACH ROW EXECUTE PROCEDURE check_geom();

I have an error message when I create a feature with a wrong geometry (self intersection,...). I can correct and validate it.
But I can't update a second times : QGIS mini dump each time I select a vertex again.
Any idea ?


texts: QGIS count how many lines cross a polygon

texts:
In QGIS I create some lines that cross three polygons. By using intersection, it will create a new shapefile that gives me as a result the lines inside the three polygons.

If I open the attribute of the intersection layer I can see that different line ID have same polygon ID (obviously because multiple lines cross the same polygon).  By looking I can tell how many lines are inside a polygon but how can I do this automatically?

I would like to modify the polygon layer attributes by adding a new column with the number of lines crossing each polygon (so the polygon layer will have only two column: "ID" and "number of lines").


texts: Watershed-segmentation with orfeo toolbox

texts:
I work at the moment with the watershed-segmentation from the orfeo-toolbox. After the segmentation the pictures is split into squares. Specially the square in the middle seems to be deeper the the rest of the picture, which doesn't confirm with the real topographic situation.

Does anyone have an idea what the problem is?


texts: Offseting point perpendicular to line direction in PostGIS

texts:
I have a PostGIS table of points, that are sitting on road centerlines (I've linear referenced them). And the second table with road centerlines themselves.
The point table also has a field roadid and a field "offst", that represents the offset in meters (negative values = offest to the left, positive = offset to the right side of the road).
The question is how can I calculate a point geometry, that is offsetted perpendicular to road direction, according to the value of the "offst" field?


texts: Reclassify a Vector Landuse file using QGIS

texts:
I am trying to reclassify a landuse dataset using v.reclass in QGIS but my algoritm seems not to be working. It is a simple reclassification:
1 = 5, 2 = 0, 3 = 10, 4 = 3, 5 = 9, 6 = 4, 7 = 0
Does anyone has any experience using v.reclass or has another idea how to reclassify the data by maybe using another tool.


texts:
I generated a Lansat 5 image NDVI in a rectangle of interest in Google Earth Engine.

After that, I tried to export the image to my Google Drive, but the final result does not load the usual values of an NDVI (-1 to +1), as is usually done in a GIS program (QGIS, ArcGIS, etc.). The DN has a single value (attached image).

How could I solve this?

My code is:

var imgL5 = ee.Image('LANDSAT/LT05/C01/T2_TOA/LT05_226074_19850719');

var poly = ee.Geometry.Polygon([[
[-56.298454803020135,-20.398292995561647],
[-56.298454803020135,-19.8799857811524],
[-57.037285369426385,-19.8799857811524],
[-57.037285369426385,-20.398292995561647]
]])

var cortado = imgL5.clip(poly);

//Calc NDVI

var nir = cortado.select('B5');
var red = cortado.select('B4');
var NDVI = nir.subtract(red).divide(nir.add(red));

Map.addLayer(NDVI);

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: NDVI.float(),
  description: 'NDVI_9',
  scale: 30,
  region: geometry,
  folder: 'GEE',
});

Map.setCenter(-56.67318953949383,-20.121019484901392, 10);


texts: I have a simple web page where I show my CartoDB map and a dataset. I want to know if there is any way to insert new rows into my CartoDB dataset from code, so external users can add points to the map by clicking on it.


texts: Unknown spatial reference in "Create Fishnet" files

texts:
I am working on an unprojected dataset (lat/lon) on a WGS84 DATUM, and for that region I have used the "create fishnet" tool from "Data Management Tools" in ArcMap 10.2, in order to apply a grid over my region. Although, the coordinates I give as an origin of the fishnet are in lat/lon (in decimal degrees), the output of "create fishnet" is completely void of any spatial reference. When I import these shapefiles into a new project in ArcMap the message it appears is: "Unknown Spatial Reference: The following data sources you added are missing spatial reference information. This data can be drawn in ArcMap, but cannot be projected.". When I click "OK" the shapefile is imported and drawn correctly (correctly overlayed on a raster that is in lat/lon) but the outputs of "create fishnet" do not have any spatial reference.
How could I give a spatial reference to these shapefiles (both lines and label points) that are produced by "create fishnet" tool. I need the spatial reference to be in Geographic Coordinate System and DATUM WGS84


texts:
I want to access annotation layers using arcpy (the purpose is to verify few text strings in text annotations). By annotation I mean draw tool annotations, not dynamic labels. These annotation layers are under a groups layer by the name of 'Default'.
The Esri documentation talks of draw annotation as being layers themselves :

There are essentially three categories of layers in a map document: feature layers, group layers, and raster layers. The isFeatureLayer, isGroupLayer, and isRasterLayer properties allow you to identify or isolate the majority of layer types but not all layer types. There are a few specialized layers and datasets that don't fall into one of these three categories: annotation subclasses, dimension features, network datasets, terrain datasets, topology datasets, and so on. In these cases you may need to test other properties to isolate a layer of interest before doing something to it.

One thing to be aware of is how ArcGIS handles draw annotations, two mean configurations are possible :

ArcGIS fully supports two types of annotation: geodatabase and map document

In my case, it's map document.

In light of all this documentation, it so appears that draw annotations are mere layers, thus a simple code such as this should be working :
mxd = arcpy.mapping.MapDocument("some\mxdpath\here")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")
lyr = arcpy.mapping.ListLayers(mxd,"current_annotation_layer",df)

But the problem is the annotations do not seem obtainable via dataframe (which is weird considering each of my annotation layers is associated to a feature layer).

Is there any specific ArcPy class to access annotations (one comparable to IAnnotationFeature2 of ArcObject SDK)?


texts: QGIS Transfer Geomerty from Shapefile (loaded in QGIS) to PostgreSQL Table

texts:
I have seen a few postings relating to this question but the solutions are described with minimal detail, so I am looking for a MORE descriptive response. The solution I desire is via QGIS.
My procedure is as follows;

Connect to the Shapefile (source table)
Connect to the cloud PostgreSQL (PG) table (target table)
Enter a record (metadata only) in the PostgreSQL. The "geom" column is null for this record.
Save the edit in the PG table
Now open the attribute table for the PG Table and select the row
On the QGIS map enable the Shapefile table and query for the corresponding record and select the feature on the map which has geometry.
Enable the Advanced Digitizing Tools tool bar
Enable editing on the PG Table???
Select the "Add Part" button

Steps 8 and 9 is where I am unsure what else needs to be done, as well as the remainder of the process to save the geometry in the PG Table (target).


texts:
I´m trying to get some GRASS geoprocessig tools running in a python script outside GRASS.
I´m able to import grass.script, but I´m not able to get the tools running yet. Does anybody know what I´m doing wrong?

import os, sys
import subprocess as subp

gisbase = os.environ['GISBASE'] = "C:/GRASS GIS 7"
gisdbase = os.path.join(os.environ['HOME'])
location = "Martin"
mapset   = "PERMANENT"

sys.path.append(os.path.join(os.environ['GISBASE'], "etc", "python"))

import grass.script as grass
import grass.script.setup as gsetup

gsetup.init(gisbase, gisdbase, location, mapset)
grass.run_command("v.hull", input="D:\Python_Test\1.shp", output="D:\Python_Test\5.shp")

It runs up till grass.run_command, but as a result I only get -1073741515 in pythonwin. Once I get this running, I want to combine ArcPy and GRASS tools in one script. Could it have to do something with the fact that I´m using the ArcGIS Python installation (2.6)?

Sorry if this is a stupid question, I´m still a python beginner with not too much programming experience.


texts: What are the differences between Append and load data in ArcGIS?

texts:
I confused about differences between Append tool in Arctoolbox and Load data in Geodatabase.
I think both of them are the same in the result feature class.Both of them use a default output feature class and compare the target feature geometry and schema with another feature but there are differences in some aspects such as performance that I don't know .Another question is that where can i use Append instead of Load data ?

Thanks for advance


texts: Mapbox select road feature, how to return uncomposited features?

texts:
I'm trying to build an interface in which the user can click to select a road. You can see a simple version here: https://codepen.io/lizzie-malk/pen/eYBWbjm
What I have noticed is that the features selected in the Mapbox Streets style, are often composite and not individual roads, particularly when selecting a smaller road, see photo below:
I'm assuming this is Mapbox simplifying the OSM information into their tileset. Is there any way around this using a Mapbox tileset, or is there a way to use the uncomposited OSM features in my own style?


texts: Symbolize points in the same shapefile differently

texts:
I was wondering if there is a way to edit the appearance of the points in the same shapefile. I have a shapefile with fifty points in which dead wild animals were found but 5 of them were positive to some diseases. I would like to put 2/50 points with another symbol, 1/50 with a different and so on for the 5 points.

Is there a way to do this?


texts: How to join two LinestringM into one without dropping M value?

texts:
I have, for example, two Linestrings with M value:
LINESTRINGM(1 0 2, 2 0 4, 4 0 7)
and
LINESTRINGM(4 0 7, 6 0 9, 10 0 10)

How can I merge them into one Linestring and preserve the M value?

Expected result: LINESTRINGM(1 0 2, 2 0 4, 4 0 7, 6 0 9, 10 0 10)


texts: In MySQL are there unidirectional and directional graphs?

texts:
I'm new to MySQL and studying about geospatial data. I have read the documentation of MySQL and learned about the ST_Geometry type, its subtypes and its functions; but I didn't find anything about graphs (nodes, links or edges and paths)?

So my question is: Does this database have these datatypes?


texts: QGIS Identify tool Browse button not working

texts:
I am trying to add photo files to specific points in QGIS 2.18  using the instructions found in the training manual section 6.4
https://docs.qgis.org/2.2/en/docs/training_manual/create_vector_data/actions.html.

I have followed the direction exactly and created a field for images.  The problem is that when I try to add the file to the image field under filename  I can see the browse button but clicking it does not open a file browsing dialogue.  It does nothing.

Is there another way to specify my file path or another way to add non-geotagged photos to specific points?


texts:
Bit of a challenge here which I've been grappling with for some time. I'll explain my full work flow so you can reproduce if needed.

I'm creating virtual landscapes in Google SketchUp which I ultimately would like to use in Netlogo to examine how turtles interact with them.

My problem is that by the time I get the landscapes into Netlogo the units don't seem to relate to the original 3D model.

Step 1: Create simple hill on a 50m by 50m square in Sketchup using the Toposhaper extension.

Step 2: Export to .dae file and import into Meshlab, ensure the Meshlab model has the same dimensions as the Sketchup model by adjusting the units with the assistance of the measuring tool. Export from meshlab as .xyz file.

Step 3: Import .xyz file into QGis as points by adding a new layer from delimited file. Selecting field_1 and field_2 as X and Y fields.

Step 4: Create raster of points using Raster > Interpolation > Interpolation. Add field_3 as interpolation attribute, set number of columns to 50 by 50 (to correspond to the 50m x 50m 3D model). The cell size of X and Y should match exactly to ensure Netlogo will read the resulting .asc file.

Step 5: Finally, I setup a model in Netlogo to receive the raster. Firstly, in model settings I set the the min and max pxor and pycor to 0 and 49 (creating 50 x 50 patches). Then, using the Gis Extension, I import the raster apply the z-value to a patch variable called elevation:

to load-gis
  set elevation gis:load-dataset "cone_50.asc"
  gis:set-world-envelope-ds gis:envelope-of elevation
  gis:apply-raster  elevation target-elev
end

Now, each patch of my 50 by 50 Netlogo world should have an elevation value taken from my 50 by 50 raster. In theory, adding all the elevation values together should (roughly) give me the total volume of the raised area of the hill? The figure I get is higher however and the problem gets worse with larger volumes.

Can anyone help?


texts: Seeking finer points of slope statistics within set of polygons using ArcMap

texts:
I have a few questions related to finding slope statistics (min, max, mean, standard deviation) of slopes  (1-meter LiDAR DEM derived) within a set of polygons (census blocks).
First, is it valid to average slopes (degree or percent)?
This is a big question I guess but could someone tell me if there is a way to use slope percent in the averages.
Are there any other methods other than simply adding up the values of raster cells and dividing by the number of cells?   (not addressed elsewhere)
Would it make more sense to present every slope value available (within each polygon) at the highest resolution (1-meter) and present the range within 1 or 2 standard deviations of the mean slope?
One of the key statistics that I am looking for is the range of slopes within each polygon but perhaps presenting slope values within one standard deviation would suffice.  Currently my method is to resample a 1-meter DEM to 10-meters to remove the extreme values and then use zonal stats to get the mean, min and max (this question is not addressed elsewhere)  ...
Is there any other way to average the slope that takes into account partial cells?
One method I thought of was to resample the raster down to a finer resolution.. e.g. 30-meters to 1-meter and then averaging the slope using zonal stats.   This is kind of a rough integration.
Is there a more elegant/perfect way to do this? (not addressed elsewhere)
I'm using ArcGIS 10.2 for Desktop.


texts: Sort selected points by latitude/northing with field calculator

texts:
I have bathimetry data in .shp format. The points (picture below shows one of multiple groups/profiles) do not have ID numbers and are randomly dispersed through the attribute table.

I have selected each group and added group ID to it. Now each point has "Sort" attribute that serves as the group ID. I still need each individual point to have ID number within the group.

I need these points numbered starting with the northernmost point in group (Field2 in attribute table is northing) and ending with the southernmost. So, a point with highest "Field2" value gets value of 1 in "point_num" field.

Basically, what I need to  do is sort selected points according to value of "Field2" attribute. How can I do this with (preferably) QGIS field calculator or some other tool?


texts: Legend symbol colors display incorrectly using data-defined override for colors with categorized style

texts: I am able to style the colors of a geology polygon shapefile according to the Red, Green and Blue fields of the attribute table using data defined override with the expression: "color_rgb("RED" , "GREEN" , "BLUE")", which works great and the colors on the map canvas render properly. However the colors in the legend symbols are totally random and don't correspond to the colors on the map. I see a similar issue used to arise where the legend symbols were all black in older versions of QGIS...


texts: I've got a pop-up for a layer with more than one attribute for the layer in question. I would like to make the first attribute in the popup appear with a bold font and a different color. Anyone have any idea?


texts: QGIS table does not dynamically update after applying a filter

texts:
Using QGIS 2.10.1, if I open a table and then apply a filter using Layer Properties > General > Provider Feature Filter, the table view does not update.  In other words, the initial list of records remains displayed in the table window.  However, if I close and then reopen the table window the filtered records are correctly displayed.

Is this table behavior correct?  I would expect that the table would automatically update after a filter is applied.


texts:
This post couldnt help my problem. I try to find out how many countries start with a 'd' and have a 'n' in their name using the query.
When I type in this, WHERE CNTRY_NAME LIKE 'd%n%' I get the following Error
SQL Expression Parsing Error: syntax error, unexpected WHERE. Occurred around :
WHERE CNTRY_NAME LIKE 'd%n%'
^

Can someone tell me where the problem lies?


texts: Qgis-How to make multiple print composers with one drawing using python

texts:
I am trying to build multiple print composers but most of the sheets that I have to work with use the same layers with an exception to one unique layer. Is it possible to create one map document with all the layers I need and then distribute the layers to various print composers that I specify in code?

I think I am creating way too many map documents when I can just create one and a python script would automatically set up print composers that corresponds to a standard I specify. I am new to the qgis library for python. Is there any base code I can work off of or any useful links that I should see?


texts: How to convert QGIS 3.0 project to use auth_db ids

texts:
I am working on a Mac and running QGIS 3.0.3 and I am exploring the authentication db setup.

I have a bunch of layers almost all of which are based on PostGIS queries and I want to move all my authentication into the managed auth_db and it is not clear to me how to proceed.

I assume I will have to edit the project file and replace the database connection parameters with a reference to the id but how should I do this?

One thing I want to achieve with this exercise is to make the project files so they are independent of the user so I can share project files as a starting point for people .


texts: What is the expected SRID value in Teradata ST_Geometry Methods?

texts:
In Teradata database there's SPATIAL_REF_SYS table in SYSSPATIAL schema. In my Teradata production database the table lists 3,167 different Spatial Reference Systems (SRSs).

Now, it seems that Teradata Geospatial has its own numbering for specifying SRSs. Let's check the below query.

SELECT *
FROM SYSSPATIAL.SPATIAL_REF_SYS
WHERE AUTH_SRID=4326;

Result:

SRID  AUTH_NAME  AUTH_SRID  SRTEXT
----  ---------  ---------  -----------
1619  EPSG            4326  (see below)

the content of SRTEXT above is:

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.01745329251994328]]

Now the question is:

Which one is the valid SRID expected in all ST_Geometry Methods that take SRID as their parameter (eg: ST_SRID, ST_Geometry constructors, ST_Transform): the one as in SRID or AUTH_SRID column?
Similarly, which one is the valid SRID expected in the TDGeoImport Command (the tool for loading geospatial data into Teradata database)? Also how to type it as SRS parameter (with "epsg:" or simply just the number? The example in Setting a Spatial Index in the PDF "GeoServer_Teradata_v1.0" here does not include SRID parameter during TDGeoImport.


texts: Consistent label sizing with Data Driven Pages

texts:
I've set up a map document with a number of maps, using DDP with a polygon feature class as an index layer. Two of the layers have labels and I need to manually adjust the position of many of the labels. The only way I know to do this is by converting the labels to annotation, however because the index layer polygons and map data frames are all different sizes this means that the maps are drawn at different scales and so the labels are different sizes on each map. I would like them to all to be font size 10 on all maps.

Is there a way to do this?


texts: ArcGIS Enterprise 10.8.1: Trying to create feature layer with referenced data in Python

texts:
I am trying to write a script that creates feature layers in ArcGIS Enterprise based on layers in an ArcGIS Pro Project. I do not want to copy the data to the server, only reference it. The data comes out of a SDE database that is also registered with the respective ArcGIS Server.
While doing the process manually works well (right-click layer --> Share as Weblayer --> Reference registered data), if I try to do the same with Python, I always end up with a hosted feature layer, which means that the data has been copied to the server instead of referencing it.
This is the code I use. Am I missing something?
            featDraft = Map.getWebLayerSharingDraft("HOSTING_SERVER", "FEATURE", "REF_" + service6, layer)
            featDraft.copyDataToServer = "False"
            featDraft.description = metadata_description
            featDraft.allowExporting="True"
            featDraft.overwriteExistingService="True"
            featDraft.portalFolder="REF_" + Mapname
            featDraft.tags=Tags
            mapDoc.save()
            featDraft.exportToSDDraft(sddraft_output_filename)

            # Stage Service
            sd_filename = service6 + ".sd"
            sd_output_filename = os.path.join(outdir, sd_filename)
            arcpy.StageService_server(sddraft_output_filename, sd_output_filename, 5)

            # Share to portal
            print("Uploading Service Definition...")
            arcpy.UploadServiceDefinition_server(sd_output_filename, "My Hosted Services")

Follow up:
Putting a arcpy.GetMessages(1) into my code gives me the following error:
WARNING 000710: Unable to connect to the server

I assume that this means that the script is not able to access the database where the data resides to reference it. However, I am not sure where to put this information. Any ideas?


texts:
I have styled my point data with circle and red color symbology. Is there any way to retain the points color when saving layer to KML?

Currently when saved and open in Google Earth, it only shows as a pin.


texts: Mapping Longitude and Latitude to US state

texts:
Is there any database (Gov data or private) or web service that map Longitude and Latitude to US state?

I have the specific Longitude and Latitude of a location and would want to identify the state that this location is in. Then repeat this for thousands and millions of records.


texts: run r.random through script in python shell

texts:
Trying to run r.random using the script but receiving an error message. It works when the module is run outside of the python shell so it must be te code, I'm new to this so I can't figure out what is wrong.

input:

grass.run_command("r.random", overwrite=True, input="island",npoints=15, raster="randomsample", vector="randomsample")

error:

Process ended with non-zero return code 1. See errors in the (error) output.


texts: Data source not found when opening .qgs file

texts:
I am using QGIS 3.6.1 with the new geopackage (.gpkg) format. When some of my co-workers open my .qgs-file, QGIS (same version as mine) does to not find the data source for most of the layers. Strangely, most of the layers appear 3 times. I get the same when I try to add a layer to my geopackage database. In my database, I can see only 1 layer with the same name.

When I open the .qgs-file, I don't have any problems. My co-workers have read-only access to the folder where the .qgs-file is located. Maybe it has something to do with it?

How can I fix the data source problem?


texts: Is there any available way to access raster images through GeoServer without WMS? (to a large number of images)


texts: Connectivity rule regarding two-level intersections

texts:
in the sDNA documentation here: https://www.cardiff.ac.uk/sdna/wp-content/downloads/documentation/manual/sDNA_manual_v4_0_1/network_preparation.html
is written that the network needs to be in shape as shown in Figure 5. Also at the bottom is written:

  Summary: ensure you understand spatial network connectivity rules, and
  (for sDNA) that your network uses coincident endpoint connectivity.

Now comes my problem: I have a crossing exactly as shown in Figure 4, which has two different levels. One street A runs at elevation 0m from ground and runs under street B, which is at elevation 5m from ground. Street A is a usual road on the ground, while street B is a bridge running over street A, without having any possibility to get from street A to street B.
When I want to prepare my network for the sDNA tool, I came across this problem and I am not sure how to solve it, since the documentation states I would need an intersection, while in reality there is no possibility to get from A to B or vice versa. I am a bit lost on how I should tackle this problem, since there is no way to change from A to B, while the documentation suggests I should add an intersection which is not available in real life.

(edited to correct manual url - manual 4.0.0 is deprecated)


texts: I uploaded some shapefiles to GeoServer and I added them to OpenLayers map using WFS protocol. Now I want to style them by some properties. How can I do this?


texts: Plot a Raster of Weather data in Leaflet JavaScript Map without generating server-side tiles?

texts: I want to show different weather data I have as binary arrays on a server using leaflet. But I do not want to generate tiles on the server. For each point in an a array I have the coordinates in another array. I also have the proj-string for the projection. Is it possible to retrieve the arrays from the server and dynamically generate the map for leaflet in the client? If so, how?


texts: arcpy.da.Editor undoOperation(): does it actually work

texts:
Ok, so I have looked over the ESRI documentation thoroughly and I'm still baffled by the arcpy data access module and the function undo operation.

So far I have tried:

workspace = pathtofile
edit = arcpy.da.Editor(workspace)
edit.startEditing(with_undo = True)
edit.startOperation()
edit.isEditing
arcpy.CalculateField_management("MyTable", "LIFEFORM", 5000)
edit.isEditing
edit.undoOperation()

The undo operation does nothing and the undo button in the user interface indicates that it will remove the added layer, not the field calculate that I just performed.

I was hoping that I could implement the da.Editor into an addin toolbar as a way of avoiding the user from having to turn on an extension or start an edit session on a separate toolbar.


texts: Converting lat/long to EPSG:3857 coordinates

texts:
The map widget I'm using gives me click locations in lat/long.  I need to convert those coordinates to EPSG:3857 to query a database.

Most of the solutions I've found use the formula described in this gist: convert from long/lat to google mercator (or EPSG:4326 to EPSG:900913) .  But I'm in doubt about the constant 90 used in the y calculation:

var y = Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);

The formula described on Wikipedia seems to suggest that instead of 90, I should use 85.051129.  I've found other references to that value in questions such as Transforming EPSG:3857 to EPSG:4326.

Which is correct, and why?


texts: Creating Dialog with Graphic Modeler for Attribute Creation

texts:
I have a shapefile with an attribute that I need to convert based on categorical rules. Hereby the TARGET / new values always stay the same and only the input values change. For example:

"Residential Area" assign if value "TYP" = 1
"Industrial Area" assign if value "TYP" = 2
"Mixed Area" assign if value "TYP" = 3 OR 4

I would like to create a dialog where the user is prompted to assign each original value to a set of predefined possible options.

I know I could do this in the field calculator but I would like it to be in the form of a menu mask where you set:

Input Shapefile
Value Assignment based on a table where all target variables are set and te input has to be select from a pulldown menu.

I have looked in the value map widget but have not found what I am looking for.

Any ideas?


texts: Can I import another model into my model in Modelbuilder?

texts:
I am new here and new to Modelbuilder.
I have been trying to create a toolbox. I have managed to branch down the tool into simple models in Modelbuilder but in most cases I can not combine my models together. I mean, they models work properly in separate models but I can not find a way to have them all in one model. Sometimes it is impossible to join them together and sometimes it doesn't work properly. I was wondering what is the best way to combine them together. should I convert all the models to script and then combine them together and use the tool as a script or there is a  way to export/import the models.
Thank you all in advance for helping me on this.

Arash.


texts: Is it possible to create Thiessen polygons around nodes in QGIS? e.g. http://en.wikipedia.org/wiki/Voronoi_diagram


texts: Using Intersect function gives error code Intersect [175528_07022018]?

texts:
I´m doing field selection studies for red deer in ArcGIS. I have several years of data and have made buffer polygons base on the red deer gps positions. The buffer polygons  should be integrated with the field polygons with cultivated crops for each year (to be able to create random positions to represent field availability).

My problem is that the intersect function is not working for some few years and I repeatedly get error messages. Several of the years have been doing fine with the intersect function, but is suddenly not working for some of the years. I have tried so many things to solve it, like subsetting the dataset to make it less heavy (worked for some years), and redoing the buffer polygons, but I cannot get it to work.

I get the error code Intersect [175528_07022018] and I have tried to google it without any hits.

Do anyone have a solution to this problem or know how to interpret the error code?


texts: Restore a PostGIS database cluster properly

texts:
Context

In order to migrate from PostgreSQL 10.12 to 12.2 on a Ubuntu 18.04.4 LTS server, I first uninstall and purge all postgres related stuff then I followed this;

https://www.postgresql.org/download/linux/ubuntu/

As every major operation on things around system or databases, I took the habit to make some backups.

To do this, I used the content that I've found there;

https://www.postgresql.org/docs/current/upgrading.html

Especially this command in order to dump all my databases at once:

$ pg_dumpall > outputfile.dump

Errors

After migrating to the new postgres version, I first faced the error that the user postgres already exists, which was fixed by commenting the two following lines in the dump file (outputfile.dump):

$ sed -ne 's/CREATE ROLE postgres/--&/g' outputfile.dump
$ sed -ne 's/ALTER ROLE postgres/--&/g' outputfile.dump

Postgis related error

Of course, then I was confronted with another error, otherwise it would be too easy.
This one is as follows;

$ psql -U postgres -f outputfile.dump -v ON_ERROR_STOP=ON
SET
SET
SET
ALTER DATABASE
CREATE DATABASE
ALTER DATABASE
REVOKE
GRANT
You are now connected to database "postgres" as user "postgres".
SET
SET
SET
SET
SET
SET
 set_config
------------

(1 row)

SET
SET
SET
SET
psql:outputfile.dump:70: ERROR:  schema "topology" already exists

So, ... I tried to think about the same thing that before, i.e. commenting the line 70 which is;

70:CREATE SCHEMA topology;

but... this line actually creates the topology schema for PostGIS, which I obviously need as my databases are spatial ones. And indeed, if I comment out those lines, e.g. by:

$ sed -ne 's/CREATE SCHEMA topology/--&/g' outputfile.dump

I face a schema "topology" does not exist error.

The line CREATE SCHEMA topology which creates the topology, is followed by these:

73:ALTER SCHEMA topology OWNER TO postgres;
76:-- Name: SCHEMA topology; Type: COMMENT; Schema: -; Owner: postgres
79:COMMENT ON SCHEMA topology IS 'PostGIS Topology schema';

But what is really interesting here, is that the postgis_topology extension itself is actually created after the previously described lines, e.g. in my case at line 153 for the very first database in the cluster:

153:-- Name: postgis_topology; Type: EXTENSION; Schema: -; Owner:
156:CREATE EXTENSION IF NOT EXISTS postgis_topology WITH SCHEMA topology;
163:COMMENT ON EXTENSION postgis_topology IS 'PostGIS topology spatial types and functions';
175:-- Data for Name: topology; Type: TABLE DATA; Schema: topology; Owner: postgres
178:COPY topology.topology (id, name, srid, "precision", hasz) FROM stdin;
183:-- Data for Name: layer; Type: TABLE DATA; Schema: topology; Owner: postgres
186:COPY topology.layer (topology_id, layer_id, schema_name, table_name, feature_column, feature_type, level, child_id) FROM stdin;
217:-- Name: topology; Type: SCHEMA; Schema: -; Owner: postgres
220:CREATE SCHEMA topology;
223:ALTER SCHEMA topology OWNER TO postgres;
226:-- Name: SCHEMA topology; Type: COMMENT; Schema: -; Owner: postgres
229:COMMENT ON SCHEMA topology IS 'PostGIS Topology schema';

These lines obviously repeats each time a PostGIS database has to be restored.

Question

I wonder why the error about the topology schema is raised before (@line 70) the extension is created (@line 156)?
This seems strange to me because I thought the backup file was supposed to be read sequentially.
Could it be possible the dump tool itself (pg_dumapp) doesn't work in a proper way, in which case I'd better ask at dba.SE or elsewhere?


texts: Connecting to ArcSDE with arcpy is slow (ArcSDESQLExecute)

texts:
I'm creating an API which will allow queries to be performed on my GeoDatabase. I'm using ArcSDESQLExecute to connect to ArcSDE (backed with an Oracle database). Creating the connection object is very slow (~1.5 s):

sde_conn = arcpy.ArcSDESQLExecute(r"data\Connection to GPSERVER3.sde")

So doing this before while handling an API request makes response times very slow.

However, if I create sde_conn once, for the life of the API, and re-use it indefinitely, queries are much faster (~20ms in some cases), since it's not recreating sde_conn each time.

Is it bad practice to keep sde_conn alive? I suspect that creating the object creates a persistent connection. Does anyone know what it's doing under the hood? (is it opening and closing the connection during each call to execute ?


texts:
I'm trying to create image pyramid  of multiple smaller images for GeoServer use. I have done this by first creating virtual dataset with gdalbuildvrt and then with gdal_retile.py.

The problem I'm having is, that the original images are disjoint, and thus the virtual dataset has a lot of empty space with no data (rendered as black). I would want to render this area transparent.

The images are all JPEgs with YCbCr colorspace thus they do not contain alpha channel, nor does the original data have nodata value set. Also the images contain all pixel values in 0-255 range (i.e. I cannot use -a_nodata of gdal_translate without getting holes in the picture). I have narrowed it down to two solutions: One would be to use external mask, but I'm wondering the performance impact of this solution. Not to mention the need to store an extra image that is almost completely empty. The second solution would be to edit every single pixel in the original images so that none would have value 0. then use 0 as nodata value.

Any ideas which method I should use? Or neither?


texts:
I'm newbie in GIS.
I want use the shapefiles of USA map for my Geoserver. Where can I download shapes for USA? I know about the TIGER\Line shapes but it is separate shapes of counties and thematic maps (e.g. american indian areas). I want something like OpenStreetMap.
P.S. Sorry for my english )


texts:
I have a table with 6 columns (Integer) containing licence numbers.

I would like to label the data using only the columns containing a licence number and not have a label for the columns containing a 0 (zero) attribute.

The column PRIMARY_LICENCE_NUMBER contains only licence numbers but the other columns have a mixture of licence numbers and a zero, reducing until License_6 has no licence numbers and only zeros.
The largest licence number is 118001 and the smallest is 11001.

I want to label in this order:
PRIMARY_LICENCE_NUMBER,License_2,License_3,License_4,License_5,License_6

I have read your other answers with expressions containing Left$ and Int and have not had any luck working out an expression that works.

My bumbling attempt is at this stage:

"L1: "+PRIMARY_LICENCE_NUMBER+Chr$(13)+"L2: "+Left$(License_2, Len(License_2)*Int(License_2<>0))

and all variations in between, and I can't figure out how to fix this issue.


texts: Aerial picture - how can I extract water drainages on a crop (to vectorize them)

texts: I have a RGB aerial picture of an agricultural crop. The drainages (hundreds of them) can be barely seen because they a very slightly greener then the surrounding crop area (picture is from April). How is it possible to make these drainages more visible in order to be able to vectorize them?


texts: Unique value symbology not loading although values exist in ArcGIS raster attribute table?

texts:
I have a weird problem. I have a feature class generated from processing multiple rasters and converted by Raster to polygon tool in ArcGIS 10.4.

I would like to apply Unique value symbology based on field gridcode in ArcMap.
The individual gridcode values (from 0 to 14) exist in the field gridcode in Attribute table, I can also Label them, as here:

However, when I want to apply Unique value symbology, only the value "1" gets load after "Add all values".

Do you know why the other values do not get loaded?


texts: gdal created tiles for the U.S. are displayed south of equator

texts:
I am using gdal_translate to convert a map in spheroid meters to EPSG:900913.  I then use gdal2tiles.py to create tiles.
I am centering my googlemaps using Lat and Long in the central U.S.. I am overlaying the tiles onto Googlemaps.   The tiles are displayed south of the Equator?  How do I get the tiles to center over the U.S.


texts: Why does the QgsFeatureIterator crash in my stand alone script?

texts:
I'm trying to write a standalone script in qgis 2.10 in ubuntu
I have written a script which crashes when I run from bash but works fine from inside the qgis desktop environment. The script is:

#!/usr/bin/python2.7
from qgis.core import *
from platform import python_version
print python_version()
# supply path to where is your qgis installed
QgsApplication.setPrefixPath("/usr", True)
# load providers
QgsApplication.initQgis()

input_file = "/path/to/shape/file.shp"
input_layer = QgsVectorLayer(input_file, "input layer", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(input_layer)
print str(input_layer.featureCount())
print input_layer.dataProvider().crs().authid()
for feature in input_layer.getFeatures():
    pass
print "hello"

QgsApplication.exitQgis()

When I run from within the qgis DE I get:

2.7.6
105777
EPSG:2193
hello

When I run from bash I get:

2.7.6
105777
EPSG:2193
Segmentation fault (core dumped)

The only thing I can guess is my interpreter choice in my shebang but none of the python executables I can find make any difference? and given the response to the version function that seems unlikely. For the shebang I also tried

#!/usr/bin/env python

with the same result.

I've seen the docs on standalone scripts here
and I figured out the correct(?) path to the installation from here but I'm stumped as to the cause of this crash.


texts: My osm2po java program compiles but gives a runtime error

texts:
I need road distances and I'm trying to use osm2po to do this.

After downloading the OSM files and converting them to gph format, using the osm2po program, I start at the very beginning with the following very simple java program:

import java.io.File;
import java.util.Arrays;

import de.cm.osm2po.logging.Log;
import de.cm.osm2po.logging.Log2poConsoleWriter;
import de.cm.osm2po.model.LatLon;
import de.cm.osm2po.routing.Graph;
import de.cm.osm2po.routing.MultiTargetRouter;
import de.cm.osm2po.routing.PoiRouter;

public class TestOSM {

    public static void main(String[] args) {

        File graphFile = new File("C:/osm2po/osm2po-5.0.0/hh/hh_2po.gph");
        Graph graph = new Graph(graphFile);

        // Somewhere in Hamburg
        int sourceId = graph.findClosestVertexId(53.5f, 10.0f);
        int targetId = graph.findClosestVertexId(53.4f, 10.1f);
        graph.close();

    }
}

I get it to compile using:

javac -d . -classpath osm2po-core-5.0.0-signed.jar TestOSM.java

It compiles and creates the TestOSM.class file.

I then use the following to run the program:

java TestOSM

It gives a runtime error (I'll show the error message below)

I suspected the program does not read the hh_2po.gph file, so I used another java program to first use "normal" java file read to confirm successful file reading operation.

My java program to try to confirm successful file reading looks like this:

import java.io.File;
import java.util.Arrays;

import java.io.Console;
import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.io.FileInputStream;
import java.io.InputStream;

import de.cm.osm2po.logging.Log;
import de.cm.osm2po.logging.Log2poConsoleWriter;
import de.cm.osm2po.model.LatLon;
import de.cm.osm2po.routing.Graph;
import de.cm.osm2po.routing.MultiTargetRouter;
import de.cm.osm2po.routing.PoiRouter;

public class TestOSM {

    public static void main(String[] args) {

        // Testing take1
        File graphFile = new File("C:/osm2po/osm2po-5.0.0/hh/hh_2po.gph");
        //File graphFile = new File(args[0]);  // tried this, and then supply file name when calling the program - does not work
        //File graphFile = new File("hh_2po.gph");  // tried this after copying hh_2po.gph in the same directory - does not work
        Graph graph = new Graph(graphFile);

        // Somewhere in Hamburg
        int sourceId = graph.findClosestVertexId(53.5f, 10.0f);
        int targetId = graph.findClosestVertexId(53.4f, 10.1f);
        graph.close();

/*
        // Testing take2
        // Part 1, to confirm that the program can read the hh_2po.gph file:
        try{
            BufferedReader reader = new BufferedReader(new FileReader(new File("C:/osm2po/osm2po-5.0.0/hh/hh_2po.gph")));
            String  thisLine = thisLine = reader.readLine();
            System.out.println(thisLine);
            System.out.println("file read good");
            reader.close();
        }catch(IOException ioe)
        {
            System.out.println("could not read file");
        }
        System.out.println("enter to continue");
        Console console = System.console();
        String line = console.readLine();
        // end of test

        // Part 2, this time try to convert file to InputStream before converting to a Graph:
        try{
            File initialFile = new File("C:/osm2po/osm2po-5.0.0/hh/hh_2po.gph");
            InputStream targetStream = new FileInputStream(initialFile);
            Graph graph = new Graph(targetStream);
        }catch(IOException ioe)
        {
            System.out.println("Heuston we have a problem");
        }
*/

    }
}

When I run it, again using java TestOSM, I get the following runtime error:

I use Windows10 and java jdk1.8.0_60


texts: If statement in Python QGIS Raster Calculator

texts:
I would like to perform a raster calculation with the if condition and I wrote this code:
import qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

inputpath='/Users/macbook/Desktop/TESI/PROGETTO/rasterUfinali/U_1990.tif'
inputpath2='/Users/macbook/Desktop/TESI/PROGETTO/SIGMAmax/SIGMAmax,har.tif'
inputpath3 = '/Users/macbook/Desktop/TESI/PROGETTO/sigm/sigmahar90.tif'
inputpath4 = '/Users/macbook/Desktop/TESI/PROGETTO/lambdamin/lambminhar90.tif'
outputpath = '/Users/macbook/Desktop/TESI/PROGETTO/alfamax/alphamaxhar90.tif'

inputrasterfile = QgsRasterLayer(inputpath)

entries = []
ras = QgsRasterCalculatorEntry()
ras.ref = 'ras@1'
ras.raster = inputrasterfile
ras.bandNumber = 1
entries.append( ras )

inputrasterfile2 = QgsRasterLayer(inputpath2)

ras2 = QgsRasterCalculatorEntry()
ras2.ref = 'ras@2'
ras2.raster = inputrasterfile2
ras2.bandNumber = 1
entries.append( ras2 )

inputrasterfile3 = QgsRasterLayer(inputpath3)

ras3 = QgsRasterCalculatorEntry()
ras3.ref = 'ras@3'
ras3.raster = inputrasterfile3
ras3.bandNumber = 1
entries.append( ras3 )

inputrasterfile4 = QgsRasterLayer(inputpath4)

ras4 = QgsRasterCalculatorEntry()
ras4.ref = 'ras@4'
ras4.raster = inputrasterfile4
ras4.bandNumber = 1
entries.append( ras4 )

if '(ras@2 - ras@1) > 0':
    alpha_max='1'
else: alpha_max = 'ras@3 / (1-2.71828^(ras@4))'

new_path = outputpath
print(alpha_max)

alfamax= QgsRasterCalculator(alpha_max, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )
alfamax.processCalculation()
iface.addRasterLayer(new_path)

The code works and print a valid output raster, but it doesn't respect the "if condition", because it prints out alpha_max = 1 even where (ras@1 - ras@1)<0.
Can someone tell me where is the mistake?


texts:
I'm having three PostGIS layers (points, lines, polygons) that I'd like to export to a single DWG file.
So far, I could export the layers as an AutoCAD DXF file with a right-click on the layer and choosing Export--> Save Selected Features As ... .
I could read about the export to dwg with a special gdal compilation, but can't use this as we only use a QGIS3 LTR version under windows10.
Does somebody have a workaround for creating the DWG file?


texts: Extracting raster value to specific range using ArcPy?

texts:
I have multiple rasters with range 0 - 255. I want to extract range from 0 to 100 using Arc GIS python script and after extraction I want to multiply all raster by 0.01 and save it with same name as it input using loop function. I have made one script till raster extraction, but I am getting an error while running the script. Below I have shown my script.

import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace = r"D:\MODIS-NDVI\FPAR-2005\MASKED-FPAR-05"
outFolder = r"D:\MODIS-NDVI\FPAR-2005\MASKED-FPAR-05A"
arcpy.CheckOutExtension("Spatial")
list = arcpy.ListRasters
for inRaster in rasters:
    outRaster = outFolder + "\\" + inRaster
    arcpy.sa.con(inRaster >3000,inRaster)


texts: How to add a topographic map to a shapefile?

texts:
I wonder how it's possible to add this topographic information (elevation) to the map of Santa Cruz Island and all of the Galápagos Island that I have here.

library(rgdal)
library(rgeos)
library(ggplot2)
library(ggthemes)
library(ggsn)
myGPS = data.frame(lat=c( -0.6850556,-0.6854722,  -0.6857778  ),lon=c(-90.22275,-90.22261,  -90.22272))
coord.deg = myGPS

class(coord.deg)
## "data.frame"
coordinates(coord.deg)<-~lon+lat
class(coord.deg)
## "SpatialPointsDataFrame"
## attr(,"package")
## "sp"

# does it have a projection/coordinate system assigned?
proj4string(coord.deg) # nope
## NA

# Manually tell R what the coordinate system is
proj4string(coord.deg)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

# now we can use the spTransform function to project. We will project
# the mapdata and for coordinate reference system (CRS) we will
# assign the projection from counties

coord.deg<-spTransform(coord.deg, CRS(proj4string(gal)))
# double check that they match
identical(proj4string(coord.deg),proj4string(gal))
## [1] TRUE
my_pts <- SpatialPointsDataFrame(coords = coord.deg, data=data.frame(id=1:length(coord.deg)))

URL <- "https://osm2.cartodb.com/api/v2/sql?filename=public.galapagos_islands&q=select+*+from+public.galapagos_islands&format=geojson&bounds=&api_key="
fil <- "gal.json"
if (!file.exists(fil)) download.file(URL, fil)
gal <- readOGR(fil, "OGRGeoJSON")
gal <- gSimplify(gUnaryUnion(spTransform(gal, CRS("+init=epsg:31983")), id=NULL), tol=0.001)
gal_map <- fortify(gal)
rand_pts <- spTransform(my_pts, CRS("+init=epsg:31983"))

# ggplot can't deal with a SpatialPointsDataFrame so we can convert back to a data.frame
my_pts <- data.frame(my_pts)
my_pts.final = my_pts[,2:3]
# we're not dealing with lat/long but with x/y
# this is not necessary but for clarity change variable names
names(my_pts.final)[names(my_pts.final)=="lat"]<-"y"
names(my_pts.final)[names(my_pts.final)=="lon"]<-"x"

gg <- ggplot()
gg <- gg + geom_map(map=gal_map, data=gal_map,
                    aes(x=long, y=lat, map_id=id),
                    color="black", fill="#FFFFFF", size=.5)
gg <- gg + coord_equal()
gg <- gg + ggsn:::scalebar(gal_map, dist=50, location="bottomleft", st.size=5)
gg <- gg + theme_map()
gg <- gg + geom_point(data=my_pts.final, aes(x=x, y=y), color="red")
gg

Do you have any idea? Do you have a package in R that does this?

I added this:

topo <- readOGR("public.geojson.json", layer ="OGRGeoJSON")
topo <- spTransform(topo, CRS("+init=epsg:31983"))
gg + geom_map(map = topo_map, data = topo_map, aes(map_id=id))

But it creates a map with a black shape in it...


texts: Qgis is slow with complex 'postgis query' layer. How to speed up?

texts:
I have a postgis layer in Qgis that is a result of a complex query. Every screen update, Qgis runs the Query again and that's slow.
Any suggestion on how to speed up?

EDIT - Here is the query:

 SELECT
 parcela.id AS id,
 st_union(areas.the_geom) AS geom,
  (((((('p'::text || round(avg(to_number(dado.valor -> 'pequena'::text, '99'::text)) * 10::numeric, 2)) || ',m'::text) || round(avg(to_number(dado.valor -> 'media'::text, '99'::text)) * 10::numeric, 2)) || ',g'::text) || round(avg(to_number(dado.valor -> 'grande'::text, '99'::text)) * 10::numeric, 2)) || ',mo'::text) || round(avg(to_number(dado.valor -> 'morta'::text, '99'::text)) * 10::numeric, 2) AS "legenda",
  round(avg(to_number(dado.valor -> 'pequena'::text, '99'::text)) * 10::numeric, 2) AS pequenas,
  round(avg(to_number(dado.valor -> 'media'::text, '99'::text)) * 10::numeric, 2) AS medias,
  round(avg(to_number(dado.valor -> 'grande'::text, '99'::text)) * 10::numeric, 2) AS grandes,
  round(avg(to_number(dado.valor -> 'morta'::text, '99'::text)) / 10::numeric, 2) AS mortas

  FROM petr.dado, petr.parcela, petr.parcela_rel, petr.areas
  WHERE dado.parcela_id = parcela.id AND parcela_rel.area_id = areas.gid AND parcela_rel.parcela_id = parcela.id
  GROUP BY parcela.id