texts: Is there a public API that can be accessed with Python and will geocode Points of Interests?

texts:
I recently downloaded tagged hospital data from OpenStreetMaps.org using the API.  Some of the data have a 'node' tag which are attributed with a lat/long while other data are tagged as 'way' format without an address attribute.  This leaves me with several hospitals that only have a 'name' attribute and no address/lat-long geolocation.

I am wondering if there are public APIs that can be accessed using Python and allow me to pass just a hospital name like, 'CityA Hospital', and either get an exact address or a lat/long.


texts: How to assign raster layer's colors to polygons in QGIS?

texts: I have a georeferenced raster file showing parcels of a town with different colors (yellow for the category single-household residential, green for agricultural areas, etc.). I also have the shapes of the parcels as a separate vector layer. I would like to assign the color codes of the raster layer to the vector layer as an attribute.


texts: How to extract tiles corresponding to a geometry from a dataset?

texts:
I'm using the GEE API in a Jupyter notebook.
I'd like to use the following dataset to carracterise the tree cover of a specific area. This dataset is unavailable on Gee  but is available for download. The map is cut in small tiles of 10x10° and there top-left corner position. I'd like to extract them programmatically on a specific area (here a country).
I created a geometry and extracted all it's coordinates :
country_code = 'CG'
rdc = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(ee.Filter.eq('country_co', country_code))

def get_coords(coordinates, array_coord=[]):
    """get all the coordinates and set them in a single table without knowing in advance the depth of the tab"""
    if isinstance(coordinates[0], float):
        array_coord.append((coordinates[0], coordinates[1]))
        return array_coord

    for item in coordinates:
        get_coords(item,array_coord)

#transformer cette horreur en tableau 1D avec toutes les coordonnées
#les points extremes sont forcements des frontières
geom = rdc.geometry()
coordinates = geom.getInfo()["coordinates"]
shape = []
get_coords(coordinates, shape)

Starting for this shape structure, is there a way to recognize the useful tiles ? Is there an already existing algorithm that I can reuse ?


texts: How to abort a layer event in Openlayers

texts:
I have a problem ( How to stop features to redraw unwillingly in OpenLayers ) that might be solved by aborting the beforefeatureremoved or beforefeaturesremoved event.

What I would like to to is to stop the layer from refreshing (removing and the adding again) features from a vector layer.
If I for example make a javscript error in those methods the features isn't removed or refreshed. But causing a javascripts error causes other issues so I hope that there is a controlled way of doing this in Openlayers

I have tried the

layer.events.on({
   'beforefeatureremoved':function(evt){
       OpenLayers.Event.stop(evt);
    }
});

But it doesnt seem to do anything, Im I

Br
Mike


texts: QGIS: how to identify a visible raster layer among large view (composition of rasters)?

texts:
I have QGIS project loaded with many raster files.

For a given area, (ex: under my mouse pointer), I am looking for way to identify which (raster)layer provides the visible data.

To start, I will already be happy with a simple solution using the GUI.

I am working with QGIS 2.4 on Ubuntu


texts: Using Arcpy.AddMessage() to list out all shapefiles that were converted in a folder

texts:
I currently have code that re-projects all my shapefiles in a folder based off an input shapefile that I select. I will attach the code below. The code runs fine, but I would like the message on the geoprocessor window to say "Projected feature1, feature2, feature3, etc"

I currently can only get the script to list:
"Projected feature1"
"Projected feature2"
"Projected feature3"

It's basically printing the message for each iteration of the loop but I need it to print each feature class selected.

Attached is my example script below.

import arcpy
folderLocation = arcpy.GetParameterAsText(0)
targetShapefile = arcpy.GetParameterAsText(1)

try:

    targetDescribe = arcpy.Describe(targetShapefile)
    targetSpatialRef= targetDescribe.spatialReference
    targetSpatialName = targetSpatialRef.Name
    arcpy.env.workspace = folderLocation
    featureList = arcpy.ListFeatureClasses()
    for fc in featureList:
        desc = arcpy.Describe(fc)
        spatialRef = desc.spatialReference
        fcName = spatialRef.Name

        if fcName != targetSpatialName:
            rootName = fc.replace(".shp", "")
            outfc = folderLocation + '\\' + rootName + "_Projected"
            outCS = targetSpatialRef
            arcpy.Project_management(fc, outfc, outCS)
            arcpy.AddMessage("Projected "+ fc)

except:
    print arcpy.GetMessages()


texts: How do I reread a catalog after applying a spatial filter in readlasCatalog, to represent the reduced set of las files?

texts:
I would like to normalize a set of LAS files produced by Digital Aerial Photogrammetry (DAP) against a set of LAS files acquired by Lidar. The idea is to create the DEM raster from Lidar once and to reuse it when processing chunks for the DAP points.
The code creates two catalogues: one for the DAP, and one for the Lidar. Before entering into the function with catalog_apply(), I would like to create a DEM from the Lidar using the extent of the DAP. My code does reduce the LAS files referenced by the LAScatalog, but when applying the grid_terrain() function the now empty tiles are still processed but result in empty tiles (i.e. the catalogue recognizes the input option of the filter but the catalogue itself remains the original size). If I could reindex the catalogue it would only process LAS files that I need - is this possible? Or should I run grid_terrain() within the catalog_apply() function separately for each chunk.
library(lidR)

DAP_LASCatalogue <- readLAScatalog("haspeaks/LASDAP/", select = "xyzc")

# eventually I would like to input extents directly from the DAP Catalogue
DAPExtent <- extent(DAP_LASCatalogue)
min_x <- xmin(DAPExtent)
max_x <- xmax(DAPExtent)
min_y <- ymin(DAPExtent)
max_y <- ymax(DAPExtent)
print(c(min_x, min_y, max_x, max_y))

LIDAR_LASCatalogue <- readLAScatalog("J:/2018/BlockE_2017/LAS",
                                      select = "xyzc",
                                      filter = "-keep_xy 539052.1 5050814.6 556597.7 5059891.5",
                                      progress = TRUE)
plot(DAP_LASCatalogue)
plot(LIDAR_LASCatalogue)

LIDAR_DTMGrid <- grid_terrain(LIDAR_LASCatalogue, 1, knnidw(), use_class = c(2L, 9L), Wdegenerated = TRUE)


texts:
I am a completely new in PostGIS.

I tried to solve a problem I had in QGIS (see Simplify multiple lines: create central axis?).

I just need to know how the correct syntax works to use the function ST_ApproximateMedialAxis(geometry geom);

I succesfully added the geometry to the pgAdmin4 Geometry Viewer with the query

select geom from public.bufferrail; - so that seems to be correct.

But how to include this geometry as argument to the function?

SELECT ST_ApproximateMedialAxis (select geom from public.bufferrail);
does not work. I tried several variations, with no success.


texts: How to refresh graduated symbology legend in layers panel using python in QGIS 3?

texts:
In QGIS 2.x, I used the following code to refresh a graduated symbology legend in the layers panel:

self.iface.legendInterface().refreshLayerSymbology(self.myLayer)

The new QGIS 3 API advises to use QgsLayerTreeNode instead, but the only way I found to get it working is the following code :

root = QgsProject.instance().layerTreeRoot()
model = self.iface.layerTreeView().model()
model.refreshLayerLegend(root.findLayer(self.myLayer.id()))

... which looks rather cumbersome. What am I missing ?


texts: Getting latLng from google maps api to store that in PostGIS?

texts:
I am using Postgres with postGIS extension.  I just installed everything a while ago.

I'm looking at the documentation (it is really difficult for beginners in my opinion - is there any other tutorial or something you recommend?) and it seems I have to create a point from my lat and long, and that in a geospatial database it is (x,y) so I have to flip my lat and long to long, lat.

So, I'm creating markers in my front end application, but I'm confused as to how I can store them in my database.

Here's an example of the marker data that my web Api (using node and a postgres client) is sending over to Postgres (I am able to insert it into my table, so I have some test data):

{
  "lat": 40.7128,
  "lng": 74.0059,
  "title": "Big Apple",
  "description": "New York City, NY, USA.  The best city in the world!",
  "createdBy": 33823,
  "type": 4
}

My goal with using postGIS is to be able to do something like this:

I am using my app and I have my current location.  Find stuff that's within 1 mile of me.

Is there any guidance you can provide to me on how to get started towards that?

From the documentation, I need to create a point.  But then it delves into using a type of coordinate system to use.  At this point I have no idea what anything means.  I am trying to put together a proof of concept for an idea I've had for a while, but this is the biggest roadblock I've hit.


texts: Apply Symbology from Layer in ModelBuilder?

texts:
I have a model in ModelBuilder that takes a shapefile, runs some tools on it, adds it to a geodatabase, applies symbology from a layer file and converts it to a KMZ.  It all works fine except for the "Apply Symbology from Layer" tool.

I created a smaller model for testing the problem sections which is shown below.  The original input feature class is a geodatabase.  When I run this model, it runs and says it completes successfully, however my symbology is not applied to the layer.

The Input Layer and Symbology Layer are both point files, so that should not be the problem.  When I manually import the Symbology Layer (outside of ModelBuilder), it displays fine, so I don't think it's an issue with the layer file.

I tried to revise my model based on suggestions from Mapperz and MollyS.  Here is my revised model and screenshots of the two Properties windows but I am still unable to get a KMZ with the symbology.  Any ideas what I am doing wrong?


texts: How to uncheck "Close tool dialog" option in ArcGIS

texts:
I have to re-run some of my tools repeatedly. In order to save time I check "Close this dialog" option:

It is a struggle to uncheck this option for tool that takes a couple of seconds to complete. To stop dialog from closing, I sometimes have to run another tool, that takes longer to complete. In this case I have enough time to uncheck an option while tool is still running.

Are you aware of some better way to do so?


texts: show some data on a map by OpenStreetMap from C# VS2013

texts:
I am very new to Open Street Map.
I have some data that need to be shown on a map by Open Street Map from C# VS2013. It can be WPF or windows forms.

The data is about some statistics results by location.

For example,

 location,  latitude, longitude, address, population, household, income of household, tax of household

The data was got from a database by running SQL query from C# code.
The data volume may have 1 million rows. I need to do some aggregation so that they can be displayed on the map clearly. For example, aggregate all houshold (located in zip code 10080) information in a circle displayed on the map.

Moreover, I need to provide the following functions the on map so that the users can do:

(1) when zoom in to a location, more detailed data (from zip code to street level) will be shown.

(2) when zoom out, only the aggregated data was shown. For example, from street level to zip code level.

(3) users can get detailed data by selecting any location on the map as long as the location has data to show. For example, if a location in zip code = 10080 was selected, all household information in the selected location can be shown.

(4) users can hold control and click any location so that multiple locations household information can be show.

I have checked OSM wiki website and API, but I still have no idea how to add the data to a OSM map and edit the map to show it clearly.


texts: Creating geodatabase feature class using geodatabase table schema with ArcPy?

texts:
Is there a way in ArcPy to create an empty polygon geodatabase feature class that will acquire the schema of a geodatabase table?

I know that it's possible in ArcCatalog (New -> Feature Class and then import all the fields from a table) but I wonder if it can be written in a Python script.

The CreateFeatureclass_management tool allows only Feature Layers as schema templates.


texts: Creating menu item and adding more than 1 plugin in QGIS?

texts:
Because if I execute these lines for each plugin, I created a menu item for each plugin.

I created many plugins and wish to place them in a new menu item (not traditionals)


texts: ArcGIS server javascript export feature to shapefile identified by unique id

texts:
I am working on an ASP.NET MVC 4 application which has a map (using ArcGIS server 10.2 and JavaScript) as integral part of it. In that application, user draws a polygon object on the map and puts necessary attributes. Now I want to get that feature in archived shapefile which was drawn by the user (based on GUID or OBJECTID which I already have).

To get the shapefile, I published a geo-processing service using "Featureclass to Featureclass" but it fails all the time:

function exportShape() {
    alert("job started");
    var gp = new Geoprocessor("http://myserver/arcgis/rest/services/ExportToShape/GPServer/ExportObjToShape");
    var params = { where_clause: "OBJECTID = 231" };
    gp.submitJob(params, gpJobComplete, gpJobStatus, gpJobFailed);
    alert("job submitted");
}
function gpJobStatus(jobInfo) {
    console.log(jobInfo.jobStatus);
}
function gpJobComplete(jobInfo) {
    alert("completed");
}
function gpJobFailed(jobInfo) {
    alert("failed");
}

I get messages "job submitted" and "completed" but I never get the message "Failed", however when I check the console then I see the message esriGpFailed


texts: Creating a mask with transparent pixels in Python

texts:
I have a GeoTIFF image with pixel intensities ranging between certain values. Let's say between -2 and 2 (the image is normalized). I have a certain threshold which interests me, let's sat 1. I would like to create a NEW image that makes all the pixels that where previously larger than 1 to now have the maximum value, 2, or just some constant value. And all other pixels to be transparent.

I'm very new to this so I hope my question is understandable.


texts: QGIS program to search qgs / qgz with specific data source layer

texts:
I need to search for all layers pointing to a specific datasource across all the project files (.qgs / .qgz) in a directory or drive.

It's necessary in a myriad of situations, e.g. when deciding to either directly edit or copy and then edit a datasource, I usually need to check and make sure the edit is not rendering related layers useless or buggy; if so, I'd need to edit a copy of the datasource instead.

Another situation is when I know I have an already set up layer and need to reuse it but can't remember exactly which of the several hundred projects has it, and narrowing down the search to the ones using that particular datasource would come very handy.

I already have a script that does this for .mxd files. I just don't know enough Python for QGIS to be able to tell if it's possible. Nevertheless, it absolutely should be, I see no reason to think otherwise.

Solutions in any other programming language would be welcome as well (provided it's a language with tutorials online, please).


texts:
I have a process in which I need to use coverages. I can create the coverages just fine using FeatureClassToCoverage_conversion within arcpy but trying to union 2 coverages together does not work in arcpy with the identified method "Union_arc". I suspect it has something to do with creating an ArcINFO workspace, but I do not know how to use this for _arc commands.

The error is "Object: Tool or environment  not found"
I have ArcINFO installed and can run these tools successfully in Modelbuilder.
I have attached a test script to get Union working.
Running 10.1 ArcGIS python 2.7.

import arcpy, os, sys
from arcpy import env
print 'Starting....'
env.OverWriteOutput = True
arcpy.SetProduct("ArcInfo")

env.workspace = "G:\\Projects\\P747\\3_Landbase\\LB1\\temp\\test2\\"
TEMP = "\\\\silver\\clients\\Projects\\P747\\3_Landbase\\LB1\\TEMP\\"

print 'Loaded Data....'
if arcpy.Exists(TEMP + "airworkspace3"):
    arcpy.Delete_management(TEMP + "airworkspace3")
arcpy.CreateArcInfoWorkspace_management(TEMP, "airworkspace3")

r = str(1)
for fc in arcpy.ListFeatureClasses():
    print 'Working on ' + fc

    CovName = TEMP + "airworkspace3\\" + "Layer1"
    CovName3 = TEMP + "airworkspace3\\" + "Success"

    if arcpy.Exists(CovName + str(r)):
        arcpy.Delete_management(CovName + str(r))

    arcpy.FeatureclassToCoverage_conversion(fc + " POLYGON", CovName + r, "", "DOUBLE")
    r += str(1)

print 'Working on Union'
arcpy.Union_arc(CovName + "1", CovName2+ "2", CovName3)

print 'Completed ' + fc + ' processing....'


texts: SAGA Upslope Area Calculation from Vector Point Layer resolves in no output and no error message

texts:
I'm using the SAGA Upslope Area algorithm within pyqgis, but i'm stranded, since my script doesn't create any output at all. I'm trying to use geometry from points in a vector point layer as input x and y in the saga algorithm, but it doesn't give me any error messages and i can't seem to find the problem. Can you see any problem with my script below?

##points=vector
##dtm=raster
##upslope_area=output raster

from qgis.core import *
from PyQt4.QtCore import *
from processing import *

upslope_area = processing.getObject(upslope_area)
iterator_points = processing.getObject(points)
dtm = processing.getObject(dtm)

for feature in iterator_points.getFeatures():
    geometry = feature.geometry()
    point = geometry.asPoint()
    x_point = point[0]
    y_point = point[1]
    Processing.initialize()
    processing.runalg('saga:upslopearea', None,x_point,y_point,dtm,None,0,1.1,upslope_area)


texts: ArcGIS Report designer Attachment placement

texts:
I would like to know if it is possible for me to manually place each attachment where I would like within report designer?

The issue that I'm having when using the _ATTACHREL within the related report window is that when the report is run the attachments are showing up one after another.

I would like to be able to place them horizontal instead of vertical.


texts: Find traveling distance between line and point

texts:
Is there a clever way to traverse a line string up until a point and return the distance to that point? As an example:

a = LineString([(1, 1), (9, 1), (9, 6), (4, 6), (4, 3)])
p = Point((4,4,))

The distance from the start of a, to point p would be:

  8  # 9-1
  5  # 6-1
  5  # 9-4
+ 2  # 2, p intesects line segment after 2 units
= 20

I suppose I can just traverse the line and summing up the distance until p intersects the line segment, and calculate the last bit, but I get a feeling that there is something more elegant in the API that I haven't yet discovered.

Never heard of shapely until today, and I'm quite impressed by what I have seen so far.


texts: GPS tracking data string format for GRTQ device

texts:
I get GPS data from tcp (php socket), but for some devices (GRTQ) I will get an unreadable string.

First devices sends:

*HQ,XXXXXX,41,4#V1,time,A,**Lat**,N/S,**Lng**,W/E,000.00,000,date,FFFFFBFF,432,35,32448,334

then sends:

*HQ,XXXXXX,V9,time,A,**Lat**,N/S,**Lng**,W/E,000.10,000,date,,Caller ID,FFFFFBFF#

then sends:

$A  grQ05Ah@‘)���ÿÿûÿÿ����°#~À‚¡U

What's the meaning of this last string?

The model of the device is GRTQ.


texts: Filtering out Annotation Feature Classes using arcpy.da.Walk()?

texts:
Is there a way to use:

for dirpath, dirnames, filenames in arcpy.da.Walk():

to filter out annotation feature classes?

The help doesn't specify them as a type, which is strange because using the Feature Class Properties returned by the Describe function, you can filter out annotation.

I'm hoping to walk through a geodatabase and grab everything but the annotation FC's.


texts: Duplicate display of annotation in SLD script

texts:
I am using one SLD file in GeoServer of village parcel display.  It works fine, but the issue is that the annotation(parcelno) display duplicates/replicates when I zoom in.

What can I do?


texts: Summing up values of neighbouring polygons using QGIS?

texts:
I hope that you can help me with the following problem: I have a vector (polygon) layer. I would like to add an attribute to the layer that - for every polygon - sums the values of a specific field of all its neighboring polygons.

To give a more concrete example: I have a polygon layer of districts containing information on population. Now, for every district I would like to know how many people live in all its neighboring districts.

As I have more than 300 districts I cannot do this by hand for each and every district.

Is there any way do this more efficiently in QGIS?


texts: Can we access the ArcMap Image Analysis tools using Python?

texts:
Esri recently added a new Image Analysis window to ArcMap - it's found under > Windows in the main toolbar:

  The Image Analysis window's Processing section simplifies the
  experience of applying complex processing and analysis techniques to
  layers of image and raster data in ArcMap. The tools provide one-click
  options to apply processes, such as clipping, normalized difference
  vegetation index (NDVI) creation, mosaicking, and exporting

Can we access these tools using Python?

Specifically, I want to automate the process of applying a Clip, then adding Functions to the clip, before exporting to a new layer.


texts: Export raster from geodatabase to tif failure in ArcMap?

texts:
I have a raster in a file geodatabase and I am trying to export it out of the gdb to tif format.  Every time I try it Arc produces a tif file that no other programs seem to be able to read.

I've tried doing it with the Copy to Raster GP tool, right click on raster in TOC -> export data, and the Raster To Other Format tool.  The tools run successfully per Arc, and ArcMap can read the resulting raster.  But when I try to load in into AutoCAD Map3D 2017 (via MAPIINSERT or the Map Workspace connection manager) I get the message:

  (file path) was not found or is not valid

Photoshop CS5.1 just throws an error window saying:

  could not complete request because it is not right kind of document

File path is less than 132 characters & pixel depth/type is FLOAT 32 (it doesn't seem to mater if I specify FLOAT 32 or not when I run the export tools)


texts:
Is there a way to append only newly added records from a database to another database?  I looked online and couldn't find an answer.  I figured this would be somewhat straight forward.    We have a working database that I want only append new records to our 'live' database.  Both of them are File GDB.

Update 1
I have two File GDBs.  One is a working GDB where users can make edits.  The other GDB serves as a 'live' GDB which needs to get updated when the 'working' has changes made.  Right now I have to stop the service remove the live database and copy over the working database in it's place, then starting the service again for every client.  I was looking at a way to just copy the updated records so I don't have to stop the service and know that only new updates went to the live database.  Should I put the live database on a SDE?


texts: ArcGIS API for Python not plotting Spatially Enabled Dataframe

texts:
I'm using Python 3.6.8, Pandas 0.25.3, and ArcGIS API for Python 1.6.2. I'm pulling data from an API and viewing it in a pandas dataframe using a Jupyter Notebook. The data has several attributes including latitude and longitude. I can successfully convert the Pandas dataframe to a spatially enabled dataframe, call the .head() method and see the new SHAPE column. I can create a map object but when I try to plot my dataframe on the map object, the map does not update and display the points. I've included my code below:

# Define parameters for API call
headers = {
    'accept': 'application/json',
    'Api-Token': 'Key',
    'Agency-ID': 'ID'
}

# Call API and save to a varialbe
response = requests.get('URL', headers=headers).json()

# Convert json to dataframe
df = pd.DataFrame(response)

# Convert dataframe to spatially enabled dataframe
df = pd.DataFrame.spatial.from_xy(df, x_column = 'longitude', y_column = 'latitude')

# Create map object and display map
gis = GIS()
m = GIS().map()
m.center = [36.081179, -79.556955]
m.zoom = 9
m

# Plot dataframe
df.spatial.plot(map_widget= m)

Line 22 returns a map and Line 25 returns 'True' but no updates to the map are made. I've also tried using the renderer and symbol_type and symbol_style attributes.

df.spatial.plot(map_widget=m,
        symbol_type='simple',
        symbol_style='d', # d - for diamonds
        colors='Reds_r',
        cstep=10,
        outline_color='Blues',
        marker_size=10)

This block of code replaced Line 25 in the first section but did not display in the map either.

I've been through the documentation and the guide on GitHub but cannot find anything to make it work.


texts:
How would I create points along a line on multiple lines?

We are mapping our vineyards, and in some cases the rows are straight, and parallel. So for example, I've got 86 rows, and each row has 93 vines. How do I first create the 86 rows without having to copy parallel every time, and then how do I create the 93 vines along each one of those rows without needing to construct points every time?

I've been messing around, and can't figure it out.


texts: How to do nearest neighbor search more quickly?

texts:
I am currently using geopandas and scipy.spatial to do a nearest neighbor query. I have a data set of close to a million points and I need to find all of the closest points from another data set that has about 2,000 points. Right now I am making my tree like so

self.shop_tree = cKDTree(self.tree[['LONGITUDE', 'LATITUDE']].values)

and my query like so

dist, ind = self.shop_tree.query(self.pnts[['longitude', 'latitude']])

After I complete the query I need to also know the distance between each point, so would it be best to do this as a separate operation after, or is there a way to compute these simultaneously?


texts: When labelling polygons, not all labels are displayed - QGIS 1.7.0

texts:
I am trying to create a simple map of Mexico showing the state polygons and their label.

The state shapefile is from here, and the specific file is the first one from the top (Áreas Geoestadísticas Estatales (7.43 Mb)).

I edited the attribute table to create a new column with shorter state names called NOM_SHORT.

To add labels I used the "ABC" dialog as below:

However not all states were labelled, including Oaxaca and San Luis Potosí, as shown below (painted in light grey).  Why?

Here is attribute table:


texts: GEOS_ERROR: TopologyException: Input geom 1 is invalid: Self-intersection at or near point

texts:
I'm using GeoDjango (but I guess this is from the geos library).
I have a set of polygons that compose a city surface and I want to merge them into a single one.

The code looks like this:

polygons = [zone.poly for zone in self.zone_set.all()]
poly = polygons[0]
for pol in polygons[1:]:
    poly = poly.union(pol)
self.poly = poly
self.save()

The problem with this is I get a TopologyException after the first three mergings, stating there is a self-intersection

Shouldn't overlapping points be merged as a 'union'?


texts:
I've been trying to install the Orfeo Toolbox but this error keeps appearing:
RuntimeError: wrapped C/C++ object of type AttachRasterValuesToLineVertices has been deleted
Traceback (most recent call last):
  File "C:/Users/Nikoleta/AppData/Roaming/QGIS/QGIS3\profiles\name/python/plugins\TlugProcessing\ThToolBox_algorithms_provider.py", line 101, in loadAlgorithms
    self.addAlgorithm( alg )
RuntimeError: wrapped C/C++ object of type AttachRasterValuesToLineVertices has been deleted

Does anyone know how it is solved?


texts: Select By Attribute - Get Parameter as Text - User input SQL ArcPy

texts:
I am trying to create a script that will search for suitable site locations.

I want to let the user input an area size (e.g. >20) to be searched but each time I run the script it comes up as

  Error 000358: Invalid expression. Failed to execute
  (SelectLayerByAttribute).

This section of the code looks like:

inputvalue =arcpy.GetParameterAsText(0)
query = "POLY_AREA = " + inputvalue
arcpy.SelectLayerByAttribute_management(featurelayer2, "NEW_SELECTION", query)


texts: New Snapping while editing vertex's in ArcEngine 10

texts:
Is it possible to use the new (not classic) snapping functionality when editing the vertices of a polygon feature in ArcEngine 10. I would like editing vertices to behave exactly the same as editing vertices in ArcMap.

I know that with ArcEngine you pretty much have to control everything. So when we are creating the polygon we use the following function to capture the new snapped point using the snapping feedback object.

Is there some events I can tap into that'll allow me to tap into the editing vertices and update the mouse pointer to a snapped location?

    /// <summary>
    /// Updates a point to a snapped point
    /// </summary>
    /// <param name="currentPoint"></param>
    /// <returns></returns>
    private IPoint UpdateSnappedPoint(IPoint currentPoint)
    {
        ISnappingResult snapResult = null;
        IPoint point = currentPoint;

        if (_snapEnabled && _snappingFeedBack != null)
        {
            snapResult = _snappingEnvironment.PointSnapper.Snap(point);
            _snappingFeedBack.Update(snapResult, 0);
        }

        if (snapResult != null)
            point = snapResult.Location;

        return point;
    }

Niroshan


texts: I have a raster file that extends longitudinally from 0° to 360°. Using R, I'd like to wrap the raster at 180°, so that longitudes between 180° and 360° are translated to the interval [-180°, 0°]. Is there an existing function or package that can do this?


texts: Transfering content from Portal to ArcGIS Online?

texts:
I am assisting a coworker with developing a Story Map for his organization.  He created all of his content, but mistakenly created and saved everything within his organizational Portal.  Although his content is ready to be shared, it needs to be in ArcGIS Online due to the URL format.  Our web content system does not support the URL provided by Portal.

Is there a way to transfer content from a Portal account to an ArcGIS Online account without using a python script or having to recreate the Story Map and associated content?


texts: Point coordinate shift after python toolbox processing

texts:
I am working on custom clustering in ArcMap. For it I have created a python toolbox, where I do the processing. However, as a result I get slightly shifted points (6-20 meters), while they shall stay the same and just be grouped:

    point = arcpy.Point()

    # Spatial reference definition
    desc = arcpy.Describe(nodes)
    spatial_ref = desc.spatialReference

    for key in clustering.keys():
        cluster = 'CMPMT_Cluster_' + str(key)
        out_cluster = os.path.join(output_dir, 'DataTest.gdb', cluster)
        if len(clustering[key].keys()) > 1:
            point_list = []
            for k in clustering[key].keys():
                if k != 'tot_cost':
                    point.X = clustering[key][k]['XY'][0]
                    point.Y = clustering[key][k]['XY'][1]
                    point_geometry = arcpy.PointGeometry(point)
                    point_list.append(point_geometry)

            arcpy.CopyFeatures_management(point_list, out_cluster)
            arcpy.DefineProjection_management(out_cluster, spatial_ref)

I have checked the coordinates that are passed and they are fine, also I save another point file with the same spatial reference and it is okay. What could be the problem?

UPD: as Vince in comments suggested, I have examined how I pass the spatial reference to the geometry object. With the current updated code it works as needed:

point_geometry = arcpy.PointGeometry(point, spatial_ref)
point_list.append(point_geometry)
arcpy.CopyFeatures_management(point_list, out_cluster)


texts: Language Translation: map outputs to XLIFF or similar for the purpose of translation

texts:
A series of maps that i'll be producing with English text and labels will need to be sent to a certified translator. I'll be exporting the maps from GIS as PDF format.
When the output PDFs were sent to the translators for the previous phase of the project, the cost of translation was high because the text on map PDFs could not be read into their language translation tools (i believe the type of software is called Computer-Aided Translation, or CAT), so the process was very labour-intensive, and was a slow process. With this in mind, the project team have identified how we provide geospatial items to the translation consultant as something to be improved on.
I've been told that if i could provide the map figures in XLIFF format or similar, then the text elements in the file can be picked up by the translation tool. My understanding is that an XLIFF is a specifically formatted .xml for the purpose of CAT. Given that maps do have text elements on them but are graphical outputs by their nature, i'm not sure how this would work.
Is there a way to convert PDF maps to XLIFF or another format that plays nice with CAT? Better yet, can a map be exported in a suitable format directly from GIS? I suspect that we will need to settle with a workaround and then reintegrate the translated textual elements within GIS.


texts:
I'm working with QGIS 3.8 and I'm facing an error that I don't know how to fix.

I have download shapefile of all the countries in the world. When I display it on QGIS it has no problem at all:

I wanted to edit the table of this shapefile, so I went to the attribute table, used ctrl a and ctrl c in order to copy the table and paste it in Excel.
Then I have removed columns that I don't need (still have the WKT geometry, hasn't touched it), saved it as CSV and then I tried to load this layer to QGIS but I don't get all the countries, some have dissapeared even though I haven't touched them or edit them in Excel:

I don't understand what can lead to this kind of error and why, I have uploaded many CSV geometry files after editing them, so what can cause such an error?

My end goal is to be able to edit my layer in Excel and then upload as CSV without losing polygons.


texts: Minimum width of an irregular polygon (pyqgis)

texts: I would like to determinate minimum width of an irregular polygon (as those in red in picture) with pyqgis. Does anyone have any advice?


texts: PyQGIS: How to store the result of cliprasterbymasklayer

texts:
This Question is a follow up Question to this one I recently asked.

In my plugin I have a function that clips rasters by mask and I would like that the function returns the raster layer object. And then in another function I can get the Object. That's the code so far:

def clipRaster(self):
    wasserstand = self.getRasterLayer()
    outFile = os.path.join(QgsProject.instance().homePath(), "rtemp.tiff")
    gebeude = self.getVectorLayer()
    clip = processing.run('gdal:cliprasterbymasklayer',{'INPUT': wasserstand,'MASK': gebeude,
                   'NODATA': -1, 'ALPHA_BAND': False, 'CROP_TO_CUTLINE': True, 'KEEP_RESOLUTION': True,
                   'OPTIONS': "", 'DATA_TYPE': 5, 'OUTPUT': outFile})

    rlayer = QgsRasterLayer(clip['OUTPUT'])

    #rlayer = QgsRasterLayer(path_to_tif, "SRTM layer name")

    return rlayer

def getclip(self):
    clip = self.clipRaster()

    self.iface.addRasterLayer(clip)

    #iface.addRasterLayer("/path/to/raster/file.tif", "layer name you like")

Unfortunately, the error message I get is the following:

  TypeError QgisInterface.addRasterLayer() arguments did not match any overloaded call   overload 1 argument 1 has unexpected type 'QgsRasterLayer'   overload 2 argument 1 has unexpected type 'QgsRasterLayer'

Why isn't it working? "outFile" defines the path and name of the raster, the processing algorithm creates a raster layer and puts it into the "outFile" directory. To create a new raster the function QgsRasterLayer() needs a path and layer name which is both in outFile, thus the raster is now saved in the rlayer variable. Finally, iface.addRasterLayer(clip) adds this object. It makes sense doesn't it? :).

Can someone explain my mistake?


texts: "ERROR 000840 | The value is not a Relationship Class"

texts:
I'm stumped here. I've seen this error in other posts but not quite in the same context. I'm trying to run the arcpy.management.CreateRelationshipClass() tool in ArcGIS 10.0. Absolutely everything I've tried throws the 000840 error. I ultimately want to run the tool from a Python script, but for now I'm trying in the Python window in ArcMap. Here's my code:

import arcpy
from arcpy import env
env.workspace = 'C:/mulleady'
arcpy.management.CreateRelationshipClass('temp/SOIL_HYDROGROUP.dbf',
    'temp/WUH_Soil_Statistics.dbf', 'temp/relclass', 'SIMPLE', 'forward',
    'backward', 'NONE', 'ONE_TO_MANY', 'NONE', 'IL_SOIL', 'IL_SOIL')

I tried \\ instead of / in the pathnames, I tried to name the directory temp.gdb instead of temp, I ran arcpy.ValidateFieldName(), I tried absolute paths, I've tried a bunch of stuff that hasn't worked.

An even stranger thing to me is, is that I'm getting the same error even when I try to run the tool from ArcToolbox > Create Relationship Class:

Any thoughts?

Edit

I created a file geodatabase and copied my two .dbf files into it. Now when I try to create a relationship class I get this error. I changed my workspace and scratch workspace to the path of the fgbd folder. What am I doing wrong now?


texts: Longest flow path in QGIS from single feature vector

texts:
I have got a vector layer with one feature (MultiLineString) representing a river network (s. below) without bifurcations. In addition there is a point shapefile containing the outlets of different subcatchments and one polygon for each of them.

Example: in the following image, the length should be from A (outlet, point in shapefile) to B (should be found automatically). As the point at A is not exactly on the line, an ideal solution would replace the outlet onto the line.

Is there a simple way to calculate the longest flow path (one single branch) in QGIS? The final goal is to have one length for each subcatchment.


texts:
Which algorithm does QGIS use for clipping shapefiles?
When I use ogr2ogr to clip a large shapefile (approx. 3gb) it takes 77sec. When using QGIS it takes approx 5 sec, clipping the same dataset (and same mask-layer). New clipped layer is approx 30mb, and the 6sec. includes saving in QGIS
Any suggestions for a better algorithms or modifications of the code below?
Here is the code I use in Python:
import subprocess
clipping_shp = r"C:\Python\buffer.shp"
output_shp = r"C:\Python\result.shp"
input_shp = r'C:\Python\data.shp'
subprocess.call(["ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", clipping_shp, output_shp, input_shp], shell=True)


texts: Building query for query attribute task in arcgis server manager for a layer related to a table

texts:
I my map document (.mxd), one layer (BUILDING with spatial data) is RELATED to a stand-alone table (EMPLOYEE) with one-to-many relationship having a common field call dept number (one building with a dept number can have many employees with same dept number in EMPLOYEE table). How do I build a query for Query Attribute Task through Arcgis Server Manager because when I choose BUILDING as my layer, it doesn't show EMPLOYEE table attribute. What I intended to write a equivalent sql query as "select * from building, employee where building.dept_no = employee.dept_no". I don't see a way to put this in arcgis server manager. Can anybody please help.

Thanks,
Archana


texts:
I am attempting to use a bash script to reproject wrf output using the following:
${gdal} -of GTiff -a_srs '+proj=longlat +datum=WGS84 +ellps=WGS84 +lat_1=30 +lat_2=60 +lat_0=42.1200065612793 +lon_0=-72.05000305175781 +x_0=0 +y_0=0 +units=degrees +no_defs' -a_ullr -74.98987 45.74738 -69.98447 40.22901 NETCDF:${data}"wrfout_d03_${loopdate}00.nc":${array[$index]} ${data}${modeldate}/wrf_${modeldate}_${1}_${array3[$index]}.tif &
While this code worked perfectly well with previous versions of PROJ, I am using version 6.3.0 and am met with the following error:
ERROR 1: PROJ: proj_create: Error -7: unknown unit conversion id
ERROR 1: Failed to process SRS definition: +proj=longlat +datum=WGS84 +ellps=WGS84 +lat_1=30 +lat_2=60 +lat_0=42.1200065612793 +lon_0=-72.05000305175781 +x_0=0 +y_0=0 +units=degrees +no_defs
Has anyone else run into this issue?


texts: Action: Move pointfeature to other layer - Point becomes Multipoint

texts:
I'm trying to write an action that moves a pointfeature from one layer to another.
It works fine until i try to save the layer with the new feature on. He says he can't save a MULTIPOINT to a point-layer.

Any idea why it suddenly becomes a Multipoint?

I based this code on the following question: How to create a QGIS action which copies features between layers?

Here is my code:

def getVectorLayerByName( layerName ):
  layerMap = QgsMapLayerRegistry.instance().mapLayers()
  for name, layer in layerMap.iteritems():
    if layer.type() == QgsMapLayer.VectorLayer and layer.name() == layerName:
      if layer.isValid():
        return layer
      else:
        return None

iface = qgis.utils.iface
sourcelayer = getVectorLayerByName('plane')
iface.setActiveLayer(sourcelayer)
sourcelayer.startEditing()
iface.actionCutFeatures().trigger()
sourcelayer.commitChanges()
targetlayer = getVectorLayerByName('dots')
iface.setActiveLayer(targetlayer)
targetlayer.startEditing()
iface.actionPasteFeatures().trigger()
targetlayer.commitChanges()


texts:
Can anyone help me what's the algorythm to calculate points of intersection whet is given one point p1 [lat,long] , bearing x (let say 137 degrees), arc center [lat,long] and radius?
I woulde like to get p2[lat,long] and p3[lat,long].

For better imagination what I mean something like on this picture but on sphere with only one point and bearing:


texts: Applying Snap_edit to polygon Geometry Object

texts:
I have a project where I am modifying over 1 million polygons (representing tracts of land) to make them line up with the PLSS land grid. My process works (mostly) by cutting up polygons based on a grid of the landgrid, determining whether to toss the split part or apply the grid polygon, the reassembles the shape to be inserted into the final dataset. Unfortunately this process over-modifies many irregular polygons, which I can identify by determining the result's percentage change in acreage, perimeter, and movement of the centroid.

In these instances (of over-modification) I would like to attempt to simply snap the original input polygon to the landgrid, so that at least the irregular polygons line up with their modified neighbors. Due to the nature of the process, I would have to apply this snap to the geometry object itself, not the entire dataset. The result of the snap will also be checked to determine whether the modification is too extreme.

My question: Can arcpy.Snap_edit be applied to a geometry object? Is there an alternative tool/process that could snap an input polygon to a line dataset?

I have tried using snap on the polygon shape object, with 0 results. I have also tried converting the vertices to point geometries (in the hopes that I could snap the points to the grid, then rebuild the polygon based on the points) but the resulting coordinates are the same.

Here's my attempt at getting this to work.

with arcpy.da.UpdateCursor(Polygons, [u'SHAPE@']) as cursor:
for row in cursor:
    if row[0].getArea('GEODESIC', 'ACRES') > 35.0:
        count +=1
        points = json.loads(row[0].JSON)["rings"][0]
        #print json.loads(row[0].JSON)["rings"][0]
        newPoints = arcpy.Array()
        for pt in points:
            print pt
            thePoint = arcpy.PointGeometry(arcpy.Point(pt[0], pt[1]), arcpy.SpatialReference(5070))
            print thePoint.JSON
            arcpy.Snap_edit(thePoint, [[grid, "VERTEX", "100 METERS"], [grid, 'EDGE', '100 METERS']])
            snappedXY = [json.loads(thePoint.JSON)["x"], json.loads(thePoint.JSON)["y"]]
            print snappedXY
            newPoints.append(arcpy.Point(snappedXY[0], snappedXY[1]))
        shape = row[0]
        NewShape = arcpy.Polygon(newPoints)
        #print NewShape[0].WKT == shape.WKT
        snap_env = [[grid, "VERTEX", "60 METERS"], [grid, 'EDGE', '50 METERS']]
        geom = arcpy.CopyFeatures_management(arcpy.Snap_edit(NewShape, snap_env), arcpy.Geometry())
        row[0] = NewShape
        cursor.updateRow(row)
    else:
        break
del cursor

UPDATE: Here's the working function:

def ShapeSnap(shape, grid):
    #Snaps input shape to grid, removes extra boundary lines
    if shape.getArea('GEODESIC', 'ACRES') > 29.0:

        points = json.loads(shape.JSON)["rings"][0]
        #print json.loads(row[0].JSON)["rings"][0]
        newPoints = arcpy.Array()
        for pt in points:
            pointFeature = "in_memory" + "//" + "point"
            if arcpy.Exists(pointFeature): arcpy.Delete_management(pointFeature)
            thePoint = arcpy.PointGeometry(arcpy.Point(pt[0], pt[1]), arcpy.SpatialReference(5070))
            arcpy.CopyFeatures_management(thePoint, pointFeature)
            arcpy.Snap_edit(pointFeature, [[grid, "VERTEX", "20 METERS"], [grid, 'EDGE', '20 METERS']])
            snappedXY = [json.loads(thePoint.JSON)["x"], json.loads(thePoint.JSON)["y"]]
            newPoints.append(arcpy.Point(snappedXY[0], snappedXY[1]))
        shapeobj = shape
        NewShape = arcpy.Polygon(newPoints)
        #print NewShape[0].WKT == shape.WKT
        newPolyFeature = 'in_memory' + '//' + 'newPolygon'
        arcpy.CopyFeatures_management(NewShape, newPolyFeature)
        snap_env = [[grid, "VERTEX", "40 METERS"], [grid, 'EDGE', '40 METERS']]
        arcpy.Snap_edit(newPolyFeature, snap_env)
        geom = arcpy.CopyFeatures_management(newPolyFeature, arcpy.Geometry())
        return geom[0]

    else:
        return shape


texts: QGIS: new install, but 8 icons to choose from. Which to use?

texts:
A very naive and basic question.  I have just downloaded QGIS 2.14.0 and now have 8 new & related icons on my desktop.  As a student with no practical GIS work experience, I have worked with ArcGIS 10, and now for it's proclaimed benefits would like to learn to drive QGIS, but am a little confused over which of the choices is appropriate.  Also, is one of these icons a similar application to ENVI?

I am using Windows 8.1.


texts: Get different SRIDs of geometry column in a table

texts: It is possible that a geometry column in a table of PostGIS has different SRIDs. I've used Find_SRID to get the SRID of this column but seems it only returns one value for SRID which is the first row. What is the reason?


texts: Rasterization of polygons, calculation of the area covered

texts:
In R environment I have (1) a SpatialPolygonDataFrame with number of polygons: polys; (2) raster grid of resolution 10'000 m: grid. Both of them are in the LAEA coordination system.

What I need is to get the area covered by polygons in each individual grid cell. I tried the following function:

r_area_covered <- rasterize(polys, grid, getCover = T)

It basically works and the result is roughly ok. However, in more detailed view it is quite inconsistent. Some polygons were ignored while other polygons of smaller sizes were considered for the calculation.

My questions are: where can be the mistake and how to prevent it? Do I have other options to do this in R.


texts:
I have created a arctoolbox with scripts in it. In the script, I set the source script file (a python file), my question is when I run this python using a batch file, it won't cause error, unfortunately when I use the arctoolbox script to run the python script, I got the following error:
Tool or environment  not found, Failed to execute (AttributeSearch).

The following is my python script details:

'''----------------------------------------------------------------------------------
 Tool Name:     SelectLayerAttribute
 Source Name:   SelectLayerAttribute.py
 Version:       ArcGIS 10.2
----------------------------------------------------------------------------------'''
import arcpy
# Set the workspace
arcpy.env.workspace = "c:/Users/lyao/Desktop/PUB_GIS_Tagging/GIS_DMS.gdb"
# get user supplied path, layers and fields
Value = arcpy.GetParameterAsText(3)
# Make a layer from the feature class
feature = "c:/Users/lyao/Desktop/PUB_GIS_Tagging/GIS_DMS.gdb/BQ_Segment_test"
arcpy.MakeFeatureLayer_management(feature, 'lyr')
where = '"DrainName" = ' + "'%s'" %Value
#arcpy.SelectLayerByAttribute_management("lyr", "NEW_SELECTION"," \"ROAD_NAME\" = "'%s'" % Value ")
arcpy.SelectLayerByAttribute_management("lyr", "NEW_SELECTION",where)


texts: Frequency analysis of linestring using QGIS

texts:
I have a data-set that consists of linestrings across different timestamps. I want to identify the segments that occur the most frequently.

I am using R to import the data from JSON files, create spatial objects to then import into QGIS. I would prefer to do the spatial analysis using QGIS, but I am facing a few hurdles irrespective of the approach I try. These are the approaches I am trying:

Use another layer that contains polygons as the baseline to identify the polygons that have the most number of intersections with the linestrings. [this will help me identify polygons, but once that has been done, what will be the optimal way to identify the linestring that is most frequent within these polygon(s)]
Use another layer of points/lines as a base layer to identify the ones that have the linestrings close to them. [But again, I am stumbling at identifying the actual linestring]

It feels like I may able to identify zones/points in the base layer but how do I find the linestrings (or section of linestrings) that occur most frequently using QGIS?


texts: How to implement a spatial query on polygon data for a specific point?

texts:
I have vector data fed into a database which returns a report at province level upon querying. This is implemented in a web tool.

I want to add an extra functionality at point level, where I can input lat/lon and run a spatial query to calculate some attribute, say distance to coast, based on that point.
Let us assume I have data which contains the distance from centroid of the province.

Please suggest me an open source approach to implement this. How feasible would it be to divide the province into grids and dynamically take up the values?


texts: Improving QGIS Server postgres layers performance?

texts:
I'm currently exploring QGIS-server 3.2.3 and QWC2.

Currently, I'm experiencing some issues with my postgres based vector layers, so I seek hints or pointers to my questions:

PostgreSQL Connection pooling: Accessing the server via one QWC2 session I got nearly 100 postgres session. How is connection pooling working within QGIS-server? How many open connections is QGIS-server using and is this configurable?
Tracing the database, I see that SELECT count(*) FROM "mytable" is issued very often. This might be an alive-check query for pooled connections(?) As this is a costly query on my large tables: where does it come from and is it configurable?
The inital GetCapabilities request takes very long, again because my tables contain millions of records and for every layer a SELECT ST_Extent("wkb_geometry") FROM "mytable" AS ogrpgextent is issued. I wonder, why st_extent is requested and not st_estimated_extent, as the table is analyzed and SELECT count(*) FROM pg_stats WHERE schemaname='public' AND tablename='mytable' AND attname='wkb_geometry'; returns 1. And where does "ogrpgextent" come from, as my layers are defined with providerKey 'postgres' and not 'ogr'? It seems to me, that the database is not accessed via qgspostgresprovider but ogr(?)
How can I enable verbose logging? I set FcgidInitialEnv QGIS_DEBUG 1, FcgidInitialEnv QGIS_SERVER_LOG_FILE /tmp/qgis-000.log,   FcgidInitialEnv QGIS_SERVER_LOG_LEVEL 0 but I don't see any log QgsDebugMsg-statements which I'd expect to see...
The docs state that QGIS has to be compiled with Postgres support. What does it take to compile QGIS with Postgres support? How can I figure out, if docker image qgis/qgis:final-3.2.3 is compiled with Postgres support? Docker container build-deps contains postgres-client , but unlike other databases the cmake call in the dockerfile  does not include a flag like WITH_POSTGRES which I'd expect (in analogy to e.g. WITH_QSPATIALITE)?

The layers in the .qgs file are defined like this:

<layer-tree-layer expanded="0" id="mytable_20180912153643064"
source="dbname='mypassword' host=192.168.33.10 port=5432 user='myuser'
sslmode=disable key='objectid' srid=31467 type=MultiLineString
table=&quot;public&quot;.&quot;mytable&quot; (wkb_geometry) sql="
checked="Qt::Checked" name="mytable" providerKey="postgres">...</layer-tree-layer>


texts:
I´ve a .gpkg file in Mollweide projection. When I open it in geopandas, it shows me  the correct crs information:
df1.crs
>>>>> {'init': 'epsg:54009'}

I´d like to re-project it to EPSG 4326, in a way to cross this layer with a global map in lat/long. If I try:
df1.to_crs(epsg=4326)

I get the error message:

pyproj.exceptions.CRSError: Invalid projection: +init=epsg:54009 +type=crs

How can I proceed to reproject those layers in a way I could intersect them?


texts:
I'd like to put maps inside of a wordpress CMS, but I'm running into troubles with the permalink structure ignoring the maps all together...

I have a .htaccess file that will force WordPress to ignore a directory, and I can have the page render, but not inside the pages...

I've tried embedding an iframe, both from an .htaccess allowed directory, and even from a subdomain that has nothing but a test map on it.

I'm fairly new to GIS, so any help in connecting a MapServer / MapScript system into WordPress would be genuinely appreciated.

If there's any other data you need, please let me know.

Matt


texts:
I have successfully opened a grib2 file from NCEP and I am having trouble being able to convert the coordinates to plot them using matplotlib, using the custom convertXY function from this post Plot GDAL raster using matplotlib Basemap.

I got what I expect, but only for half of the world, I can solve it by subtracting 180.0 from my xmin and xmax, but then I lose the coordinate conversion, I guess the problem is that I am not shifting the data, possibly using shiftgrid from mpl_toolkits, but I can not get the function to work either, any suggestions?

Here is an image of the map without the subtraction:

Here is what I get when I subtract 180.0 from the xmin and xmax variables:

You can download the grib2 file I am using from:
https://drive.google.com/open?id=1RsuiznRMbJNpNsrQeXEunvVsWZJ0tL2d

from mpl_toolkits.basemap import Basemap
import osr, gdal
import matplotlib.pyplot as plt
import numpy as np

def convertXY(xy_source, inproj, outproj):
# function to convert coordinates

    shape = xy_source[0,:,:].shape
    size = xy_source[0,:,:].size

    # the ct object takes and returns pairs of x,y, not 2d grids
    # so the the grid needs to be reshaped (flattened) and back.
    ct = osr.CoordinateTransformation(inproj, outproj)
    xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T))

    xx = xy_target[:,0].reshape(shape)
    yy = xy_target[:,1].reshape(shape)

    return xx, yy

# Read the data and metadata
ds = gdal.Open(r'D:\Downloads\flxf2018101912.01.2018101912.grb2')

data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()

xres = gt[1]
yres = gt[5]

# get the edge coordinates and add half the resolution
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmin -= 180.0
xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
xmax -= 180.0
ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
ymax = gt[3] - yres * 0.5

ds = None

# create a grid of xy coordinates in the original projection
xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]

# Create the figure and basemap object
fig = plt.figure(figsize=(12, 6))
m = Basemap(projection='robin', lon_0=0, resolution='c')

# Create the projection objects for the convertion
# original (Albers)
inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)

# Get the target projection from the basemap object
outproj = osr.SpatialReference()
outproj.ImportFromProj4(m.proj4string)

# Convert from source projection to basemap projection
xx, yy = convertXY(xy_source, inproj, outproj)

# plot the data (first layer)
im1 = m.pcolormesh(xx, yy, data[0,:,:].T, cmap=plt.cm.jet)

# annotate
m.drawcountries()
m.drawcoastlines(linewidth=.5)

plt.show()


texts: Changing color of each selected feature based on expression using PyQGIS?

texts:
Code below changes the color of all the selected features based on the last condition it satisfied instead of changing the color of the feature which satisfied the condition.

for f in FEATURES:

    total = f['parkeer_approx']
    alloted = parking_capacity[f.id()-int(1)]['parking_capacity']

    print(total,alloted)

    if total == alloted:
        print('in red')
        iface.mapCanvas().setSelectionColor(QColor("red"))

    elif alloted > 0 and alloted < total:
        print('in orange')
        iface.mapCanvas().setSelectionColor(QColor("yellow"))

    else:
        print('in green')
        iface.mapCanvas().setSelectionColor(QColor("green"))

What I want is to change the color of each selected feature separately based on the conditions defined in the code above.

How can I achieve that?


texts: Including labels in map legend and point symbology of ArcGIS for Desktop?

texts:
I've been using gis.stackexchange for some time already and time has come for me to ask my first question here.

I have point feature class with habitat types. There are two important attribute fields. One with the name of the habitat and second with a corresponding number.

I've placed labels inside the point symbology to represent the number on the map but I have trouble showing them in the legend. I need the legend to include both description and the number of the habitat.

Any suggestions?


texts: Interpolating within a boundary using QGIS

texts:
I am interpolating species density in a tidal flat however I would like to limit my interpolation to the sampled area in the tidal flat.

I tried making a shapefile of the sampled tidal flat area (barrier) and in Toolbox -> Thin plate spline -> output extent set to the polygon.

This did not produce the result I am looking for.
What is the correct process to limit the interpolation to my sampled area boundary?
I only have access to QGIS.


texts: Same python code works in arcgis python window, but fails in pyScripter

texts:
This code works in the arcgis python window but in PyScripter, this warning message had been written:

  Traceback (most recent call last): File "", line 16, in
   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy_init_.py", line 1167, in
  SearchCursor return gp.searchCursor(dataset, where_clause,
  spatial_reference, fields, sort_fields) File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing_base.py", line 359,
  in searchCursor self._gp.SearchCursor(*gp_fixargs(args, True)))
  IOError:

import arcpy
arcpy.env.workspace = r"C:\Desktop\spatial statistics analysis tutorials\PatternAnalysis\PatternAnalysisExercise\DengueData.gdb"

# place all the rows from the feature class into a search cursor
cursor = arcpy.SearchCursor("IndCases")

# iterate through the cursor and print the attributes to the interactive window
# print it in a standard three line format
for row in cur:
    print row.AGE

I was wondering why this situation happened, any suggestions?


texts: Clustering population data respecting given regional divisions?

texts:
I am familiar with clustering algorithms from machine learning and signal processing, for example for training vector quantisers. Recently, I have come across a problem in the GIS area that I am unfamiliar with. What if you want to cluster a land area according to population density. I.e., you have a fine-grained grid of population numbers across the area. You want to cluster this area into regions with approximately homogeneous density and you want clusters with fairly few inhabitants (but > 150) and you want the clusters to respect some pre-defined borders. Especially the latter constraint is important. Which approaches (if any) exist for doing this?

I have come across DBSCAN here: https://www.r-bloggers.com/geographic-clustering-of-uk-cities/. Is that a possible approach?

If possible it would be great if there is a Python implementation and/or an ArcGIS-compatible solution.


texts: Overlap custom TileLayer with satellite layer in Goolge MAP API

texts:
I've added a new mapTypes TileLayer (png tiles, transparent background) into my mapOptions that ables me to switch from others mapTypes to the new one, using button option. At the moment, when I change of mapTypes the new layer appears and the old layer disappears ==> I can't use 2 layers at the same time.

How can I overlap 2 layers into one mapTypes ? To understand the method I would like to overlap Satellite (from google) and the new one.

The code I use ables me to changes of layer, but it doesn't overlap: Note: I use an API KEY that I can't provide..


texts: Using raster data from vector and raster files in PyQGIS?

texts:
I am a complete beginner to GIS. I am using QGIS 2.6 32bit on Windows 7.

I have large set of raster data files, each file contains a 1km x 1km gridded scalar quantity X for a different year/month.  I have loaded one file into QGIS as a raster layer, set the correct coordinate reference system and set the raster layer to a blue banded rendering.

I have loaded three vector files, one is a 'country' outline, on contains 'county' boundaries while the last contains 'height contours' within the country. They all have the correct coordinate reference system set.

I can see that quantity X is varying in relation to the vector contour heights as I would expected. So far so good.

I want to obtain data relating X and average contour height from some predefined areas of the 'map'. I expect these to be a set of smaller areas along a line eg a section across a mountain.

For each raster data file I want to obtain an output data file that contains, for each of the 1km x 1km grid locations defined above, both the value of X and the average contour height in that grid.

From each output files I can then work out any correlations, variations over time and location etc by using a spreadsheet and a scatter graph.

For example the output file may contain 20 areas along a line in increasing height. Thus 20 sets of X and height.

I assume I will need to create a raster version of average heights from the contour data [but aligned with the X data 1km x 1km grids].

How do I do this by selecting a region or line as described ?

If this can be done, how could it automated in python ?


texts: Read an INSPIRE GML in QGIS2.0 : gml:xlink resolving

texts:
I'm trying to read an INSPIRE GML, which contains internal and external links between the elements. I saw on the GDAL site that we can set environment variable in order to load this link in QGIS (http://www.gdal.org/ogr/drv_gml.html) but it seems it doesn't work on this QGIS version.

The impacted lines look like this one:

Does someone have some experience in integration of this type of model?


texts: win10/CREATE EXTENSION Pgrouting--could not access file "$libdir/libpgrouting-2.3"

texts:
My system version is win10 10.0.14393.

I install the binaries-win64 into my C:/PostgreSQL/pg96.

At the end of the installing, there is a Warning "could not create database file "$libdir/libpgrouting-2.3": No such file or directory".

I tried to use psql to install using "CREATE EXTENSION Pgrouting", I got the same bad return: "could not access file "$libdir/libpgrouting-2.3": No such file or directory".

How do I deal with this problem.


texts: Any clue how to read a .style (ArcGIS) file in Python, for converting into another platform-specific format?

texts:
The problem I am face with is that I have all my data as shapefiles but my organization works with only 1 ArcGIS license. Migrating to QGIS to prominent amount of work, I need to convert my ArcGIS .style file to the .sld format of QGIS, to ensure uniformity of data capture and maps.

I am working with Python for this, Could anyone give me any suggestions?


texts: How to use gdal warp with modified crop_to_cutline?

texts:
I know that this question might be confusing, so I'll try my best to explain it.

I have a shapefile that I will be using to clip my raster:

This is my shapefile (whole island of Palawan, Philippines):

And this is my input raster to be clipped/extracted:

Now using gdalwarp, my result looks like this:

My question is, can I produce an output where its extent is based on the input raster only (excluding the upper part covered by shapefile) and does not include the whole shapefile?

This is the ouput raster extent that I want to have (as covered by green bounding box).


texts: Loop For to export polygons from a Feature Class

texts:
I need to export 150 polygons from a feature Class. I have a python snippet code to export the polygons one by one.
I think a for loop code can do it fast, but I don't know how to write the code into the loop. Can anyone tell me how to do it?

arcpy.FeatureClassToFeatureClass_conversion(in_features="D:/ASSUKKARGIS3/ty_plg_gglimgs_georef/plg_clip_imgs.shp", out_path="D:/ASSUKKARGIS3/plg_clips", out_name="001_clip_plg.shp", where_clause='"plg" = 1', field_mapping='plg "plg" true true false 19 Double 0 0 ,First,#,D:\ASSUKKARGIS3\ty_plg_gglimgs_georef\plg_clip_imgs.shp,plg,-1,-1', config_keyword="")


texts: Calculating area of MultiPolygons within 250m buffer of points

texts:
I am using using PostgreSQL 9.2 and PostGIS 2.0. I have two tables: 1) a geocoded table of home addresses (address_list) and 2) a table containing data on tree density in the same geographical area. I have included examples of the relevant columns in each table below.

address_list:

house_id (serial) |  geom(Point, 4283)
------------------+-----------------------
        1         |
        2         |

tree_density:

  gid (serial)  |  tree_dens |  geom(MultiPolygon, 4283)
    ------------+------------+----------------------------
          1     |  scattered |
          2     |  medium    |
          3     |  dense     |

What I would like to do is compute the area of forest (m2) by tree_density categories within 250m of each home address for my whole address table (approx 6,000 entries). (E.g for house_id = 1: area_scattered_within250m = 50m2,  area_med_within250m = 0m2, area_dense_within250m = 20m2)

Can this be done easily? Even if in 3 or more seperate queries? Does anyone have a moment to assist me with the syntax? I guess I need to combine ST_Area, ST_Buffer and perhaps ST_Contains or St_Within, but I'm new to PostGIS and although I've been ok with other (simpler) spatial joins, I'm lost here and haven't found any similar examples. For the area of my study I transform to SRID 3112 (meters) for measurements.

Many thanks
Mike


texts:
I'm using gdal_grid to convert a shapefile with points and I would like to set the pixel size of the output. I've only found the -outsize flag that lets me set the size of the output raster, but I would prefer to just say that I want 2 x 2 m pixels and let the tool figure out the size of the output raster.

With gdal_rasterize you can set the -tr flag, but that tool won't do interpolation.


texts:
I have a Raster representing climate data. I converted the raster to vector using the tool raster to vector (Arc Map) and i got something like below:

Can somebody help me out by guiding me in actual smoothing the edges to make it look like a polygon. Im using Arc Map 10.2 and also except answers for QGIS and ILWIS Software.


texts: Regarding spatial relationships, what is the Dimensionally Extended 9 Intersection Model (DE-9IM)?

texts:
In plain English what is the Dimensionally Extended nine-Intersection Model (DE-9IM)?

What is the meaning of each row/column in the matrix?

How can I use that information to build a relation matrix?


texts: Using PostGIS ST_Clip to crop raster in Django webapp

texts:
I'm currently building a geospatial webapp with Django and PostGIS on the back. Within the Django webapp, users can draw a polygon on a map, which is saved to a Django PolygonField() field. An example entry in this field may look something like:
SRID=4326;POLYGON ((0.9063720703125029 52.32023207609735, 0.8239746093749998 52.10819209746323, 1.170043945312496 52.14191683166823, 1.170043945312496 52.31351619974807, 0.9063720703125029 52.32023207609735))

On my server I have a reasonably large raster file (a .tif that's ~200Mb) of landcover data which covers the extent of the British Isles. What I would like to do is take the large raster file, clip/crop it to the dimensions of the user-created polygon and then save the cropped .tif file to a GeoDjango RasterField().
I was wondering what the best way of approaching this problem would be. For example, would it be best to handle the crop in a task/view within Django (e.g. with libraries such as rasterio or shapely) or somehow leave this for PostGIS to handle with a call to ST_Clip?
My feeling here would be to use ST_Clip but, given my large landcover raster isn't itself within a model, I don't think I could do something like:
myModel.objects.raw('SELECT ST_Clip(landcover_raster, 1, user_polygon)
    from rasters
    WHERE id = pk;')

How can I crop and save to a Django RasterField()?


texts: The spatial index grid size is invalid [feature set]

texts:
I am building a model in ArcMap ModelBuilder and I want to implement Feature Set as a model parameter so the user can select areas but I get this error:

  Create: An unexpected failure occurred. The spatial index grid size is
  invalid. [feature set] The spatial index grid size is invalid.

I even created another simple model that was just a Clip tool with the Feature Set as a model parameter but I got the same error when running it. The feature class I am using as an input is a polygon fc with all waterbodies in a county.

I don't really understand the error and what's causing it and Google isn't helping.

I am using ArcMap 10.4.1.


texts: QGIS Modeller - Are there system variables

texts:
In a QGIS 3.2 model can we use system variables?

Basically I want the input filename to be used in the output filename.

So

inputName.shp -->run process-->inputName_processed.shp

In arcGIS it's %inputName% and the value comes from Get Field Value - http://pro.arcgis.com/en/pro-app/tool-reference/modelbuilder-toolbox/get-field-value.htm

===update

based on the answer from Erik -how do i use layer_name in the final column 'calculated' in this? i want it be the input name with '_processed' at the end.


texts: Drawing circle of specified radius in meters using OpenLayers 3

texts:
I need to present a mask where the user put the radius (in meters) to draw a circle. How can I do this? Do I have to do geodesic calculations to get a circle that follows the projection of the map (EPSG:3857)?

new ol.geom.Circle([lon, lat], radius);


texts: Referencing Dataframe with ArcPy in Geoprocessing Service?

texts:
I have a geoprocessing service that was successfully published to our Server instance. I'm using 10.3 which means we could not publish unless it was successfully run as a tool in ArcMap. But now when I try calling the GP service form my client I am getting an error that the dataframe object was not created properly. I am not sure why the script works in ArcMap and not on the server.  Here is the script

import sys, string, arcpy, os

# get default script location
scriptPath = sys.path[0]

# Construct paths to mxd
mxdFilePath = os.path.join(scriptPath,"Mymxd.mxd")
mxd = arcpy.mapping.MapDocument(mxdFilePath)
df = arcpy.mapping.ListDataFrames(mxd)[0]
shpPath = arcpy.GetParameterAsText(0)
newLayer = arcpy.mapping.Layer(shpPath)
arcpy.mapping.AddLayer(df,newLayer,"TOP")
arcpy.RefreshActiveView()
del newLayer
del mxd

And here is the error

esriJobMessageTypeInformative: Submitted.
esriJobMessageTypeInformative: Executing...
esriJobMessageTypeInformative: Executing (MyService): MyService "Feature Set" C:\Users\~\AppData\Local\Temp\MyService\MyService_122215_1630_gpserver\j25ae6595a46840369acebfbd9d3fd12f\scratch\Printout1.pdf
esriJobMessageTypeInformative: Start Time: Tue Dec 22 16:55:18 2015
esriJobMessageTypeInformative: Executing (MyService): MyService "Feature Set" C:\Users\~\AppData\Local\Temp\MyService\MyService_122215_1630_gpserver\j25ae6595a46840369acebfbd9d3fd12f\scratch\Printout1.pdf
esriJobMessageTypeInformative: Start Time: Tue Dec 22 16:55:18 2015
esriJobMessageTypeInformative: Running script MyService...
esriJobMessageTypeInformative: Success reading the mxd from \\myserver\myshare\arcgisserver\directories\arcgissystem\arcgisinput\MyService\MyService_122215_1630.GPServer\extracted\v101\folder\MyService.mxd
esriJobMessageTypeError: Traceback (most recent call last): File "\\myserver\myshare\arcgisserver\directories\arcgissystem\arcgisinput\MyService\MyService_122215_1630.GPServer\extracted\v101\folder\GeoprocessingService_Test.py", line 20, in df = arcpy.mapping.ListDataFrames(mxd)[0] File "e:\arcgis\server\arcpy\arcpy\utils.py", line 181, in fn_ return fn(*args, **kw) File "e:\arcgis\server\arcpy\arcpy\mapping.py", line 1479, in ListDataFrames result = mixins.MapDocumentMixin(map_document).listDataFrames(wildcard) File "e:\arcgis\server\arcpy\arcpy\arcobjects\mixins.py", line 728, in listDataFrames return list(reversed(list(self.dataFrames))) File "e:\arcgis\server\arcpy\arcpy\arcobjects\mixins.py", line 695, in dataFrames return map(convertArcObjectToPythonObject, self.pageLayout.dataFrames) AttributeError: 'NoneType' object has no attribute 'dataFrames'
esriJobMessageTypeError: Failed to execute (MyService).
esriJobMessageTypeInformative: Failed at Tue Dec 22 16:55:19 2015 (Elapsed Time: 1.00 seconds)
esriJobMessageTypeError: Failed to execute (MyService).
esriJobMessageTypeInformative: Failed at Tue Dec 22 16:55:19 2015 (Elapsed Time: 1.00 seconds)
esriJobMessageTypeError: Failed.


texts: How can I import geometries from an invalid GML quickly?


texts: How can I add a timestamp in a Image Collection for producing images subtracted by a number in a monthly aggregation?

texts:
Following code produces a monthly aggregation for maximum 2 meters air temperature, in Kelvin degrees, by using a point reduce region with Image Collection "ECMWF/ERA5/DAILY". Results are identical to those obtained with daily series running an external developed python code. However, below code also has a commented portion destined to convert Kelvin to Celsius degrees. It was adapted based in Create time series and export it to csv GEE answer.

var fg_points = ee.Geometry.Point (-70.2166985,-38.5275513); //p18

// Set years and month
var startYear = 1981;
var endYear = 2014;
var years = ee.List.sequence(startYear, endYear);
var months = ee.List.sequence(1,12);
// load the image collection
var Daily = ee.ImageCollection("ECMWF/ERA5/DAILY")
    .select('maximum_2m_air_temperature');

// Convert  Kelvin to Celsius degrees
//Daily = Daily.map(function(img){
//  var date = img.get('system:time_start');
//  return img.subtract(273.15).set('system_time_start', date);
//});

// make monthly mean mosaics
// loop over the years and months to get summed monthly images
var byMonth = ee.ImageCollection(ee.FeatureCollection(years.map(function(y){
  var yearCollection = Daily.filter(ee.Filter.calendarRange(y, y, 'year'));
  var byYear = ee.ImageCollection.fromImages(
    months.map(function(m) {
      var maxImage = yearCollection.filter(ee.Filter.calendarRange(m, m, 'month'))
        .reduce(ee.Reducer.mean());
      var date = ee.Date.fromYMD(y, m, 1).format("MM_dd_YYYY");
      return maxImage.set('system:time_start', ee.Date.fromYMD(y, m, 1)).rename(date);
      //.set('month', m).set('year', y); // eventually set year and month
  }));
  return byYear;
})).flatten());

print(byMonth);

// filter the empty one out
var outputMonthly = byMonth.filter(ee.Filter.listContains('system:band_names', 'constant').not())
                    .sort('system:time_start').toBands();
print(outputMonthly);

var features = outputMonthly.reduceRegions(fg_points, ee.Reducer.first(), 30);
//print(features);

Export.table.toDrive(features,
"temp_max_mean",
"GEE_Folder",
"temp_max_mean_monthly_p");

When I uncomment code portion above referred:

.
.
.
// Convert  Kelvin to Celsius degrees
Daily = Daily.map(function(img){
  var date = img.get('system:time_start');
  return img.subtract(273.15).set('system_time_start', date);
});
.
.
.

for producing monthly series in Celsius degrees, I got following error message:

  ImageCollection.reduce: Can't apply calendarRange filter to objects
  without a timestamp

A similar message is obtained in a tasks Tab:

  Error: ImageCollection.reduce: Can't apply calendarRange filter to
  objects without a timestamp.

Recently, I found out this function:

// Convert  Kelvin to Celsius degrees
Daily = Daily.map(function(img){
  var start = ee.Date(img.get('system:time_start'));
  var end = ee.Date(img.get('system:time_end'));
  var label = start.format('YYYY-MM-dd').cat(' - ').cat(end.format('YYYY-MM-dd'));

  return img.subtract(273.15).set({label: label});
});

and it produces same error message. Nor does it add a timestamp.

So, How can I add a timestamp in a new Image Collection obtained from a native GEE Image Collection for producing images subtracted by a determined factor in a monthly aggregation?


texts: Python script to open many small images in QGIS

texts: I have access to about 1000 satellite images that covers my country. Each of them is associated to a smaller image that is used for browsing the real deal. I would like to write a script that will sweep through the folders and opens all the "browsing" images. I can program in python but I am not sure on how to connect a script like this to QGIS. Can anyone help me?


texts: Polyline to points ordered in cumulative distance from starting point

texts:
I have a problem with converting a line to points, where the points are ordered as the cumulative distance from the start of the line. I originally had 3 individual lines which I merged to make 1 continuous line using the Dissolve tool in ArcMap. However, it appears that the new line kept the features of the individual lines, because when I use the Divide function in Editor to create evenly spaced points along my 'merged' line, the points are ordered hap-hapzardly.

I am seeking a way to create points from that 'merged' line so that the points are ordered as incremental distances from the start of the line. It would be rather easy if the line were straight, however it meanders, therefore just ordering the points in ascending order (based on X or Y value) would not suffice.

Any ideas?


texts: Confusion regarding distance calculation in R (euclidean distance, "great circle distance")

texts:
OK so considering these two cases:

ln1<-SpatialLinesDataFrame(SpatialLines(list(Lines(Line(matrix(c(53.3604464,53.36062,-6.2424442, -6.242413),ncol=2)),ID="a"))),data=data.frame(dummy="a"),match.ID=F)
proj4string(pt1) <- CRS("+init=epsg:4326")

SpatialLinesLengths(ln1,longlat=T)*1000
SpatialLinesLengths(spTransform(ln1, CRS("+init=epsg:3857")),longlat=F)

ln2<-SpatialLinesDataFrame(SpatialLines(list(Lines(Line(matrix(c(15.43911,15.43914,47.00849, 47.00837),ncol=2)),ID="a"))),data=data.frame(dummy="a"),match.ID=F)
proj4string(ln2) <- CRS("+init=epsg:4326")

SpatialLinesLengths(ln2,longlat=T)*1000
SpatialLinesLengths(spTransform(ln2, CRS("+init=epsg:3857")),longlat=F)

I calculate the lengths of the lines (ln1 and ln2) in meters.

The first calculation being the "great circle distance" the second Euclidean distance.
Well I read that those distances should lie pretty close to each other when calculated for small distances.
That is true for the first case:

  Great Circle:

SpatialLinesLengths(ln1,longlat=T)*1000
[1] 19.51758

  Euclidean

SpatialLinesLengths(spTransform(ln1, CRS("+init=epsg:3857")),longlat=F)
[1] 19.63836

But in the second case the length difference is pretty great. I mean its over 40%...

  Great Circle:

SpatialLinesLengths(ln2,longlat=T)*1000
[1] 13.52404

  Euclidean

SpatialLinesLengths(spTransform(ln2, CRS("+init=epsg:3857")),longlat=F)
[1] 19.87276

Well I understand the difference between both methods (straight line vs. "as the crow flies" etc.) but reading (and understanding so) that the difference on small scale should not be to big. I worry seeing something like that...

Is it just because of the distance to the Equator? (What I can't imagine)
Is it a rounding issue?
Is my code wrong? (Well the same effect takes place using gLength(rgeos) or spDists/spDistsN1(sp) or any other distance calculation out there for R)

So what's going on here?


texts: Digitising animal positions in aerial photos

texts:
Digitizing objects from orthophotos can be done using QGIS. It is just a matter of establishing a points vector layer with the right attributes and working through the QGIS user interface.
Counting wildlife can involve looking at orthophotos and noting down the number of animals within a series of pictures.
In a survey we are planning we will fly transects with a drone while taking vertical images. The images are stitched together to several long GeoTIFFs. The next step is to have wildlife experts marking with a point every occurrence of up to 10 different species.
The wildlife experts need a straightforward tool which allows them to capture species type and position. Using QGIS or other high level GIS tools will require too much from our experts. The wildlife expert does not need to know where the observations were made. The use-case on their side would ideally look something like this:

Open images
Click on wildlife occurrence
Select species from drop-down list

This will generate a file/database with the resulting positions and species type.
We are looking for a survey tool to do this, so that out workflow is efficient and uncomplicated for our wildlife experts.


texts: PostGIS geom line string to display in Django map

texts:
I would like the WKB PostgreSQL geom data (extracted from a shapefile) in the below table to display in a django web app.  I have previously successfully mapped points to a Django web app from a Points-based PostgreSQL table, so I do have some rather basic knowledge.  My question is how to deal with the geom data in WKB format.  In particular, how to write an HTML template to display the geom WKB data onto a Django map?  I have studied some Python but have limited knowledge of HTML and JavaScript.
What additional info would be needed to answer my request? I was unable to find another question that asks this question (probably too simple).
My end goal is to display a map similar to the one on this page:
http://ryanpeek.org/2016-07-01-pctmap/


texts: Is there any way to make an ArcMap mxd doc into a simple exe file so that one who does not have ArcGIS can view my work?


texts: Two different layers of the same area have to be united to form a single layer

texts:
I have two different layers containing different polygons and I need to combine them to see the differences over time.

So the output has to be a layer containing the attribute of interest of layer one, the attribute of interest of layer two and the attribute of the change (example : entity 1 - entity 2 OR entity 1 - entity 1 IN THE CASE OF NO CHANGE)

How can I do that using QGIS?


texts:
Is it possible to order grid codes in a shapefile visually by building a query in ArcMap?  For example in symbology I have grid codes 1, 2, 4, 5 in a shapefile and I want to make it so that when the symbology is displayed on the map itself, grid code 2 is always onto grid code four.

I can currently make grid code 4 go away by
"gridcode" = 1 OR "gridcode" = 2 OR "gridcode" = 5

But I want to maintain grid code 4, and if it over laps grid code 2 I want grid code 2 to be on top of grid code 4 when I am visually looking at the map

I'm using the most recent version of ArcMap.

Ordering the grid codes and applying it in symbology currently does not solve the problem.


texts: Recursive query PostGIS - quadgrid function

texts:
I'm just getting started on the WITH RECURSIVE function in PostGIS. I'm trying to recurse through a list of ids, passing them to through to the quadgrid function (credit goes to Mark Wynter at dimensionaledge). The quadgrid function only accepts a single geometry row as input.

The ids are returned via the query:

SELECT parent.gid FROM proc.gfn_rdsid78 parent

The ids are then called in the WHERE clause:

SELECT the_geom::geometry(Polygon, 4326) as geom
FROM DE_RegularQuadGrid((SELECT child.geom FROM proc.gfn_rdsid78 child WHERE child.gid = 1276643),
                    'rds_autogenerated.rdsid78_pointseed_urbanflood_singlepoint',
                    'geom',
                    3,
                    10)

What would be the optimal of getting back my unioned geometries?


texts:
I am running this tool: https://www.arcgis.com/home/item.html?id=a269d03aa1c840638680e2902dadecac Mxdperfstat- basically it records the loading time for drawing the layers so you can optimize data and map docs before publishing them.

Everything seems to work fine but my output file is just a .txt with a bunch of source code which means nothing to me (I can't understand). I tried to open it by double clicking and it opens with notepad Like this img:

and when I right click and open with a web browser it just opens an empty page nothing is displayed (all white).

The results I am expecting are referred to this page:
http://egis3.lacounty.gov/eGIS/2012/02/21/tune-your-map-data-before-publishing-mxdperfstat/

Anyone have any ideas why my output is coming out like this?


texts: Use field calculator to reclassify values from three different fields to populate a new field AND determine overlap with different raster

texts:
I want to reclassify values of 3 different fields AND determine if one of those fields is located inside a different polygon to populate the values of my new field. Everything is raster.

I created a new field (WScore) that I want to populate with numeric values (double) based on the values of 3 different fields - Age_Databa, Riparian_S, and SITE_PHY.

My parameters:

If Age_Databa is 200 or >200 and Riparian_S is not 1 (any number but 1), then WScore should = 5

If Age_Databa is >200 and Riparian_S = 1, then WScore should = 2

If SITE_PHY is <70 and it is in a hydric soil polygon (a raster in my gdb), then WScore = 4

Before realizing the last part about the needing SITE_PHY to be inside the hydric soil polygon, this is what I started with. Pretty sure my syntax is wrong. Any help would be fantastic.

def Reclass(Soils):
global out
 if (Age_Databa <= 200) and (Riparian_S > 1 or Riparian_S < 1):
  out 5
 elif (Age_Databa <= 200) and (Riparian_S == 1):
  out 2


texts: QGIS 1.8 on Ubuntu 12.04 doesn't transform CRS properly

texts:
I have a weird QGIS-Problem that only occurs in Ubuntu 12.04, but not in Windows 7. I've also tried this on different Ubuntu-Computers.

When I try to transform a layer with WGS84 to a different projection type (example Austrian Gauss Kruger) the result is shifted exactly 300m to the East to where it should be. This happens with both "on the fly transformation" and when I save the Layer in a different CRS!

When I save a WGS84-Layer in Gauss Kruger and I open both the new and the old layer in Ubuntu, they overlap perfectly (with on the fly-enabled). When I switch to Win7 and open both layers (the original and the one I created in Ubuntu) I can even see the 300m gap!

And ideas? Could this be a bug? Thanks!

(I am using the same EPSGs in Windows and in Ubuntu - I can also see/not see the problem, when I open the same project in Ubuntu/Windows.)


texts: QGIS not recognizing special characters and problem with intersection?

texts:
Started playing with QGIS recently and I have encountered a problem. I have a layer with administrative areas of my country and when adding it the encoding is set to UTF-8. When I open its attribute table the feature names are jumbled and special characters are replaced with symbols.
Also, when I try to use intersection with another layer i get the following message: Input layer A contains invalid geometries (feature 56). Unable to complete intersection algorithm.

I tried with 32 and 64 bit version of QGIS 2.18.9. I'm running Windows 10.
The catch is, a friend of mine is using QGIS 2.14.2 on Windows 10, and didn't have any problems with encoding (her attribute tables are fine) or intersection (works without any error messages) while using the same SHP files.

What could be the problem?


texts: Projecting shapefile with transformation using OGR in python

texts:
I want to change coordinate system of a shapefile from British National Grid to WGS84 (in order to be able to convert it to KML later on). To get the best results of projecting I would like to use OSGB_1936_To_WGS_1984_NGA_7PAR transformation (or later on define different transformation depending on location in UK).

The only description of this transformation that I know comes from ArcMap and looks like this:

Coordinate Frame - dx=446,000000 dy=-99,000000 dz=54,000000 rx=-0,9450000 ry=-0,261000 rz=0,435000 s=-20,892700

.. and I have no clue how to use it with osr.CoordinateTransformation(BNG, WGS84)

Any help highly appreciated!!

Jarek


texts: Why does the layer list context menu not open anymore?

texts:
I cannot access layers properties (i.e. change color, filter by query, anything).
The mouse right click doesn't do anything.

Is there a way restoring qgis to default?

Ubuntu 12.04, QGIS 1.8


texts: How to fix 'GeoSeries' object has no attribute '_geom'

texts:
I'm working on a python script whose goal is to detect if a point is out of a row of points (gps statement from an agricultural machine). To do this, my idea is to make a buffer around the 10 points around considered point (five before and five after). After that detect if my point is in the buffer. I'm stuck with  'GeoSeries' object has no attribute '_geom'

Here is my code for this part:

for i, row in src.iterrows():

    # Tague les points trop éloignés du rang en cours
    imin, imax = getMinMax(i, nl, 5) # Retourne les 5eme voisins avant et apres

    buddies = src.loc[imin:imax, 'geometry'].drop([i], axis=0) # DF des voisins de i avec uniquement la geometrie
    buddies_buf = buddies.buffer(ird*1.2) # Tampon avec un rayon egal a la distance inter rang + 20%

    buddies_buf = buddies_buf.unary_union # Fusionne tous les tampons
    buddies_buf = gpd.GeoDataFrame(gpd.GeoSeries(buddies_buf)).rename(columns={0: 'geometry'}).set_geometry('geometry') # Formate le tampon en GDF
    buddies_buf.crs = crs_work # Attribut le crs

    print((row.geometry).within(buddies_buf.geometry))

Somebody can help me?

I don't know if my idea is the best for this kind of detection, so if you have another best idea...


texts: Aligning Web Map Layers from Bing Maps and Google Maps in QGIS?

texts:
Why doesn't Bing and Google imagery "line up" in QGIS.

Basic steps to replicate the problem:

Open new QGIS project and set project CRS to EPSG 3857 (Pseudo Mercator).
Open Bing and Google Satellite layers using QuickMapServices plugin.
Toggle layers over an urban area, and watch the buildings shift.

The Google imagery seems to line up better with vector layers projected on the fly, so I am assuming that there is a problem with the Bing layer, and how I am using it.

Should I use a different CRS with the QuickMapServices Bing satellite layer?


texts: How can I change the pixel value from 0 to 10 on QGIS?


texts: How to average daily data to create a weekly estimate in Google Earth Engine?

texts:
I need to get weekly averages of AOD for CA for 15 years. The filter.calenderRange does not have a "week" options. I have filtered the data to the day level but am unsure how to average for the week. Here is an example of my code for the year 2000 : `

         var filtAOD00=AOD.filterDate('2000-03-01','2000-12-31');
         var CA=US.filterMetadata('NAME_1','equals','California');

          //function to filter out images with missing data

         var imgColPerday = ee.ImageCollection.fromImages(
                  ee.List.sequence(1,366).map(function(day){
         var image = filtAOD00.filter(ee.Filter.calendarRange(day,day,'day_of_year'))
         .select('Optical_Depth_047').mean();
         image=image.addBands(ee.Image.pixelLonLat());
         return image.set('day_of_year', day)// set the day
         // filter out empty images
         })).filter(ee.Filter.listContains('system:band_names', 'Optical_Depth_047'));

         // then apply sampleRegions
           var data =  ee.FeatureCollection(imgColPerday.map(function(image){

         // calculate data
          var AOD00Census = image.sampleRegions({
          collection: CA,
          scale: 1000,
          }).map(function(feat){
          return feat.set('day_of_year',image.get('day_of_year'));
          });

           return AOD00Census;
            })).flatten();

       var polyOut = data.select(['.*'],null,false);
       Export.table.toDrive(polyOut, 'AOD_CACensus1km_00');

`


texts: java.lang.NumberFormatException: For input string: "-?"

texts:
I'm getting an error in Maxent which arises from one of the output files created by Maxent, in this case the _ASCII_novel.asc at row 159 column 444. However, the data in this position appears to be OK.

I saw a similar question where the error noted is For input string: "Z^,,," which seemed to be caused by an invalid file type (GeoTiff not .asc), but Maxent didn't seem to have a problem using my .asc files before the output file produced an error, so I don't think this must be the problem.

I have run Maxent without my bias file, and this avoids the error.

If it's of any help the rest of the error reads:

at sun.misc.FloatingDecimal.readJavaFormatString(Unknown Source)
at sun.misc.FloatingDecimal.parseFloat(Unknown Source)
at java.lang.Float.parseFloat(Unknown Source)
at density.MyScanner.nextFloat(MyScanner.java:46)
at density.LazyGrid.readRowAsc(LazyGrid.java:166)
at density.LazyGrid.readRow(LazyGrid.java:74)
at density.LazyGrid.eval(LazyGrid.java:87)
at density.LazyGrid.hasData(LazyGrid.java:81)
at density.Grid.minmax(Grid.java:156)
at density.tools.Novel.makeNovel(Novel.java:134)
at density.tools.Novel.go(Novel.java:99)
at density.Runner.makeNovel(Runner.java:618)
at density.Runner.start(Runner.java:563)
at density.GUI$7.construct(GUI.java:342)
at density.SwingWorker$2.run(SwingWorker.java:131)
at java.lang.Thread.run(Unknown Source)


texts: Performing route analysis in R based on OSM transport network?

texts:
I'm analysing the likely paths taken by commuters on their way to work, based on Census flow data, job density distributions and the transport network.

At the moment I'm using pgRouting and the pgRouting add-on for QGIS (described here) for my analysis:

However, I'd like all the steps I take to be full replicable (see excellent article on this here), then I'd prefer to do my analysis from the command-line and save the script files for future reference.

R-geo packages fulfil this requirement well, as they are free, lightweight, and stable: download a tiny script file and you should be able to re-do all the things I have done (see here for an unrelated example of this on Rpubs).

The R packages I've identified for doing this are:

gdistance published December 2012, looks promising, but no sign of compatibility with osm transport networks
e1071 seems rather abstract and possibly overkill for what I want to do
pairdist.lpp function in spatstat package looks promising, but has little documentation

Before diving in to these, thought it worth asking people who may have faced similar dilemmas.

Is it possible to do my routing analysis in R, or would it be best to learn PostGIS properly and do it all in pgRouting?


texts: How to clip a raster and exclude pixels that overlap the edge of the mask layer

texts:
I recently posted an answer to Zonal stats for complete pixels qgis that required clipping a raster in a way that excludes pixels that fall partially outside the mask layer. Here's an illustration:

The best method I could figure out is as follows:

Polygonize the raster.
Use Select by location tool to select features in Vectorized layer
that fall within the polygon layer.
Clip  raster by mask layer using the selected features of the
Vectorized layer as the mask layer. Don't crop the extent.

But this method seems like a lot of steps, and polygonizing a large raster can take a long time.

Is there a better method of clipping a raster and excluding cells that overlap the edge of the mask?

By "better" I mean fewer steps and/or requiring less processing time.


texts: How to dissolve polygons that overlap both on each other in QGIS?

texts:
I have two layers of polygons with landcover classes(meadows, cropland, forest, ect) They overlap in both directions, i.e like a zipper. I need to dissolve them by one attribute (Class) so I can have the total area of each polygon... but in the overlapping sections sometimes the landcover classes don´t match.

If I try UNION, the algorithm creates new features for this overlapping regions, but also creates another AREA column and in some cases I have values for both columns simultaneously (I have now 2 AREA columns: one for each layer I "unioned").

I also tried the merge layers of the MMQGIS plugin, but the attributes are not summed up. And dissolving by class just creates an empty layer.

(The image is the union of the two layers. The arrows indicate which layer is overlapping over the other. Both images are corrected in terms of projection)

Does someone knows how to do something like this?
I work on QGIS 2.12.3 on a PC.