texts: Using SAGA-GIS topographic wetness index in R?

texts:
I want to create a topographic wetness index (TWI) in R and I am trying to set an environment according to the help of RSAGA-package but I am always facing with an error appears in a different term.

my code is as follows:

library("RSAGA")
library("envirem")
myenv <- rsaga.env()
myenv <- rsaga.env(workspace = ".",path ="C:/OSGeo4W64/apps/saga-ltr"
               ,modules ="C:/OSGeo4W64/apps/saga-ltr/modules" )
rsaga.hillshade("dem1","hillshade205",env=myenv)
top<-topoWetnessIndex(dem1, myenv)

the error

Error in .getGDALtransient(x, filename = filename, options = options,  :
filename exists; use overwrite=TRUE


texts: How do I remove a portion of a vector layer line string?

texts:
This is probably basic but I am a complete newbie. I searched for editing layers but am having difficulty understanding it. Here's what I have and understand so far.

I imported a KML file as a vector line. I understand I can't edit the KML import directly. I need to first save it as an ESRI shapefile. I have done that. I also toggled edit mode for the new layer and turned on the Vertex tool, which I think (not sure) is the same as the Node tool I've seen descried elsewhere - at least the icon looks the same. But then I am completely lost.

In the graphic, the magenta line is the complete track. The yellow box represents the portion I am trying to keep. I would like to delete the rest. I'm hoping this does not require a point-at-a-time deletion. For example, the segment in the yellow box from the loop on the top left to the "racetrack" at the middle right is about 9 miles. That would be a lot of points!


texts: Generating a heatmap with the twitter search term using cartodb api

texts: I have a search box in my UI and a div for displaying maps.if i enter any term in the search box, say for example (chennai rain) does cartodb has any api or sample codes to generate a heatmap for that particular search term in the map div. Could any of you explain briefly how this functionality could be done using cartodb api?I have seen the possibility of generating geotagged maps using cartodb.


texts: loading multiple shapefiles using shp2pgsql in a bash script

texts:
I have this bash script (loadShapefiles.sh) to load multiple shapefiles at a time. The script and the shapefiles are stored in the same directory. When I run the bash file with the command of bash loadShapefiles.sh, I have this error message: syntax error near unexpected token '$'do\r'' I appreciate your suggestions to fix this problem!

for f in *.shp
do
    shp2pgsql -d -I -s 2913 $f rlis.`basename $f .shp` > `basename $f .shp`.sql
done

for f in *.sql
do
    psql -d portland -f $f
done


texts:
My test scenario is to confirm the image below, I select the points on the line in order of the direction the line was drawn.

I used the sql:

    Select
  testepoint."id"
FROM
    testepoint
inner JOIN
    teste
ON
  ST_DWithin(teste.geom, testepoint.geom, 1)
WHERE
    teste."id"= 1
ORDER BY
  ST_Line_Locate_Point(teste.geom, testepoint.geom)

I placed the order by to sort the way I want, as the line was drawn

    ORDER BY
  ST_Line_Locate_Point(teste.geom, testepoint.geom)

And the error occurs:

 ERROR: line_locate_point: 1st arg is not a line

This was because my geometry is a MultiLineString.
How can I do this sort of which is a MultiLineString?


texts: XY coords for clicked point QGIS graphical modeller

texts:
I am trying to calculate a bearing from each point in a layer to a single target using the graphical modeller as its part of some other processing I need to do. I can do the bearing calculation (below) but is there a way to get X,Y coordinates from a mouse click and use that as the target?
(atan2((targetx-$x),(targety-$y)))*57.295779513082320876798154814105


texts: Does converting polygon to raster generate more accurate results from the Zonal Statistics in Arcgis 10.2

texts:
I have heard the following for Arcgis 9.1 and wonder whether it is still an issue in Arcgis 10.2.  If so, should I convert the polygon into raster before dissolve or convert the dissolved polygon directly?

http://support.esri.com/en/knowledgebase/techarticles/detail/19184


texts:
I am in the process of learning ArcGIS Desktop, and as I was working through the tutorial provided to me, I was told to create a new layer group by right-clicking on "layers" in a new blank map, and choosing "new layer group" from the menu.

My issue is that it does not work, it does respond (screen twitches) but it doesn't show a new window for it.

This is on ArcGIS 10.5.1, and on a personal laptop.

I am looking for suggestions for a workaround for this issue.


texts: Autocad 2007 raster design loading problem

texts: I try to load Autodesk raster design 2007 from this link https://knowledge.autodesk.com/support/autocad/downloads/caas/downloads/content/autodesk-raster-design-2007-object-enabler.html, after I downloaded and instaled AutoCAD, I loaded the application but without results, I have to load some .sid, .ecw file with this application.


texts: How to keep the ids of way points of a shapefile when I convert it to a GPX file in QGIS?

texts: I tried to convert a shapefile to a GPX file in QGIS, and I do not know how to keep the ids (attributes). What is the solution?


texts: How to remove compression from jpeg2000 NAIP imagery?

texts: I downloaded a 4-band 1m spatial resolution NAIP image from Earth Explorer.  The default image format is jpeg2000. You can see the trees in the attached image appear grainy, more so than what you would expect from a TIFF image. I believe this may be a result of the compression algorithm. To test, I used the Copy Raster tool in ArcGIS to convert the jpeg2000 to TIFF. I specified "NONE" for the compression in the tool environment. Unfortunately, the image appeared the same, although I noticed that there is LZW compression associated with the TIFF. Is there a way to essentially "unpack" a jpeg2000 image to remove all compression? I am open to any software-based approach or programmatic (Python, R, or Matlab) approach.


texts: How to do a Spatial Join after opening a GDB file with Fiona?

texts:
After I successfully open up an ESRI GDB file/folder with Fiona, what open-source Python package should I use to do a spatial join?

I was told Shapely is not the appropriate package for this.

I am completely new to GIS.


texts:
I'm new to QGIS and have tried to create isochrones around points for quite some time.
Is there any option to create isochrones on the GUI in QGIS?

I know how to create them using pgAdmin (using eg underdark's explenation and example (http://anitagraser.com/2011/02/09/creating-catchment-areas-with-pgrouting-and-qgis/)), but I would like to be able to create them using only QGIS.

There is one layer which contains the geocoded addresses around which I need to create the isochrones. It would be great to be able to create different sized isochrones
I've tried to find a plugin or any other option to do this, but nothing seems to work.

I'm using QGIS 2.4.0 on Ubuntu.


texts: Recoloring geotiff for visual presentation

texts:
I have a rather big Geotiff. It is the mosaic of four RapidEye satellite images. They are all a bit bland/grayish. When combined I would like to freshen up their colors so that they are more visually appealing. The imagery will be used as a backdrop for other data.

I have tried right-clicking on the raster in the layers list and selecting Properties. Using Symbology allows for detailed adjustments. But there is no easy way to auto-adjust the image color balance.

Some image software like IrfanView has "auto adjust colors". IrfanView can not open my geotiff image, but this kind of functionality represents what I am  looking for.

To the extent I am doing analysis on the imagery I will of course not mess with the color balance.


texts:
I'm trying to apply this equation through the raster calculator in ArcGIS:

I'm using this syntax:

(  (  ( (Raster) ^ 4 )  - (1-e) *  ( Sky T ^ 4 )  )  / e )  ^ 0.25

The raster calculator create a raster with no values

I'm doing something bad?


texts: Concatenating string data plus numeric data in ArcMap field calculator?

texts:
I am using ArcMap 10.2.2.

I have three attribute fields (Address_No: 1, StreetName: Main, Street_Ty: ST) that I want to merge together. How do I concatenate string data found within the table plus numeric data?

I would like the following result: 1 attribute field - Address: 1 Main ST


texts: 405 Method not Allowed on uploading a shapefile via REST API on GeoServer

texts:
I'm using PHP cURL to upload a Shapefile to GeoServer. It works fine in a development server. But when I deploy to Stage server and I try the same I get:

405 Method Not Allowed
The requested method PUT is not allowed for the URL /geoserver/rest/workspaces/my_workspace/datastores/my_datastore/file.shp.
Apache/2.4.18 (Ubuntu) Server at my.stage-domain.com Port 80

I've tried Google about how to enable PUT requests, but nothing seems to work. I've installed both development and stage servers (both running Ubuntu 16.04) and I didn't had to enable PUT on my development server. It just worked.
Looking around here I've found this question and the answers suggests to register the webapp, but I don't even know how to do that...
Do you have any suggestions? I'm really stuck here...
UPDATE:
As suggested in comments, I replaced the Linux binary version for the geoserver.war file running over Tomcat 8. I've also reinstalled the server so now I have Ubuntu 17.04 and Apache 2.4.27 as in my development server. Unfortunately, the issue remains exactly the same. So I guess it's a web server issue.
My access.log file on Apache2 logs shows:

aaa.bbb.ccc.ddd - - [08/Jan/2018:15:02:27 +0000] "PUT /geoserver/rest/workspaces/my_workspace/datastores/my_datastore/file.shp HTTP/1.1" 405 629 "1" "-"

UPDATE 2:
Some answers and comments suggest to use POST instead of PUT method. I am uploading a zipped shapefile to an existing datastore. This:

/workspaces/{workspaceName}/datastores/{storeName}/{method}.{format}
Uploads files to the data store, creating it if necessary

In the API docs it says clearly that I must use PUT


texts: How can I size the marker (square) to let its scale value be its actual area?

texts:
I want to have a square for each administrative area.

A value tells its area in reality.

For each administrative area the square should be as big as in reality.

So if I zoom in and out it should not change in size.

If I measure the width of the square using the scale of the map, and take it to the square (^2), I end up with the area value used at the beginning.

Proportional Circles may be related to what I try to accomplish, however I need a square, not a circle.


texts: Using ArcPy cursor to concatenate row values?

texts:
I am ArcGIS 10 user and I want to create ArcPy script for simple concatenate row values like this tool Concatenate Row Value that works with ArcGIS 10.

I have three fields in my shapefile id,field1,field2.

Here the image from my shapefile .dbf after sort by id.

If you see the image you can see the id not have unique values.

I want to create the layer with unique id values but I don't want to loose some values from field2 (is text). For field2 I need to use for next conditions.

In Excel if I use some condition like this:

=IF(AND(A3=A2;C3=C2);D2;IF(AND(A3=A2;C3<>C2);D2;IF(AND(A3=A2;C3<>C2);D2&","&C3;C3)))

Then I take new field to like this:

That works in Excel, I want to do with Python or ArcPy or some tool from ArcGIS.

In QGIS, it is easy for me to do with this plugin that does dissolve with stats.

I try to follow this code:

updateFields = ["id", "field2", "field3"]

with arcpy.UpdateCursor("featurelayer", updateFields) as cursor:
    for row in cursor:
        if row[0]+1 == row[0] and row[2]+1==row[2]:
            row[3] = row[2]
        elif row[0]+1 == row[0] and row[2]+1!=row[2]:
            row[3] = row[2]
        elif row[0]+1 == row[0] and row[2]+1!=row[2]:
            row[3] = row[2]+','+row[2]+1
        cursor.updateRow(row)

That looks correct but any time to try to run it I take this error:

 Cursor AttributeError: __exit__

How do I fix the problem?


texts: PyQGIS `iface` object for layer event handling in QGIS plugin testing

texts:
I'm writing a functional test for the QGIS plugin I'm developing. There's a button in the plugin that is enabled only if there is a valid vector layer present. The thing I'm bumping into is the layer event handling from a standard lib unittest isn't triggering the button enable as expected. When I test manually, the button is enabled fine, but from the test it isn't.
My test setup is pretty straightforward:
import unittest

from qgis.core import QgsApplication, QgsVectorLayer

from namari_dockwidget import NamariDockWidget

app = QgsApplication(argv=[], GUIenabled=True)
app.initQgis()

class NamariDockWidgetTest(unittest.TestCase):
    def setUp(self) -> None:
        self.dock_widget = NamariDockWidget()

    def test_dockwidget_layer_selector(self) -> None:
        with self.subTest('When we start the widget, there is no layer set'):
            self.assertTrue(self.dock_widget.isEnabled())

            layer = self.dock_widget.mMapLayerComboBox.currentLayer()
            self.assertEqual(layer, None)

        with self.subTest('So the "Build model" button is disabled'):
            enabled = self.dock_widget.pushButtonBuildModel.isEnabled()
            self.assertFalse(enabled)

        with self.subTest('But when we load a vector data source'):
            layer = QgsVectorLayer(
                path='test/data/amersfoort-centre.gpkg',
                baseName='amersfoort-centre',
                providerLib='ogr')

            self.assertEqual(len(layer.fields()), 10)

            # Load the vector layer
            self.dock_widget.iface.addMapLayer(layer)

        with self.subTest('Then the build button is enabled'):
            enabled = self.dock_widget.pushButtonBuildModel.isEnabled()
            self.assertTrue(enabled)

However the last assertion fails: the button is not enabled once the layer is added through the iface object (that I learned from this excellent question). Any idea on what I'm missing?
To reproduce: the current state of the plugin is in https://github.com/reinvantveer/namari/tree/button_test_issue-alpha2. There's a very helpful Dockerile that can help you get up and running the test within minutes: just run docker build ., but there's also the output of the associated GitHub action build at https://github.com/reinvantveer/namari/runs/2208580233.
EDIT: I now see that I only loaded the UI elements from the NamariDockWidget, but this does not include the Namari class I wrote for the logic. Somehow, I need access to the QGIS iface object from within my unit test, but how does this work?
The plugin logic class begins like
# ...
class Namari:
    def __init__(self, iface: QgisInterface) -> None:

So a lot of the logic depends on the iface object being available. But I have currently no way of strapping a class instance to my test.


texts: Creating custom projection using provided Helmert Parameters and ArcMap?


texts: ArcGIS 10.4.1 export attribute table loses 'shape' column

texts:
I'm running ArcGIS 10.4.1 under Windows 10.

Until a couple of days ago I wasn't having any trouble with exporting attribute tables to csv, editing the resultant files in Excel (Office 365, build 1908) and reloading them. I haven't knowingly changed anything although there has been a MS upgrade to Windows.

Now, every time I download a table it loses the 'shape' attribute with the result that my polygons reload as points. I've tried it with files I've worked with previously and new ones and keep getting the same really frustrating result.

I've also just upgraded to 10.7.1 and the same thing is happening. Short of completely uninstalling and reinstalling I'm stumped about what to do.


texts: How to use Esri data in QGIS using ArcGIS REST API Connector (beta)

texts:
I'd like to include Esri's aerial images in QGIS maps.  I am using QGIS 2.18.3

I have discovered the ArcGIS REST API Connector (beta) http://giswiki.hsr.ch/ArcGIS_REST_API_Connector_QGIS_Plugin

However, I can't work out what URL to place in the box that the plugin launches.  The closest I can get is from this page http://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/3 there is a link in the top right page for 'API reference' http://services.arcgisonline.com/arcgis/sdk/rest/index.html#//02ss0000005t000000

but this page doesn't make it clear what URL to point to.


texts: Print and Export Functionality on web GIS application built with OpenLayers+GeoExt+Extjs+Heron Configuration

texts:
I have used OpenLayers+Extjs+GeoExt+Heron configuration to create a thin mapping client. I have added print icon to the client using heron.options.map.toolbar (type: printdialog, options: {printURL etc}).

I have also the mapfish printing e.g. config.yaml in the geoserver data directory and validation of this installation via xxx/geoserver/pdf/info.json shows the print options as expected. However back to the app, whenever I click on the print icon, I get an error saying:

Error getting print options from server: undefined.

Has anyone experienced this and have solved this problem?


texts: Seeking geocoding with advanced language support?

texts:
I'm building an address disambiguation page and I'm interested in having geocoding results returned in the local language of the country the result's in.

For example, if I make a reverse geocoding request for "Bucharest, Romania", I want the location returned to be "București, România" (the local name). If I have an ambiguous query, I want each result to be returned in the language of the country it's in. Say there was a Bucharest in the U.S, I'd want a search for "Bucharest" to return "București, România" and "Bucharest, U.S.A".

I've studied the Google Maps API, Yahoo! Maps and Bing Maps. Google is the only one that has some language support (see the lang parameter), but only for the entire resultset, as far as I understood. Bing and Yahoo! don't have any language support in their API.

Do you know of a way to do this? Is there any other service that supports this type of functionality?


texts: Extract metadata with arcpy from orthoimagery

texts:
I'm very new to ArcGIS, and GIS is general.  I am trying to build a script that will extract metadata from orthoimagery .tifs. Below is my test code to see if I can get the metadata from one file. It is basically the code from ESRI's docs on arcpy and the metadata toolbox.

import arcpy
from os import path
from arcpy import env
env.workspace = path.normpath("E:/RESOURCES/TEST_RESOURCES/Pasquotank/TIF/data")
dir = arcpy.GetInstallInfo("desktop")["InstallDir"]
translator = dir + path.normpath("Metadata/Translator/ESRI_ISO2ISO19139.xml")
arcpy.ExportMetadata_conversion(path.normpath("E:/RESOURCES/TEST_RESOURCES/Pasquotank/TIF/data/OC6i0_37_000_20793902_20100411R1.tif"), translator)

When I run this is get the following error.

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000816: The tool is not valid.
Failed to execute (ExportMetadata).

I've tried just using the file name as the first parameter in the ExportMetadata_conversion(), and got the same error.

When I run the tool from ArcCatalog, I get results.  So, I'm stuck. I am using ArcGIS 10.0


texts: How to solve the error: query must return columns 'id', 'source', 'target' and 'cost' in postgres 9.3

texts:
SELECT seq, id1 AS source, id2 AS target, cost FROM pgr_kdijkstraCost(
    'SELECT gid, source, target, cost FROM road_network',
            26, array[1055,2013], false, false

);

I am getting the error message:

  (ERROR:  Error, query must return columns 'id', 'source', 'target' and 'cost') when run the above query in postgres....


texts: Get closest (source or target) of an edge from a point in postgres

texts:
I have a Postgres custom function that uses pgr_dijkstra to generate a route working. It works be passing the lat/lon of the start and destination points. I then find the closest segment vertice to that lat/lon point and use that in the pgr_dijkstra call.

snippit:
EXECUTE 'SELECT id::integer FROM trail_split_segs_vertices_pgr
            ORDER BY the_geom <-> ST_GeometryFromText(''POINT('
            || x1 || ' ' || y1 || ')'',4326) LIMIT 1' INTO rec;
    source := rec.id;

The problem is even though I may be moused over a specific segment, the closest vertex might belong to a different segment, since its doing the lookup based on distance.

(the red point is the one I want, but I get the purple one since its closer)

I tried creating a different function to instead find the closest edge to the users lat/lon points rather than the closest vertex. But I have to choose if i'm going to use that edge's source or target.

snippit
SELECT
  id, source, target, geom, ST_Distance(
    geom, ST_GeometryFromText(
      'POINT(-123.036073 49.357152)', 4326
    )
  ) AS dist
FROM
  trail_split_segs
WHERE geom && ST_SetSRID(
    'BOX3D(-123.036083 49.357142, -123.036063 49.357162)'::box3d, 4326
  )
ORDER BY dist
LIMIT 1

The box created by this is very small, so pretty much will always select the correct edge, because I have the users mouse snap to the edges in the web interface.

What I want to do is select the edge source or target that is closest to the lat/lon point where the id of the edge has to be the same.

UPDATE
Zia's response below helped set me on the right track, I modified his query.

WITH
point AS
    (SELECT st_setsrid(st_makepoint(-123.0568528175354, 49.36230206100319), 4326) AS point),
line AS
    (SELECT
    geom AS line_geom
    FROM
    trail_split_segs, point
    WHERE st_dwithin(geom, point, 0.00012)
    ORDER BY geom <-> point
    LIMIT 1)
SELECT
id, ST_AsText(the_geom) AS point_geom
FROM
trail_split_segs_vertices_pgr, line, point
WHERE st_dwithin(the_geom, line_geom, 0.00001)
ORDER BY the_geom <-> point
LIMIT 1;


texts: Is there an ArcGIS Pro Style Manager like ArcMap's?

texts: Is there a way to manage all your styles within ArcGIS Pro? I tried looking under Symbology and Options but cannot seem to find a place to manage my custom styles so that I can reuse them in other projects.


texts: GRASS GIS changes shape of tiff file when imported

texts:
I do have a raster which shows 13 rows and 16 columns when viewing it with Windows Viewer and also when using imread for loading it into python.

I want to work with my dataset in GRASS. My steps:

    r.import input=RGI60-11.00083.tif output=RGI60 --overwrite
    g.region --overwrite raster=RGI60

However, when I now look at the raster in the Map Display, the raster has 12 rows and 22 columns. It seems that every third column is repeated.
I loaded the tiff file to ArcMap, and it displayed correctly. I therefore think I might be having a problem with my grass configuration, but I am unable to find any configuration variable that might be causing this.

Sample data can be downloaded from here:
https://www.research-collection.ethz.ch/handle/20.500.11850/315707
Folder to download: composite_thickness_RGI60-18.zip
And then the first image: RGI60-18.00001_thickness.tif

How could I fix this problem?


texts: Importing analytical data (CSV) with different data types into ArcGIS Desktop to generate concentration heat maps?

texts:
After reviewing all the possible related solutions to my problem, I could not piece together a solution to process and import multiple sets of groundwater analytical data that contains different data types from CSV format (example below) into ArcGIS 10.1 for heat mapping (or generating isoconcentration lines).

Ultimately, I need a CSV file that contains only numeric data types.

As you can see, I have integer, floating point, and non-numeric values (e.g. ND = not detected, NA = not analyzed).

I can substitute NDs with 0's, but I don't know how to handle the NAs since they have no value but cause ArcGIS to interpret my fields as strings.

If I substitute the NAs with no value, then when I export my data to a shapefile, it automatically populates the blank fields with 0's which is not accurate.

I want to find a solution where I don't have to manually add a new field for each constituent and run the field calculator to redefine my data types.

I was hoping to find a ModelBuilder solution, but I couldn't figure out how to avoid the string problem.

I am not opposed to using Python directly, but I am not an experienced user.


texts: I need to know if WMS of geoserver with "localhost" url linked in openlayers or leaflet code will show on public access by other users/clients or will not be shown?


texts:
I'm very new to this so bear with me.

Hopefully a relatively simple query.

I have a city map with property boundaries and I am wanting to filter the boundaries by lot size. i.e just show lots (or polygons) of a certain size, or between two sizes. i.e. lots over 700m2, or lots between 800 and 900m2.

The information I have pulled in doesn't have an area column associated with it so thinking I will need to add my own expression?

Any help would be greatly apprecaited.


texts: Asking user to select file in Python Script tool using parameter?

texts:
In a custom toolbox you can make a parameter and ask user for input where he can type something and used as arcpy.GetParameterAsText(0) in the script in the toolbox. This is good if you want to use address of a folder but what if needs to select a csv file?

One way would be:

User types the folder path and the script will add + file.csv

Example code:

import arcpy
import pandas as pd
file = arcpy.GetParameterAsText(0)
df = pd.read_csv(file)

In this case the user will be asked to select a path when using the toolbox.
To do this he just copies the path of the folder the csv is. Then he must add manually to the path the: ex. \pop2012.csv.


texts: Showing spatial coincidence using QGIS / PostGIS?

texts:
I have a database of fire hydrants (HYD). The positions aren't very exact, the other fields are very useful. I also have a database of GPS'd fire hydrants (GPS). We are manually moving the HYD points to snap to the GPS points, and will then do a join on spatially equal features.

Being manual, this gets tedious, and it can be helpful if a visual indication shows you have actually successfully snapped the HYD to the GPS point. (lets say the HYD point goes from red to yellow, for instance, or changes size/shape, etc)

I've come at this from different angles, but can't seem to find a good way with PostGIS or QGIS to make this visual change happen.

Help? Also, note, we are doing this with approx 40 features types, and some will involve lines touching points. Similar enough concept, though.


texts:
I'm using PostGIS/Tiger for a mapping solution.

For a given longitude, latitude (derived when someone clicks on a point in my Google Map), I want to retrieve the closest point from my table. My table has 5 columns on it:

id, customvalue, the_geom, longitude, latitude

(so I already have a longitude/latitude and geometry for each row)

For a given long/lat, I would like to find the closest row in the database. I achieve it in this way:

    SELECT * FROM mytable
    WHERE ST_Distance_Sphere( ST_Point(ST_X(ST_Centroid(the_geom)),
    ST_Y(ST_Centroid(the_geom))),
    (ST_MakePoint(-73.94039154052734, 40.68779755192939))) <= 300

This works and gives me my results! The 300 I just made up as it seems to give me 2-3 (but not thousands). I then just take the first record and i'm done.

The problem is this query takes about 10-13 seconds to run each time. I would love to speed this up.

I have two possible solutions I would like help on:

This query is overkill and not constructed right because I'm not taking advantage of the fact I already have a geometry, long, and lat on each row in my table, and I should construct the query differently.
The query is actually correct, but I am missing some indexes that could help me out and make things much faster.


texts:
The TPI is calculated with a 150 to 300m annulus in Google Earth Engine.
The underlying Elevation-Data from SRTM does not show such striped.
I did use raw SRTM from GEE for elevation, no reprojecting or resampling

Unfortunatly i cannot share the whole code, the TPI is not depicted well with Map.addLayer()
 The depicted example is an asset export with a Landsat 5 projection from central europe. The image is from Switzerland.

My TPI Function

var TPI = function(elevation, inner, outer){
  var inner = ee.Number(inner)
  var outer = ee.Number(outer)
  var kernelo = ee.Kernel.circle(outer,'meters',false,1)
  var kerneli = ee.Kernel.circle(inner,'meters',false,-1)
  var kernel = kernelo.add(kerneli)//created annulus
  var focalmean = elevation.reduceNeighborhood(ee.Reducer.mean(), kernel)
  var name = ee.String('TPI').cat(ee.String(inner)).cat('_').cat(ee.String(outer))
  var TPI = elevation.subtract(focalmean).regexpRename('elevation',name)
  return TPI
}

are the stripes fixable with some reprojecting? My best Gues: the TPI algorithm highlights the Data pattern from Sensorcollection. Thus all stripes are aligned,


texts:
I am using openlayer3. And using icon as a marker. So my requirement is that when I hover on marker, it should increase in size to show that it is hovered.
Till now I have success to detect the icon and change cursor pointer.
Live example for this can be found at http://codepen.io/anon/pen/Wxzbjv

Javascript Code is as follows :

  var
    sourceFeatures = new ol.source.Vector(),
      layerFeatures = new ol.layer.Vector({
          source: sourceFeatures
      })
  ;

  var map = new ol.Map({
      target: 'map',
      view: new ol.View({
          center: [-5484116.753984261,-1884416.14606312],
          zoom: 16,
          minZoom: 2,
          maxZoom: 20
      }),
      layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM(),
            opacity: 0.6
        }),
         layerFeatures
      ]
  });

          $(map.getViewport()).on('mousemove', function(e) {
            var pixel = map.getEventPixel(e.originalEvent);
            var hit = map.forEachFeatureAtPixel(pixel, function(feature, layer) {
              return true;
            });
            targetStr=map.getTarget();
            targetEle=typeof targetStr==="string"?$('#'+targetStr):$(targetStr);
            if (hit) {
              targetEle.css('cursor','pointer');
            } else {
              targetEle.css('cursor','');
            }
          });

  var
    fill = new ol.style.Fill({color:'rgba(255,255,255,1)'}),
      stroke = new ol.style.Stroke({color:'rgba(0,0,0,1)'}),
      style1 = [
          new ol.style.Style({
              image: new ol.style.Icon(({
                  scale: .7, opacity: 1,
                  rotateWithView: false, anchor: [0.5, 1],
                  anchorXUnits: 'fraction', anchorYUnits: 'fraction',
                  src: '//cdn.rawgit.com/jonataswalker/map-utils/' +                        'master/images/marker.png'
              })),
              zIndex: 5
          })
      ]
  ;

  var feature_start = new ol.Feature({
          geometry: new ol.geom.Point( [-5484116.753984261,-1884416.14606312])
      });
  feature_start.setStyle(style1);
  sourceFeatures.addFeatures([feature_start]);


texts: Problem with open data ArcGisMapServer from SITG (Geneva) in QGIS

texts:
I would like to use this open data DSM data source QGIS from Geneva, Switzerland municipality. Here is the data server

I tried to link it as an ArcGisMapServer source or as WMTS. But when I open the layer properties, QGIS 3.12 (on Linux) cannot determine raster data type.

What am I doing wrong? Thank you very much in advance.

PS I have been struggling with this for a couple of days, sorry. At this point I do not even know if it is a dumb question...


texts: Seeking world city boundary map or data from OpenStreetMap DB?

texts:
I'm seeking world city boundary data(longitude and latitude) with city name from OpenStreetMap DB.(If you have experience with other DB Please share the information to me.)
It is OK if the data is not fully accurate.

I just want to draw world city boundary in the map with city name.

Could anyone experience this kind of data?


texts: Visualisation of "dominant" point in each location

texts:
I had an idea for a kind of visualisation, and I'd like to know if it's been done before, what it's called etc.

Basically it would visualise a single numerical quantity that varies greatly, at relatively few sample points. For instance, population sizes of cities. It would probably also work well for signal strength from radio transmitters.

To calculate the visualisation using a raster scan, given n(i) is the population of city i:

Assign every city a unique colour. (Or, better, a small number using eg 5-color algorithm.)
At every (X,Y), calculate n(i)/d for every city i where d is the distance from the point to i. (Perhaps something like log(n(i))/d would be better.)
Choose the city that has the highest value, and color the point accordingly.

The end result would be that a city with a big population would have a large circle around it, interrupted when it runs into smaller towns. A very small town might be just a small dot within a much larger circle.

Obviously it would look nicer to calculate this using vectors, but I have no idea how.

Is this a known thing? Any libraries, tools, techniques to produce it?

EDIT

Here's the closest thing I can find, using the tip of searching for "influence map". (Unfortunately that term is also used to describe an unrelated meme in the deviantart online art community...):

I'm imagining something with much simpler, more geometric, borders between points though.


texts: python script to retrieve names and scale range from an ArcGIS MapServer url

texts:
i am pretty new at ArcPy.
i need to write a python script to retrieve names and scale range from an ArcGIS Map Service url.
i had a look around json and url2lib but i am lost.

I want to retrieve the layer information (name and scale range) from my ArcGIS Server map service. I am using ArcGIS Server 10.0 SP5.

I try this:

import urllib2
url="http://<myserver-name>/ArcGIS/rest/services/CU/MapServer?f=pjson"
json = urllib2.urlopen(url).read()

and i stuck at parsing tha json objetc to get the layer information i am needing.


texts:
I am using QGIS v 1.7.4 for MAC and I am trying to calculate mean coordinates for a set of points using ftools.

To calculate my result I am going to Vector -> Analysis tools -> mean coordinates and then getting a plot output to a new layer.

However, the point seems to be totally off, even if I filter out outliers in my original points set.

Furthermore, I've read it's possible to calculate both regular and weighted means but don't find this option anywhere.


texts:
I'm trying to set up QGIS Server (2.18.18) with ISS (7.5) on my local machine with Windows 7 Ultimate 64bit.

I managed to get it running with Apache (XAMPP) without any problems, but I need QGIS Server running with IIS. I could really find a tutorial or anything on google. I tried following these instructions (https://domasin.github.io/valgrande/qgisserveriis.html) but never got the positive GetCapabilities request (just the download of the qgis_mapserv.fcgi.exe or Error 500.0 with code 0xc0000135).

Also tried following this one: https://github.com/mapserver/mapserver/wiki/Installing-Mapserver-under-IIS-7.n-with-FastCGI

..with the same results as mentioned above.

What am I missing?


texts: Destination Relationship Not Found in Survey123?

texts:
The goal of my project is to create a road sign inspection survey with the ability to take inspection photos.  I have some limitations, so I can only work with survey123 NOT collector.

I create the survey from a feature service. There are global IDs, everything is editable, attachments are enabled etc etc.

The trouble arises with the photo. Since I am doing the survey on existing geopoints (aka the signs), I can't just add a photo question.  The work around is add the photo question via a "repeat."  For this to work, I created a table and related it to the feature service layer then uploaded to AGOL. I really don't know what's going wrong or where to go from here. I'm having a hard time finding documentation on this topic, so I'm hoping someone can assist me.

I get this error message when I try to publish:
The custom feature service submission url is not compatible with this survey (Destination relationship not found for table)


texts: ArcGIS Online vs JS rendering differences

texts:
I am experiencing some issues working with ArcGIS JS API v4.9. I created a map on my offline software that I imported into ArcGIS Online. The map contains polygons with labels inside them and some arrows. In ArcGIS Online they are rendered correctly.

Now, when I use this map inside JS API, some of the labels are placed incorrectly (they are far from the center, sometimes they are rendered outside their polygon). Arrows have the issue where the arrowheads are not rendered at all, only the base line that marks the direction of the whole arrow. There are no error messages in the console, so I have a hard time debugging this issue.

WebGL for Feature Layers is enabled. I initiate my WebMap like this:

let webmap = new Webmap({
    portalItem: {
        id: "id-of-map-hosted-in-arcgis-online"
    }
});

And my view like this:

let view = new MapView({
    container: "viewDiv",
    map: webmap,
    center: [...],
    zoom: 18,
    higlightOptions: {...},
    showLabels: true
});

I would like to know the source of my issues.

Why is it rendered differently in ArcGIS Online and when using ArcGIS-JS-API 4.9?


texts: Raster calculator error: Invalid Layer: GDAL provider Cannot open GDAL dataset

texts:
I'm attempting to create an NDVI raster using a multiband image from a Parrot Sequoia. It's got only 4 bands instead of 5 so I'm using the following setup for creating the NDVI .tiff:

However, I get the following error message. This does not happen to me when I use the raster calculator for 5 band images like those produced by MicaSense RedEdge. I'm following the procedure written down here on MicaSense's website.

Invalid Layer: GDAL provider Cannot open GDAL dataset ndvi.tif: `ndvi.tif' does not exist in the file system, and is not recognized as a supported dataset name. Raster layer Provider is not valid (provider: gdal, URI: ndvi.tif

This are my layer's properties:


texts: Substitute different equation parameters (e.g. slope) when batch processing raster calculator in ArcGIS Desktop?

texts:
I have 15 multiband rasters for a total of 45 total single-bands (5 dates x 3 sensors x 3 bands per sensor) for which I need to apply a unique calibration function in the form Y=A*[exp(B*X)] to for each band and date. I have a spreadsheet with all the A and B parameter values corresponding to dates and band IDs.

Is there a feasible way for me to apply these functions iteratively in model builder (I'm be willing to learn how to do so in Python as well) to my data? Essentially, how would I tell the raster calculator to substitute the parameters for each date?

My thoughters were to rename all the rasters in the format "Sensor_date.tiff", then try to split them into individual bands using the "Make Raster Layer" tool . I would then feed those temporary outputs into the raster calculator and output the fully-calibrated single-band rasters. Ideally, I also want to feed those outputs into the raster calculator multiple times to determine separate indices.


texts: How to extract biggest Polygon from MultiPolygon in GeoPandas?

texts:
Given a Shapely MultiPolygon, how can I extract the biggest Polygon?

For example, given a MultiPolygon as follows:

from shapely.geometry import Polygon, MultiPolygon

multipolygon = MultiPolygon([Polygon([(0,0), (1,0), (1,0.25), (1, 0.5), (1,0.75), (1,1), (0,1)]),
                             Polygon([(2,0), (1.2,0.25), (1.2,0.25), (1.2,0.75), (2,0.75)])])

How can I get just the largest part?


texts: Whats an appropriate map projection for georeferencing aerial photos?

texts:
I have some aerial photos that were shot from directly overhead a 5x5 mile property. They were then mosaiced to create a continuous overhead tiff of the entire property.

I am now going to georeference in ArcMap using satellite imagery as my base reference. Is there one map projection that would be more appropriate than others for using as reference?

Edit: By appropriate I mean - which projection will results in the most visually similar look to an overhead, aerial photo? Using WGS84 is too wide and stretched.


texts: Convert like raster cell values to a polygon?

texts:
I have a population density raster and would like to create 5 polygons from the 5 natural break (jenks) classification ranges. 1 polygon for each classification range. I'm using ArcMap10 and have all of the Extensions enabled in the ToolBox to my disposal.

Thanks!


texts: Measuring the intensity of raster pixels within buffers

texts:
I need to create a variable ranging between 0 and 1 to compute the intensity of a flood.
I have the raster file for the flood and the layer with the coordinates of the households that have been affected by the flood. I created for each coordinate a buffer with a radius of 5 km. I now need to compute the intensity of the flood in each buffer (i.e., take the number of pixels that are within each buffer and somehow compute its intensity).
I tried with a shapefile of the flood using the "Intersection" command between the buffer and the layer of the flood, but it is not an appropriate measure because I do not have information on the intensity, which is instead given by the pixels of the raster.
Any idea of how to do it?


texts: I can see that mapproxy advertises both image/TIFF and image/GeoTIFF as supported formats, but the output of the latter does not seem to contain geocoding information. Have I misconfigured something, or does mapproxy simply not implement this feature as you'd expect?


texts: Creating rectangles along line using ArcGIS for Desktop?

texts:
I'm trying to create rectangles along a line and can't seem to find a solution anywhere.

I've got a curved line roughly 16km (~10miles) long that I've placed points along in 5m increments.  From there I need to create a rectangular buffer around these points with dimensions of 1600 m length and 400 m width (200 m on either side of the line).

Is there a method that I could use that would allow me to create rectangular buffer around these points but have the rectangle rotate with the direction of the line?


texts:
First, I'm new to TileMill, Ubuntu, and OpenStreetMap so I may be overlooking something basic.

I have a Postgres/PostGIS database with my OSM data. I created a new project in TileMill and added a layer for my OSM data in Postgres.

The query selects everything in planet_osm_line. When I click the Features icon in the Layers box, a grid is displayed with my linear features (roads, bike paths, transmission lines, etc.) To me, this indicates data is being returned from the database to TileMill. The next step is displaying it.

I'm trying to use the roads.mss style from https://github.com/gravitystorm/openstreetmap-carto but nothing is displayed.

When I simplify my styles to the following, the screen color changes but still no linework.

@water-color: #b5d0d0;
@land-color: #f2efe9;

Map {
  background-color: @water-color;
}

What am I missing?


texts: SLD styling in QGIS, elsefilter not working

texts:
I need to figure out how to make elsefilter work when imported into QGIS (2.8.3) via sld.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd" xmlns:se="http://www.opengis.net/se">
<NamedLayer>
<se:Name>dopravni_znacky</se:Name>
<UserStyle>
  <se:Name>dopravni_znacky</se:Name>
  <se:FeatureTypeStyle>
    <se:Rule>
      <se:Name></se:Name>
      <se:ElseFilter/>
      <se:PointSymbolizer>
        <se:Graphic>
          <se:Mark>
            <se:WellKnownName>circle</se:WellKnownName>
            <se:Fill>
              <se:SvgParameter name="fill">#991321</se:SvgParameter>
            </se:Fill>
            <se:Stroke>
              <se:SvgParameter name="stroke">#000000</se:SvgParameter>
            </se:Stroke>
          </se:Mark>
          <se:Size>2</se:Size>
        </se:Graphic>
      </se:PointSymbolizer>
    </se:Rule>
    <se:Rule>
      <se:Name></se:Name>
      <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>umisteni</ogc:PropertyName>
          <ogc:Literal>v chodníku</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>
      <se:PointSymbolizer>
        <se:Graphic>
          <se:ExternalGraphic>
            <se:OnlineResource xlink:type="simple" xlink:href="transport/transport_aerodrome.svg"/>
            <se:Format>image/svg+xml</se:Format>
          </se:ExternalGraphic>
          <se:Size>10</se:Size>
        </se:Graphic>
      </se:PointSymbolizer>
    </se:Rule>
  </se:FeatureTypeStyle>
</UserStyle>

When I use this sld, these styles overwrite each other (I see two overlapping icons), I want the first style to be shown only when the conditions are not met.


texts: Creating new Python Script tool with ArcGIS Notebooks

texts:
I am trying to create a new tool in ArcGIS Pro that is based on a python script written in the Jupyter Notebook interface. I realize I can write the script elsewhere and save it as a *.py file and point the new tool to it, or simply convert the *.ipynb file into a *.py file (as explained here, for example). But these methods force me to have a *.py file for every script I write in Notebooks.
Is there a way to link the new tool directly to the *.ipynb file?


texts: Performing Project_management in batch using ArcPy?

texts:
import arcpy
import os

fc = []

walk = arcpy.da.Walk('C:/Wiley/P1/gis/OSM.gdb')

outdir = "C:/Wiley/P1/gis/HK80.gdb/osm"

for root, dirs, datasets in walk:
    for ds in datasets:
        fc.append(os.path.join(root, ds))

if fc:
    output = os.path.join(outdir, os.path.basename(ds) + "_project")
    hkfc = "C:\Wiley\P1\gis\HK80.gdb\Building_merge_dissolved"
    # Describe the feature class and get its spatial reference
    desc = arcpy.Describe(hkfc)
    spatialRef = desc.spatialReference

    out_coordinate_system = arcpy.SpatialReference('spatialRef.name')

    arcpy.Project_management(fc, output, out_coordinate_system)

I am sure with a trained eye you guys know what to do.
It  came back to me with the following error:

Runtime error
Traceback (most recent call last):
  File "", line 21, in
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 962, in init
    self._arc_object.createFromFile(item)
RuntimeError: ERROR 999999: Error executing function.

What went wrong?


texts:
I have a set of geometries in my database which I used to create a set of rasters (each raster has a unique color). Now I would like to merge all of these rasters  into one but get an error such as :

rt_raster_from_two_rasters: The two rasters provided do not have the same alignment

I tried to resample all of these rasters using ST_Resample(geometry_set, 800, 600) methods but it does not help me.

The numbers 800 and 600 are the width and height values used to generate rasters.

How could I merge all of these rasters into one?

Here is my PostgreSQL script:

SELECT
        layer."SummaryId",
        layer."Name",
        encode(
            ST_AsPNG(
                ST_Union(layer.raster)
            ),
            'base64'
        ) AS pngBase64String
 FROM

(
  SELECT
    s."Id" AS "SummaryId",
    r."Name" AS "Name",

            ST_ColorMap(
              ST_AsRaster(
                  ST_MakeValid (
                      ST_GeomFromText("WktString", 4326)
                  ),
                  800,
                  600
                ),
                1,
                mr."ColorMap"
          --)
       ) AS raster

  FROM "TableA" AS s
  INNER JOIN "TableB" AS r ON r."Id" = s."SomeId"
  INNER JOIN "TableC" AS l ON l."TableAId" = s."Id"
  INNER JOIN "TableD" AS mr ON mr."Id" = l."RangeId"
  WHERE
    s."StartDate" >= '2018-09-10 00:00:00' AND
    s."EndDate" <= '2018-09-10 23:59:59'

  GROUP BY s."Id", r."Name", l."WktString", mr."ColorMap"

) AS layer

GROUP BY
  layer."SummaryId",
  layer."Name"


texts: Grabbing value from field and input into SearchCursor then selecting results for further analysis in ArcPy?

texts:
The data:
PL_Fac -point file - facilities used in service area analysis - contains PNTID in Add_PNTID field
PD - polygon - catchment areas assigned to a facility (there can be many polygons assigned to one point. the polygon file contains a PNTID that corresponds to the PNTID of the facility. The PNTID is also in the feature name
Civic - point file - civic addresses of people who live in the PD polygon areas. The PNTID is in the feature name
The Purpose:
I'm trying to determine if the civic point fall within the service areas they area assigned to. Thus far the civic points do not contain a PNTID but I may need to add one. The script below searches for similar names of the features.
The Problem:
I'm working with the below script to iterate through  over 600 points. It finds PNTIDs that are similar if I fill in the ROW section with the PNTID however, I want the script to automatically grab the PNTID from the Add_PNTID field in the Facilities PL_Fac point file. After the search cursor finds features I then want to select the two found features (one PD_polygon and one Civic_point) and feed them into separate tools. Currently I can't get MakeFeatureLayer from the output from the search cursor. Not sure if some other method might work?

I'm using ArcGIS 10.2 Desktop

import arcpy, os, sys

### #UpdateCursor(dataset, {where_clause}, {spatial_reference}, {fields}, {sort_fields})
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"F:\PL_DriveTime\PL_DriveTime.gdb"
fcPNTID = r"F:\PL_DriveTime\PL_DriveTime.gdb\EG\PL_Facilities"
fcPNTID_lyr = "fcPNTID_lyr"
flds = ('Add_PNTID')

try:

    with arcpy.da.SearchCursor(fcPNTID, flds) as cursor:
        row = cursor.next()
        print(row)
        commonSearchText = ['ROW']
        for nameX in commonSearchText:
            foundFCList = []
            walk = arcpy.da.Walk(arcpy.env.workspace, datatype="FeatureClass")
            for dirpath, dirnames, filenames in walk:
                for filename in filenames:
                    found = filename.find(nameX)
                    if found >-1:
                        foundFCList.append(os.path.join(filename))
                        print foundFCList
except StopIteration:
    print("no rows")
except arcpy.ExecuteError:
    print arcpy.GetMessages(2)
    arcpy.AddError(arcpy.GetMessages(2))
except Exception as e:
    print e.args[0]
    arcpy.AddError(e.args[0])


texts:
I am looking for a way to plot data relating to the number of people in a given location - the location being an ecclesiastical parish. So I need software that will allow me to import a map of parishes in a diocese.  As an ancient Geography graduate (1961 degree) I have no experience of digital mapping.

Do I need to research a grid reference for each parish to add to the people data I already have?


texts:
At the moment I am having some difficulties with calculating a spatial lag in R.
I know how to calculate the lag in space-wide format but am unable to do it in long form, i.e. have repeated observations for the unit of analysis.

Below is some mock data to illustrate what I am trying to do.
Let's start by generating some observations of events that I'm interested in.

# Create observations
pts<-cbind(set.seed(2014),x=runif(30,1,5),y=runif(30,1,5),
       time=sample(1:5,30,replace=T))
require(sp)
pts<-SpatialPoints(pts)

x and y are the coordinates while time represents the time period in which the event takes place. The events need to be aggregated to polygons which is the unit of analysis. In this example the polygons are grid-cells and for simplicity the boundaries are fixed over time.

# Observations take place in different areas; create polygons for areas
X<-c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5))
Y<-c(rep(seq(1,5,1),5))
df<-data.frame(X,Y)
df$cell<-1:nrow(df) # Grid-cell identifier
require(raster)
coordinates(df)<-~X+Y
rast<-raster(extent(df),ncol=5,nrow=5)
grid<-rasterize(df,rast,df$cell,FUN=max)
grid<-rasterToPolygons(grid) # Create polygons

We can plot the data just to get an overview of the distribution:

For space-wide format I would calculate the spatial lag the following way:

pointsincell=over(SpatialPolygons(grid@polygons),SpatialPoints(pts),
              returnList=TRUE)
grid$totalcount<-unlist(lapply(pointsincell,length))
require(spdep)
neigh<-poly2nb(grid) # Create neighbour list
weights<-nb2listw(neigh,style="B",zero.policy=TRUE) # Create weights (binary)
grid$spatial.lag<-lag.listw(weights,
                            grid$totalcount,zero.policy=TRUE) # Add to raster

However, as you can see doing it this way doesn't take into account the fact that the events happen at different moments at time. It just simply aggregates everything to the polygon level. Now I want to calculate this spatial-lag taking into account this temporal dimension so aggregating the data in this case to the polygon-time level.

I wonder whether anyone has a useful suggestion on how this could be accomplished? What is the most convenient way of calculating spatial lags in long format?

I had a look at the spacetime package but was unsuccessful in applying it.


texts: Exporting data frame as PNG after adding Layer File using ArcPy in ArcMap

texts:
My question relates to the answer found at Exporting data frame as PNG in ModelBuilder using ArcMap
I've followed the steps in the answer provided by @PolyGeo, but my output PNG is just a blank image.
I also added a few extra lines of code from ExportToPNG example 2 here to control the size of the PNG, because by default it's 640 by 480 pixels:
https://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-mapping/exporttopng.htm
So my final code is:

import arcpy

lyrFile = arcpy.GetParameterAsText(0)
pngFile = arcpy.GetParameterAsText(1)

mxd = arcpy.mapping.MapDocument("D:/Old C Drive/School/Base Data/Scripts/Untitled.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
addLayer = arcpy.mapping.Layer(lyrFile)
arcpy.mapping.ExportToPNG(mxd, pngFile, df,
    df_export_width=1920,
    df_export_height=1080)
del mxd

My mxd is totally empty, so I'm not sure where I went wrong. I did see in the original question someone mentioned that using a blank document didn't work for them either so they used another mxd and it worked. I tried that but no luck.


texts: Issues with atlas generator and labelling across maps/images

texts:
I have a lot of layers in my project and I am trying to use the atlas generator in the print composer to create a tiled map. Some of the layers have labels which have different configurations with respect to placement.

So far I have found that only the placement offset from centroid with the centroid option whole polygon selected works seamlessly when exporting to tiles. In the following images four neighboring tiles has been merged and the blue lines represents the tile borders:

Offset from centroid:

Around centroid:

Free:

Using perimeter:

The labels looks fine in QGIS when viewing the different areas, but every label is placed automatically to give a good visual appearance based on the current view. In my case the current view should be the entire map and each tile should be rendered according to that view.

Question:
Is it possible to configure the label engine in the atlas generator to calculate the label placements for the entire map instead of each tile separately?


texts: What is the simplest way to host a qgis2web export?

texts: I am new to hosting web maps and would like to know how I can run a qgis2web export and embed a webmap url to sharepoint or share with others. How can this be done open source?


texts: Generate elevation profile graphs for multiple line segments in shapefile?

texts:
I'm familiar with the profile graph feature in 3D analyst and I saw the stack profile tool as well.

However, I'd like to generate elevation profile graphs for 137 line segments.

I don't want to hand draw each segment, is there a tool to do this?

Does anyone have any suggestions?


texts: QGIS 3.10 Add ring to shapefile - "Unknown error"

texts:
I have a polygon shapefile layer with some overlapping polygons, styled using categorised labelling with different styles layered on top of the others.
I would normally use the "Add ring" tool to draw around the top polygon and cut the shape out of the bottom polygon.

This has worked all but tow of the styled features, i.e. I can draw a ring around the trees, cutting a circle out of the hard standing within that layer but if I try doing the same thing and drawing a ring into the grassland (same layer), it doesn't work "unknown error has occurred".

Hard standing and grassland are both set to level 0 while trees are level 2. So I am unsure why it works for one but no the other.

What is this and how do I fix it?


texts: Arcpy Search Cursor Looping through columns to search for specific values

texts:
I am attempting to get a list of row numbers for any rows including either of  two values (-1.797693e+308 and -2147483648) in a list of columns.

I am using an Arcpy search cursor to loop through columns based on a list of headings.

Within each loop I am looking at two columns - the Index column and the column potentially containing either of the two numbers (a) - I then want the Index Column value to be returned.

At the moment the process is successful with uncomplicated numbers (ie 1 or 2) but seems to be struggling with the two numbers I am looking for - potentially as one is actually a string and the other is a minus number.

How could I adjust this so that the loop is able to pick up these numbers / strings?

Field_List = getFieldNames(Shapefile)
Field_Matchers = ['CON',"CAT"]

Search_Fields = [s for s in Field_List if any(xs in s for xs in Field_Matchers)]

for t in Search_Fields:
        FILE = Shapefile
        fields = ["INDEX_COL"]
        fields.append(t)
        cursor = arcpy.SearchCursor(FILE)
        for row in cursor:
                a = ["-1.797693e+308", "-2147483648"]
                #print row.getValue(fields[1])
                if str(row.getValue(fields[1])) in a :
                    print (row.getValue(fields[0]))

fields =["INDEX_COL"]


texts: How to update the Service Definition of a Feature service on ArcGIS Online (AGOL)?

texts:
Here's my workflow:

create a file geodatabase containing lots of data
create an MXD to display this data
> File > Share As > Service within ArcMap
publish to ArcGIS Online*
the dataset and its definitions are uploaded to ArcGIS.com

This is all working fine, and I can access the features via a feature layer in the JSAPI. In the background, ArcGIS Online creates Service Definition and Features entries under My Content.

My question:

I've now updated the symbology in the MXD. How can I get the new symbology definition into ArcGIS Online?

One option is to delete the existing service and recreate it - but this involves uploading the large dataset again, even though there are no changes to the data.
I tried > Share > Save A Service Definition but this fails since the service already exists (even if I delete the existing Service Definition while retaining the existing Feature service)
ArcMap tells me to Overwrite the existing service, but won't this re-upload the data?

Thanks

*aka ArcGIS.com aka ArcGIS for Organisations aka AGOL


texts:
I'm trying to merge two grids with different cell size. One has a cell size of 10m and the other 50m. The 10m one covers a subarea of the 50m one. I know that the most accurate one is the 10m one and hence I would like to keep its values in the overlapping areas, and not do any kind of interpolation. My problem is that the resulting merged grid (after resampling the 50m grid to 10m) has obvious step changes as seen in the picture

Is there a way I can feather out these discrepancies without editing the grids cell by cell? I'm using Qgis and Saga.

Thanks


texts: Missing feature labels stored with QGIS auxiliary storage

texts:
Using QGIS 3.x, I appended new features to a PostGIS feature table for which I had previously generated and modified labels stored in auxiliary storage. No labels are placed for these new features even when the label editor is active.

I don't want to delete existing labels because they are where I want them. What is an efficient way to place labels for these new features, and subsequently to be able to edit them?


texts: Is it possible to convert Polygon to centerline linestring?

texts:
I have a Polygon which I want to convert to centerline linestring.

IMAGE 1

I used the Vector to convert to Polyline but they are not center line as under:

VECTOR TOOL

IMAGE 2

Polyline created by the above tool is not centerline polyline.

IMAGE 3

What I am expecting to get is as under (centerline per road, instead of the lines on the side of the road[Image 3]):

IMAGE 4


texts:
I am a new user of QGIS. I was surprised how easy it is to visualize data from a SQL database in QGIS.
Now I wonder if it is possible that changes in the database can be passed on directly to QGIS. Even better would be the reverse way. Where can I find information about this?

At the moment I use MariaDB 10 and QGIS 2.18.13


texts:
I am having issue with completing download of components using the OSGeo4w 64bit installer. What could be causing this issue and how can I fix it.
I am testing this in Windows 10 20H2 64bit machine, no proxy network.
After contunuing with 'No' option, then during installation another error message box shows as shown below.


texts: Deploying Web AppBuilder Application with php

texts:
I have created my web application with the Web AppBuilder Developer edition, downloaded the zip file and customized the index.html file in my localhost environment. Everything works fine until I change the name to index.php. I am doing this because I have some custom php code that I want to include in the file. I can name it index2.html but *.php seems to not work. Can I do this? and if so, how?

I posted this question here as well with no interest so far.


texts:
I was looking at the data posted on this page
https://datadryad.org/handle/10255/dryad.83330

The readme file says that the coordinate column is given as

  Coordinates (UTM, midpoint): WGS84 32V

and in the main data file some examples of entries are

NM9890048100
NM9920047850
NM9975046050
PM0170046800
PM0110045100
PM0120046100

Coordinate entries start with either NM or PM. Any idea how I can convert these to longitude/latitude?


texts: Fixing overlapping polygons within same feature class using ArcPy?

texts:
I am trying to determine the proper arcpy function or workflow to "Integrate" the borders of overlapping polygons within the same feature class. I am looking for a scripting solution, rather than manual editing in ArcMap, as there will be tens of thousands of these polygons to deal with. (It's important they are not overlapping so I can get the proper area calculations; I don't really care if the vertices are moved, I just want to get rid of the overlap.)

In the picture above, by highlighting one feature I am showing how it is overlapping its neighbour (which is not obvious until you select it or use the Identify tool).

I have played with the Integrate tool, which works for this kind of thing but only between two separate feature classes. Unless I am missing something there ( I have played around with the ranks and XY Tolerance), you can't use Integrtae for this.

Any ideas?


texts: How do I get arcpy.Select_analysis to skip conditions that don't exist

texts:
I have data from 2007 to 2015 of various animals that I'm trying to break up by year and name. I wrote a little for loop that goes through all the items in my attribute table (gps_points).

The problem is, it will search for items that don't exist.
For example lets say I have animal1 for the years 2007 and 2008 only.
The script will make a shapefile for animal1 for year 2007, then 2008, then it will run for 2009 where there is no data and crash.

I want to find a way to skip these periods of no data so after 2008, the code will skip 2009, 2010, ect then move onto the next animal. Is that possible? I can't seem to wrap my head around how an if statement would work since I just want to know if the SQL function will return a result or not before running the select line or skipping. A mixture of using searchcursor and boolean?

for byear in xrange(2007:2016):
    for bname in animal_names:
       sav_loc = "C:\\Desktop\\temp.gdb\\GB_"+str(bname)+"_"+str(byear)
       arcpy.Select_analysis(gps_points, sav_loc, "name={} AND loc_year={} AND Age >= 5".format(str(bname),str(byear))

I'm using Python 2.7 and ArcGIS 10.3.1

extra code:

def unique_values(table, field):
    data = arcpy.da.TableToNumPyArray(table, [field])
    return numpy.unique(data[field])

animal_names = unique_values(gps_points,"name")


texts: Raster correlation: Error (arg6) in raster.modified.ttest (R)

texts:
I am performing pairwise spatial correlation of 6 rasters using the raster.modified.ttest() function in R, as per this very useful post.

I have consistently had the following error with my data;

Error in SpatialPack::modified.ttest(x.var, y.var, sp::coordinates(x[nb[[i]],  :
  NA/NaN/Inf in foreign function call (arg 6)

and yet no errors with any dummy data. The following will yield a result;

s.pix1 <- as(raster(xmn=0,xmx=10,ymn=0,ymx=10,res=1,vals=sample(1:5,100,replace=T)), "SpatialGridDataFrame")
s.pix1[sample(1:100,25)]<-NA
s.pix2 <- as(raster(xmn=0,xmx=10,ymn=0,ymx=10,res=1,vals=sample(1:5,100,replace=T)), "SpatialGridDataFrame")
s.pix2[sample(1:100,25)]<-NA

corr <- raster.modifed.ttest(x=s.pix1, y=s.pix2)

plot(raster(corr), main="spatially adjusted raster correlation")

Very small sample areas of my data would also yield results. However forcing larger contiguous areas of NA data into the dummy data recreated the error;

s.pix1 <- as(raster(xmn=0,xmx=10,ymn=0,ymx=10,res=1,vals=sample(1:5,100,replace=T)), "SpatialGridDataFrame")
s.pix1[c(1:5,11:15,21:25,31:35,41:45)]<-NA
s.pix2 <- as(raster(xmn=0,xmx=10,ymn=0,ymx=10,res=1,vals=sample(1:5,100,replace=T)), "SpatialGridDataFrame")
s.pix2[c(1:5,11:15,21:25,31:35,41:45)]<-NA

corr <- raster.modifed.ttest(x=s.pix1, y=s.pix2)

Error in SpatialPack::modified.ttest(x.var, y.var, sp::coordinates(x[nb[[i]],  :
  NA/NaN/Inf in foreign function call (arg 6)

Could be a bug? @JeffreyEvans provided the excellent answer in the linked post


texts:
I have a Lenovo Miix 2 tablet and I am developing a project with a Java servlet. I want to read GPS data in my project. Can I read this data using Java? and/or HTML5?

update:
In this project, I can't use the Internet to obtain the position of the user. And this is a local web application and the user views the application in a browser.


texts: Is it possible to select every pixel with a certain RGB combination in ArcMap 10.5? For example, from a .tif image I want to select and extract all the red pixels. In this case the RGB codes would be: 251, 118, 121. Is it possible to select these pixels and possibly also expand the range, so that pixels that are slightly darker or lighter would also be selected. Is this at all possible in ArcMap or do I have to use Photoshop (or something similar)?


texts: WebMap that automatically links to Google Drive Spreadsheet

texts: I would like to create a webmap from a contacts spreadsheet stored in Google Drive that will update the map automatically when new data is added to the spreadsheet. I know I can do this with an ArcGIS Online Organisation subscription, but are there any other ways?


texts: How to provide a WFS in multiple projections

texts:
I would like to set-up an editable web feature service (WFS-T) for serving my geospatial data. The data tables are in a PostGIS database on the same server. For security reasons my server admin doesn't allow remote connections to the Postgresql database. So I'm looking at using the WFS service. Clients will view and edit the data in QGIS or using a OpenLayers javascript client.

I succeeded in installing the TinyOWS software and setting up the transactional WFS.

However, I have an extra requirement: All data layers must be available in 3 projections: EPSG:3857, EPSG:3067, EPSG:4326.

Right now each table has a column the_geom with coordinates in the Finnish projection EPSG:3067.

What would you recommend to add the support for the other projections? Does the server like TinyOWS do "on the fly" re-projection on each request, or do I need to add extra geometry columns with coordinates in the other projections to each of my tables?


texts: QGIS how to select entries not in another layer

texts:
I have 2 layers of polygons in QGIS which are layer1 and layer2.
they both have column id.

The goal is to select items in layer1 that not appear in layer2 using unique field id.

For example:

layer 1 has id: 1, 2, 3, 4, 5

layer 2 has id: 1, 2, 3, 4

by using expressions like

id not in layer2

How can I write the expression to have item 5 in layer 1 selected?


texts: Seeking methodology in R similar to ArcGIS Combine?

texts: Is there a way/tool in R that could produce similar output as ArcGIS "Combine" tool does?


texts: PyQGIS 3 Script for PDF Exporting not showing map

texts:
I'm writing this script for exporting PDFs, I managed to get an output PDF as I wanted but it's not showing the Map window. If I get into the Composer, the map is blank but magically refreshes and shows when I click on it.

I tried refreshing the layout after I set extent and size of the map, but it didn't work. I also tried to open the layout through iface.openLayoutDesigner(layout) but that freezed QGIS and forced me to terminate it.

Here is the processing part of my script:

def processAlgorithm(self, parameters, context, feedback):
    input_featuresource = self.parameterAsSource(parameters,'INPUT',context)
    nomenclatura = '1.1.E.3.43' #str(input_featuresource)

    project = QgsProject.instance()
    QgsExpressionContextUtils.setProjectVariable(project,'Manzana_actual',nomenclatura)

    pScope = QgsExpressionContextUtils.projectScope(project)
    tx = pScope.variable('Manzana_actual')

    layer = project.instance().mapLayersByName('Mzs,Qt,Chs')[0]
    iface.setActiveLayer(layer)

    query = '"EJIDO" = ' + tx.split('.')[0] + ' AND "CIRC" = ' + tx.split('.')[1] + ' AND "RADIO" = \'' + tx.split('.')[2] + '\' AND "CC" = ' + tx.split('.')[3] + ' AND "MZNA" = \'' + tx.split('.')[4] + '\''
    layer.selectByExpression(query)

    iface.actionZoomToSelected().trigger()

    projectLayoutManager = project.layoutManager()
    layout = projectLayoutManager.layoutByName("Manzanero")
    #iface.openLayoutDesigner(layout)

    vg = layout.itemById('ventanaGrafica')

    vg.setExtent(iface.mapCanvas().extent())
    vg.attemptResize(QgsLayoutSize(20.85,18.85,QgsUnitTypes.LayoutCentimeters))
    layout.refresh()

    base_path = os.path.join(QgsProject.instance().homePath())
    pdf_path = os.path.join(base_path, str(tx.replace('.', '_') + '.pdf'))

    exporter = QgsLayoutExporter(layout)
    exporter.exportToPdf(pdf_path, QgsLayoutExporter.PdfExportSettings())

    return {'OUTPUT': pdf_path}

I took out the comments cause they were not written in English.


texts:
How can I perform a spectral unmixing fusion ?

The process consists on using a fraction image of a sensor "A" to obtain the reflectance image of a sensor "B", basically the inverse processes of spectral unmixing.

The spectral unmixing assumes that reflectance at each pixel of an image is a combination of the reflectance of each material (endmember) present within the pixel. A linear spectral unmixing, uses endmemembers in order to decompose a reflectance source spectrum into a set of given endmember spectra.

I succesfully used Envi or SPRING softwares in order to perform the spectral unmixing. Although, I don't know how to do the inverse process: use the fraction images to generate the reflectance image, since the softwares I used did not provide a transformation equation.


texts:
Does any version of QGIS support GDB 10.5? I have a GDB 10.5 file that I need to open in QGIS and I can't. I'm using QGIS 3.6.

The file is named Virginia_Parcel_Dataset_2018Q4.gdb -- it can be downloaded here: https://vgin.maps.arcgis.com/home/item.html?id=29627d7c051a47dc8ce71b4484531ab3


texts: Is there any difference between gdal_proximity.py and ArcGIS Euclidean Distance?

texts:
The gdal_proximity.py description says:

  The gdal_proximity.py script generates a raster proximity map
  indicating the distance from the center of each pixel to the center of
  the nearest pixel identified as a target pixel. Target pixels are
  those in the source raster for which the raster pixel value is in the
  set of target pixel values.

While, regarding Euclidean Distance, ArcGIS documentation says:

  Conceptually, the Euclidean algorithm works as follows: for each cell,
  the distance to each source cell is determined by calculating the
  hypotenuse with x_max and y_max as the other two legs of the triangle.
  This calculation derives the true Euclidean distance, rather than the
  cell distance. The shortest distance to a source is determined, and if
  it is less than the specified maximum distance, the value is assigned
  to the cell location on the output raster.

Is there any real difference between these algorithms?


texts: Simulate a mouse-click inside an ArcGIS Online JavaScript Application

texts: I have a JavasScript application consuming an ArcGIS Online Web Map.  I want to open a default popup at a specific coordinate.  Is there a way to simulate a mouse-click at that point, or a better alternative?


texts:
I want to split a line haveing zigzag(different angles) but bthere are many vertics.
I have an option in Arc GIS 10 to split from vertics but i need to split only where the angle or dimension is changed


texts: Matching IDs between point layer and network layer when doing network analysis using "Service area"

texts:
Update: It works!
The attempted solution shown below works in QGIS 3.10.4 LTR.
But in QGIS 3.12.3 it failed in the processing step "qgis:serviceareafromlayer" as some of the features were not added correctly to the output.

Initial question
Using QGIS 3.12 I have run Service area (from layer) to create service area from points checking the whole network.
Background: Each point have an line-identifier corresponding to a line-identifier in the network.
I would like to create the service area only using corresponding identifiers between the point layer and the network layer.
So in a junction where the a different line-identifier connect to the network, I only want the service area to continue on the same line-identifier as observed in the point layer. I have many points and many line-identifiers, and doing it manually would be a daunting task.
Is there a way to get "Service area" from point layer to network where line-identifier-point is the same as the line-identifier-network?

Solution attempted
I have made an attempt at solving this using my limited knowledge of python, and it works with limited number of IDs.
I am concerned that the method would be slow running on larger datasets (both network and point layer), and wonder if there are ways of optimizing the code?
Attempted solution:
from qgis import processing
#Setting input and output
InnPnt = 'c:/GIS/pntFile.shp'
InnNet = 'C:/GIS/netFile.shp'
final = 'C:/GIS/output/final.shp'

# Finding unique ID that matches both inputs
pnts = iface.addVectorLayer(InnPnt, '', 'ogr')
idx = pnts.fields().indexOf('ID')
UniquePnts = pnts.uniqueValues(idx)
net = iface.addVectorLayer(InnNet, '', 'ogr')
idz = net.fields().indexOf('ID')
UniqueNet = veg.uniqueValues(idz)
UniqueID = set(UniquePnts) & set(UniqueNet)
QgsProject.instance().removeMapLayer(pnts) #cleaning up environment
QgsProject.instance().removeMapLayer(net)

# Running extract features where matching ID, running ServiceArea and saving data
NetById = None
for ID in UniqueID:
    PntID = 'ID'
    PntExpr = """"{0}" = '{1}'""".format(PntID,ID)
    ParamPnt = {'INPUT': InnPnt, 'EXPRESSION': PntExpr, 'OUTPUT':'memory:'}
    TempOutPnt = processing.run("qgis:extractbyexpression",ParamPnt)['OUTPUT']
    NetID = 'ID'
    NetExpr = """"{0}" = '{1}'""".format(NetID,ID)
    ParamNet = {'INPUT': InnNet, 'EXPRESSION': NetExpr, 'OUTPUT':'memory:'}
    TempOutNet = processing.run("qgis:extractbyexpression",ParamNet)['OUTPUT']
    ParamService = {'INPUT':TempOutNet, 'START_POINTS':TempOutPnt, 'TRAVEL_COST': '500','OUTPUT_LINES':'memory:'}
    TempSa = processing.run("qgis:serviceareafromlayer",ParamService)['OUTPUT_LINES']
    if NetById is None:
        NetById = TempSa
    else:
        features = []
        for feature in NetById.getFeatures():
            features.append(feature)
        TempSa.startEditing()
        data_provider = TempSa.dataProvider()
        data_provider.addFeatures(features)
        TempSa.commitChanges()
        NetById = TempSa
_writer = QgsVectorFileWriter.writeAsVectorFormat(NetById,final,'utf-8',driverName='ESRI Shapefile')

Problem in the code
I tried running it on many IDs and a larger input layers, I think there might be a problem with appending data, as I would get several features that were empty except for the first two attributes. In the first one it would say "line" and in the next there would be a coordinate, the rest of the attributes says "NULL".
This seems to be an issue with QGIS 3.12.3, as the attempted solution executed without errors in QGIS 3.10.4 LTR


texts: Calculating cost distance in Google Earth Engine?

texts:
I need to calculate cost distance for each cell in my raster (bounding box: Sumatera island) to points (palm oil mills), with slope as its cost. The cost layer (slope) has float data type (val: 0-90), and I can't seem to get the right code to make the variable cost calculation works.

// Add features
var sumatera = ee.FeatureCollection("users/putraditama/sumatera"); // Sumatera island boundary (polygon)
var po_mills = ee.FeatureCollection("users/putraditama/po_mills"); // palm oil mills (points)

// Add raster
var elev = ee.Image("USGS/SRTMGL1_003");

// Clip elev to Sumatera boundary
var elevsmrt = elev.clip(sumatera);

// Get slope
var slopesmrt = ee.Terrain.slope(elevsmrt);

// Create a source image (palm oil mill=1, others=0)
var sources = ee.Image().toByte().paint(po_mills, 1);

// Mask the sources image with itself.
sources = sources.updateMask(sources);

// The cost data is generated from slope
var cover = slopesmrt.select(0);

// Variable cost
var cost =
  cover.eq(0-10).or(cover.eq(11-20)).or(cover.eq(21-30))
      .multiply(1).add(
  cover.eq(31-40).or(cover.eq(41-50)).or(cover.eq(51-60))
      .multiply(2).add(
  cover.eq(61-70).or(cover.eq(71-80)).or(cover.eq(81-90))
      .multiply(3)));

// Compute the cumulative cost to traverse the land cover.
var cumulativeCost = cost.cumulativeCost({
  source: sources,
  maxDistance: 1000 * 1000 // 1000 kilometers
});


texts: Cannot rollback changes when adding features to QGIS layer

texts:
I am testing a script that add features from one *.shp file to another. But also I need to wait for user to save or discard changes. But it seems that adding features  from script using dataProvider() automatically commit changes without waiting for user to do it.
Here is my code

layer_first = QgsProject.instance().mapLayersByName(u'layer_first')[0]
layer_second = QgsProject.instance().mapLayersByName(u'layer_second')[0]
feat = [n.geometry() for n in layer_first.getSelectedFeatures()][0]
f_pnt = QgsFeature()
f_pnt.setGeometry(feat)

layer_second.startEditing()
layer_second.dataProvider().addFeatures([f_pnt])

I tried:

layer_second.rollback(True) # doesn't work

qgis.utils.iface.actionRollbackEdits().trigger() # also doesn't work

with edit(layer_second): # works only from edit session which is not preferrable
    layer_second.dataProvider().addFeatures([f_pnt])


texts: Adding button to move map Direction in Leaflet?

texts:
I want direction button in my Leaflet map so user can move map in mobile or desktop browser easily.

I am new and do not know how to start. Please guide me.


texts: What are my hosting options for an app that needs tilestashe along side a django app?

texts: I'm building a GIS webapp that uses tilestache to serve MBTiles and a django app to serve geoJson. I'm having trouble understanding exactly what kind of setup I'd need to host this. I'm new to server admin and setup, so something heroku-ish would be nice, but I'd be willing to learn if the solution involved more work.


texts: PyQGIS Processing - wrong number of parameters

texts:
I'm facing a very strange issue from the module Processing using it as a standalone (not in the Qgis python console). I'm trying to create voronoi polygons from a points csv file. I can load my file without issue but the Processing.runAlgorithm always fail with " wrong number of arguments ". There are only 3 arguments described in the help.

#create uri
uri = "file://{}?delimiter={}&xField={}&yField={}".format(os.path.abspath(fichier.name),",","x","y")
#open file
layer = QgsVectorLayer(uri, "tmp", "delimitedtext")

# check
if not layer.isValid():
    print "Layer didn't load properly."
    print "Check your filename and the file integrity."

# UTF-8
layer.setProviderEncoding('UTF-8')

#feedback
print("File successfully loaded as " + layer.name() + " layer.")

Processing.runAlgorithm("qgis:voronoipolygons",layer,0.0, "test2.shp")

Traceback :

File successfully loaded as tmp layer.
Error: Wrong number of parameters
ALGORITHM: Voronoi polygons
    INPUT <ParameterVector>
    BUFFER <ParameterNumber>
    OUTPUT <OutputVector>

but it works from the Qgis python console :

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()
processing.runalg("qgis:voronoipolygons",layers[0],0.0,"test.shp")

Output :

{'OUTPUT': 'test.shp'}

Even stranger, when i'm trying to add an unknown argument like :

Processing.runAlgorithm("qgis:voronoipolygons",layer,0.0, "test2.shp",0.0)

it returns :

File successfully loaded as tmp layer.
Error: Wrong parameter value: test2.shp

With adding 2 unknown arguments :

Processing.runAlgorithm("qgis:voronoipolygons",layer,0.0, "test2.shp",0.0,0.0)

TB :

File successfully loaded as tmp layer.
Error: Wrong number of parameters
ALGORITHM: Voronoi polygons
    INPUT <ParameterVector>
    BUFFER <ParameterNumber>
    OUTPUT <OutputVector>

Looks like there is a forth argument ... Any idea ?


texts: How to save aliases and metadata in shapefile in QGIS

texts:
Is there a way to save the aliases in a QGIS shapefile?

When I save the shapefile in a different name or open it in a different QGIS project the entered aliases have disappeared.  The same happens for metadata.

Is this a bug?


texts: Getting a single point from a Mapbox polygon layer

texts:
I want to render a polygon source in Mapbox as a symbol or a single point.

I currently have the following code:

  map.addLayer({
   'id': 'polygon_layer',
   'source':
      {
        'type': 'vector',
        'url': 'mapbox://user.source'
      },
   'source-layer': 'Polygon_Layer',
   'type': 'circle',
   'paint':
      {
        'circle-color': '#EEEEEE',
      }
  } );

This returns a circle at every vertex of the polygon. Is it possible to filter the source data to get just one point, or even a centroid location?


texts:
I am producing some maps to include in a publication. However, I am not able to include a graticule on my plot such that it resembles that in the image shown below.

The image is projected but still has regular longitude and latitude gridlines. I am using the levelplot function for my mapping. Here is the code so far (insights from @Oscar Perpiñán):

library(raster)
library(rasterVis)
library(graticule)
rr <- getData('alt', country='CAN', mask=TRUE)

## Here is where the graticule routine starts
crs.longlat <- CRS("+init=epsg:4326")
ras = projectRaster(rr, crs = ('+proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675 +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80 +units=m +no_defs'))

prj <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

extLL <- projectExtent(r, crs = crs.longlat)

lons <- pretty(c(xmin(extLL), xmax(extLL)))
lats <- pretty(c(ymin(extLL), ymax(extLL)))
## optionally, specify the extents of the meridians and parallels
## here we push them out a little on each side
xl <-  range(lons) + c(-0.4, 0.4)
yl <- range(lats) + c(-0.4, 0.4)
## build the lines with our precise locations and ranges
grat <- graticule(lons, lats, proj = prj,
                  xlim = xl, ylim = yl)
## Labels
labs <- graticule_labels(lons, lats,
                         xline = lons[2],
                         yline = lats[2],
                         proj = prj)
labsLon <- labs[labs$islon,]
labsLat <- labs[!labs$islon,]

levelplot(ras,margin=F,col.regions=col,contour=F,
          par.strip.text=list(cex=0),
          par.settings=list(panel.background=list(col="gray93"),axis.line=list(lwd=1.1), strip.border=list(lwd=1.1),
                            layout.heights=list(xlab.key.padding=-1.5)),
          cex=0.8, scales = list(x=list(draw=TRUE), y=list(draw=TRUE)),
          xlab=list(label=NULL,cex=1.5),
          ylab=list(label=NULL,cex=1.5))+
 ## and the graticule
 layer(sp.lines(grat)) +
 layer(sp.text(coordinates(labsLon),
               txt = parse(text = labsLon$lab),
               adj = c(1.1, -0.25),
               cex = 0.6)) +
 layer(sp.text(coordinates(labsLat),
               txt = parse(text = labsLat$lab),
               adj = c(-0.25, -0.25),
               cex = 0.6))

The problem is that the lonlat labels are not in the conventional format as in the image. I need to keep the projection of ras

Edzer provides a nice solution (curved text label demo) here add graticule to map via sp but I am unable to implement the approach using levelplot

How do I add nice gridlines to the map?


texts: ArcMap Date values lose precision after copying feature class from fGDB to SDE GDB

texts:
I have some data in a file geodatabase feature class including some date fields.  The precision on the date fields is to milliseconds and this is important.

When I copythe feature class from the fGDB to an SDE enterprise geodatabase (MS SQL Server 2014), the date precision appears to have some odd rounding applied to it.

Here's what the difference looks like:

>>> with arcpy.da.SearchCursor("Form_1_fGDB", ['EditDate']) as rows:
    for row in rows:
        print row[0]
2018-12-03 21:22:57.565000
2018-12-03 21:24:59.411000
2018-12-03 21:29:24.494000
2018-12-03 21:32:23.298001
2018-12-03 21:36:04.656000
2018-12-03 21:42:11.316000
2018-12-03 21:47:45.409001
2018-12-06 23:40:23.387000
2018-12-06 22:54:35.083000
>>> with arcpy.da.SearchCursor("Form_1_SDE", ['EditDate']) as rows:
    for row in rows:
        print row[0]
2018-12-03 21:22:58.000001
2018-12-03 21:24:59.000001
2018-12-03 21:29:24
2018-12-03 21:32:23
2018-12-03 21:36:05
2018-12-03 21:42:11
2018-12-03 21:47:45
2018-12-06 23:40:23
2018-12-06 22:54:35.000001
>>>

It looks like they are either rounded DOWN to the nearest whole second, or rounded UP to the whole second plus 1 millisecond.

I get the same result whether using:

right-click and copy/paste in the catalog pane
arcpy.FeatureClassToFeatureClass_management()
arcpy.Append_management()

Why is this happening, and how can I preserve date value precision when copying data from an fGDB into an SDE GDB?


texts: Can I use ogr from OSGeo4W in python env?

texts:
I use QGIS and therefore also OSGeo4W.
At the moment I am writing a script that use ogr in python. I need a virtual environment for because it is going to be deployed on a server (which also has OSGeo4W). After creating a virtual env (python -m venv venv) then I cannot import ogr as I would do from the OSGeo4W shell from osgeo import ogr. I am a windows user, so is the server.

I have been reading many post about issues using/installing ogr in python but at the moment I am completely confused, I . I am not a very experienced python user, and setting up the dependencies appear to be the hardest bit yet.
This post is similar to what I want to do Use OGR and GDAL in Python outside OSGEO4W shell but the link in the post to binaries include mapserver which I am not interested in. Any advice?


texts: Populating ComboBox with unique features in selected attribute field using PyQGIS

texts:
I am creating a plugin in QGIS where I want to read a Vector layer --> Select attribute column --> list unique values within the chosen attribute column in the selected attribute from the column.
I have been successful with the two first steps, but not with the third.
This is the code used so far, but I cannot get further than this.
self.layerComboManagerPoint = QgsMapLayerComboBox(self.dlg.comboBoxVector)
self.layerComboManagerPoint.setCurrentIndex(-1)
self.layerComboManagerPoint.setFilters(QgsMapLayerProxyModel.PolygonLayer)
self.layerComboManagerPoint.setFixedWidth(175)

self.layerComboManagerPointField = QgsFieldComboBox(self.dlg.comboBoxField)
self.layerComboManagerPointField.setFilters(QgsFieldProxyModel.AllTypes)
self.layerComboManagerPoint.layerChanged.connect(self.layerComboManagerPointField.setLayer)

I have tried googling, but I might just not have a reasonable word to google on.


texts: Unfinished Layer with QGIS Contour Plugin

texts:
I am using QGIS 3.4.13 Madeira on Windows 10

I have an issue with the QGIS Contour Plugin. I have two point layers, they contain about 45'000 points in the exact same locations. The only differences between the layers are the attributes. I now want to create filled contours using the Contour Plugin.

With the exact same settings, the process works perfectly fine for one of the point layers. When I try the second one, the resulting temporary layer is left in editing mode and cannot be saved. It also displayes unfinished (see first and second screenshot below). There are no error messages showing during processing. The only message I get is when I try to save the layer (second screenshot). It seems there is a problem with one of the resulting layers.

I have already checked the following:

The point layers do not contain any points with missing values;
I have used the same classification settings for the contouring;
For some odd reason, the problem does not occur when I classify another attribute in the problematic shapefile...
The problem persists after a restart of QGIS.

Does anyone have an idea what could be the issue? Or how to fix it?

Update 24.01.2020: The plugin always seems to be creating multipolygons. Saving to polygon works fine, if the plugin has worked correctly.
In my case, the plugin failes to create all contour categories (in my example the dark blue category is missing...), thus the saving is also problematic.


texts: Reduce raster resolution and determine pixel inclusion based on % pixel overlap

texts: I'm trying to reduce the resolution of a raster from 1m to 50m using a pixel aggregation (average) method, but I need to set a threshold for inclusion (e.g. 75% of the 1m pixels must fall within the 50m pixel in order for the 50m to be included in the new raster).  How to go about doing that? I've looked into the resize tool in ENVI, the congrid tool in IDL, gdalwarp in GDAL, and the aggregate tool in ArcGIS, and haven't figured out a solution.  Any thoughts?


texts: Analyzing building proximity to restaurants using ArcPy?

texts:
I am trying to make a script that will go through buildings in a shapefile row-by-row and do a count of points within .25, .5,...1 mile of it. The other points are restaurants.

import arcpy
arcpy.env.workspace = "C:\Users\zafard\Desktop\Building.gpd"

#need to do 4 iterations for the distance radii
buffDist = [0.25,0.5,0.75,1]
for Dist in buffDist:
    print ("FOR A %s MILE RADIUS: " % (Dist))
    inTable = "C:\Users\zafard\Desktop\Building\XYBuild.shp"

    #set up cursor in table
    rows = arcpy.SearchCursor(inTable)
    row = rows.next()

    #set up loop to analyze every row of the table
    While row:
        print row.OBJECTID

        #use select by distance tool
        arcpy.SelectLayerByLocation_management("FRest","WITHIN_A_DISTANCE", row, Dist, "NEW_SELECTION")

        #count the hopefully selected points for the layer.
        count = int(arcpy.GetCount_management("FRest"))

        print count
        row = rows.next()

I am now getting an error with running the tool in the while loop:

  Object: Error in executing tool


texts:
I would like to use some GRASS functionality from a Python programme. With GRASS 6 this required the set up of various environment variables, as detailed in the GRASS Wiki.

I am trying a similar setup with GRASS 7, as in this example:

import os

os.environ['LC_ALL'] = "en_GB.UTF-8"
os.environ['GISBASE'] = "/usr/lib/grass70/"
os.environ['PATH'] = "$PATH:/usr/bin:$GISBASE/bin:$GISBASE/scripts"
os.environ['LD_LIBRARY_PATH'] = "$LD_LIBRARY_PATH:$GISBASE/lib"
# for parallel session management, we use process ID (PID) as lock file number:
os.environ['GIS_LOCK'] = "$$"
# path to GRASS settings file
os.environ['GISRC'] = "$HOME/.grass7"
os.environ['PYTHONPATH'] = "$PYTHONPATH:$GISBASE/etc/python"

import grass.script as grass
env = grass.gisenv()
print env
#List all environment variables
for key, value in os.environ.items(): print key, value
#Print region extent
r = grass.read_command("g.region", flags='p' )
print r

But the end result is an ImportError:

Traceback (most recent call last):
  File "grass_test.py", line 13, in <module>
    import grass.script as grass
ImportError: No module named grass.script

The basic sanity check shows the script library in the expected folder:

$ ls -la /usr/lib/grass70/etc/python/grass
total 44
drwxr-xr-x 9 root root 4096 Mai  3 08:44 .
drwxr-xr-x 3 root root 4096 Jan  7 08:16 ..
drwxr-xr-x 2 root root 4096 Mai  3 08:44 exceptions
drwxr-xr-x 2 root root 4096 Mai  3 08:44 imaging
-rw-r--r-- 1 root root  151 Mai  1 21:46 __init__.py
-rw-r--r-- 1 root root  363 Mai  3 08:44 __init__.pyc
drwxr-xr-x 2 root root 4096 Mai  3 08:44 lib
drwxr-xr-x 2 root root 4096 Mai  3 08:44 pydispatch
drwxr-xr-x 9 root root 4096 Mai  3 08:44 pygrass
drwxr-xr-x 2 root root 4096 Mai  3 08:44 script
drwxr-xr-x 2 root root 4096 Mai  3 08:44 temporal

What else must be added to this environment setup?

Update: I tried the setup suggested by Gene, but it produces the same ImportError exception:

import os

os.environ['GISBASE'] = "/usr/lib/grass70/"
os.environ['PYTHONPATH'] = "${PYTHONPATH}:$GISBASE/etc/python/"
os.environ['LD_LIBRARY_PATH'] = "$LD_LIBRARY_PATH:$GISBASE/lib"
os.environ['GIS_LOCK'] = "$$"
os.environ['GISRC'] = "$HOME/.grass7/rc"

import grass.script as grass


texts: Postgis using st_geometry linestring on column

texts:
I'm trying to create a linestring based on a value in a column in PostGIS. My code looks like this:

create view left_span as
SELECT st_geomfromtext('LINESTRING(line_att.wkt_left)', 27700) as geom
FROM line_att

It creates the view but without any results.

When I try to select just LINESTRING:

SELECT st_geomfromtext('LINESTRING(line_att.wkt_left)', 27700) as geom
FROM line_att

I get this error:

ERROR:  parse error - invalid geometry HINT:  "LINESTRING(li" <--
parse error at position 13 within geometry
********** Error **********

My column structure looks like this:

wkt_left
608499.156438 127925.765696, 608583.734491 127621.137605
609186.961722 128406.984816, 609287.145785 128550.322806

Calling just one row works fine.


texts: Fixing bad topology on polylines with QGIS

texts:
I've been sent a shapefile containing 12 layers of some city's road system. Each layer depicts a portion of the road system, and all of them partially overlap on the map.

I'd like to work on the sum of these 12 layers in order to have the broadest extent of the road system. The problem is that the person who digitized the analog maps to create the shapefile didn't use the snap feature, and the line segments and their vertices don't exactly overlap where they should.

How would one sensibly deal with this kind of data? Is there a way to automatically correct these messed up topologies?


texts: CRS from a planar projection around a point (pref. in Geopandas)

texts:
Is there a canonical way of projecting a set of points (for example, originally given in WGS84 latitudes/longitudes) into a planar approximation centered on a particular (latitude, longitude) pair and measured in km?

I've basically being doing some changes of coordinates by hand to get ball- and hexagon- shaped neighborhoods with a measurements given in km, but it seems simpler to have a local distortion and deal in planar (km x km) coordinates.

I'm not including any equations because I'm hoping this is a standard "named" CRS family. The Geopandas documentation claims it accepts proj.4 strings.


texts:
I have a series of satellite images in 4 bands (red, green, blue and NIR). What is the easiest way to estimate what percentage of each image is covered by clouds? The image is less than 1x1 km with 3-5 m resolution.

As far as I understand, fmask does not work here (it requires more bands), but maybe there is some simple heuristic (preferably one that can be calculated just on a pixel-by-pixel basis)?


texts: Creating domain with thousands of values for ArcGIS feature class?

texts: I am creating a feature class for the ArcGIS Collector App and want to create domain table for an attribute field.  The problem is that I have 6335 possible values that could be used for the field.  Is it possible to create this table externally and then import into domain table?


texts: Receive 'handle bad layers' and 'invalid data source' message when importing project from QField back to QGIS

texts:
using QGIS 3.14.1 and Qfield 1.5.3 with a Samsung Galaxy TAB A (2016)
I have created a Project in QGIS for use in Qfield. I am able to open the project in QField on the tablet and edit the layers with no problem. However when I import the project back to the laptop and open in QGIS I get a message to handle bad layers for any layers I have edited (the rest are fine) and then when I navigate to the source of the layer I get an 'invalid data source' error and the layer won't open.
This happens whether I use Qfield sync to package the project or not.
File paths are set to relative.
I have tried using QGS rather than QGZ - same outcome.
I have made sure versions of QGIS and Qfield are up to date.


texts: DEM error when trying to do topographic correction in GRASS

texts:
I'm new to GRASS gis and I'm trying to do a topographic correction (i.topo.corr) of 9 ASTER bands using SRTM dem. I keep getting this error message:

  ERROR: Illumination model is of CELL type

Is this error referring to my elevation raster, and if so how do I change the type to the appropriate one?