Friday 31 January 2020

references - Seeking ArcGIS Desktop tutorials?



I'm looking for a good way to learn more ArcGIS techniques than I came across in my undergraduate Geography degree to really use the power of the software.


Browsing around and asking questions on this site has helped quite a bit e.g. I did not realised that you could link dataframes together, or even that you could have multiple data frames ([see How to provide automatic overview map in ArcGIS Desktop?).


Similarly, I suspect there are lots of ArcGIS features that could really make my life easier, but that I just don't know about e.g. geodatabases.


What I'm really looking for are tutorial recommendations, in any medium (books, videos, websites, whatever). Obviously free or cheap are better, but I do have a budget for books as part of my PhD. As I probably hinted at above, I'm more interested in general usage of ArcGIS rather than specific tools (like DEM processing or route-finding), and something that is relatively structured would be good, so that I can work through it in a sensible order.


I know this is a bit of a tall order, but any ideas?



Answer






Ensure you have ArcTutor installed (comes with ArcGIS Desktop install)


Start doing the tutorials.


Id also advise paying a trip to the ESRI Virtual Campus. - Plenty of free training on there to get you started.


Go to video.esri.com for inspiration.


arcgis desktop - Grid-like lines appearing after using Curvature tool


I have been attempting to use the Spatial Analyst Curvature tool with mixed results. I originally tried to use it with a raster that was projected using a 1983 Teale Albers projection and the results were black, indicating the lowest values. I then transformed the raster into UTM projection and ended up with grid-like lines on the raster which have values that are obviously not accurate. Does anyone know why these lines have appeared or what I'm doing wrong. Thanks.


enter image description here




Hide unselected features with pyqgis


In the context of QGIS/PyQGIS, I would like to hide the unselected features after performing, for example, a simple selection.


feature = layer.getFeatures().next()
layer.setSelectedFeatures([feature.id()])

After running this commands a feature becomes highlighted. Besides, I'd like to hide the unselected features. How can I do that?


I believe this is possible because the Query Builder performs a filter that hides the undesired features.



Answer



I just made it. In fact, I should use the setSubsetString method. Then I built a a query string inside a loop in order to look like this "gid IN ('8816','8836','8839','8864')". Then, when you run the setSubsetString method it hides the features that are not specified in the query string.



#Do any spacial query which returns some features
features = lyrPnts.getFeatures(QgsFeatureRequest().setFilterRect(geomPoly.boundingBox()))

#build the query string
strsel="gid IN ("
for feat in features:
if feat.geometry().intersects(geomPoly):
strsel=strsel+ "'" + str(feat.id()) + "',"

#Closes the string with a parenthesis

strsel=strsel[:-1] + ")"
#Hides the undesired features
lyrPnts.setSubsetString(strsel)

Nice job


Alex


Seeking municipal-level administrative boundary data for Poland?



I am looking for municipal-level administrative boundary data (NUTS 5) for Poland. Since it is public information, I thought the data would be readily available, but I can't seem to find it anywhere. I have checked the usual sources. GADM has data down to the county (powiat) level, and Natural Earth seems to only have data down to the provincial level. I have looked at many of the Polish government's sites, but while I can find plenty of GIS data viewers, I can't find a way to download the administrative boundary data I am interested in. To make matters more complicated, I would really like historic data, but I would settle for current NUTS 5 level data.



Answer



You might try contacting Francis Harvey, one of the Co-Directors of the MGIS program at the University of Minnesota. He does a lot of GIS research in Poland, and is probably very well aware of what is and is not publicly available.


arcpy - Defining parameter descriptions for Python Toolbox help?


I'm trying to create some Python Toolboxes for our ArcMap application (e.g. MyTool.pyt)


I can see that the help text is defined with the classes self.description attribute.


However, once I run the program, and click into any of the parameter fields, the help/description text goes empty. I would like to be able to provide the description field for each parameter. How is this accomplished?


After some responses, I see that via the 'Item Description' right-click context menu there are many fields that may be populated. Is there a 'pythonic' way to do this? That is, just by embedding some attributes in the .pyt file classes?


For instance, in the .pyt toolbox definition you have the Toolbox class:


import arcpy


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

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



class MyNiceTool(object):
def __init__(self):
"""Define the tool (tool name is the name of the class)."""
self.label = "My Tool Class"
self.description = """
A description that shows up in the help context side pane when the tool is launched.
"""
self.canRunInBackground = True


def rest_of_required_methods....

From the self.description string the Tool dialog help window will display this text. However, what I am wanting to do is have a 'description' embedded in my code for each parameter also, so that when the tool is launched, and the user clicks into a parameter field, the parameter description is shown. If I were to do this using the 'Item Description' method referenced in the replies below, I would edit the fields Dialog Explanation under the Syntax section for each parameter... I guess.



Answer




I can see that the help text is defined with the classes self.description attribute.



This is where you are going wrong. In the help page Documenting a tool in a Python toolbox it says:




For Python toolboxes, the documentation for the toolbox and tools are stored in .xml files that are associated with the toolbox and tools by name. The help for each tool will be stored in a separate .xml file.



This means that you cannot set help text from within the .pyt file itself. This makes sense when you consider that the help text is not plain ASCII but rich text that can include formatting, bullets, and images.


Fortunately, Python supports for reading and writing XML, so you should be able to dynamically edit the help text from a separate script.


How to set the Identify Feature tool in QGIS 2.4?


The Identify Results window does not appear when using the identify feature tool in QGIS 2.4. I guess there's a check box somewhere I need to check/uncheck, but the settings of the tool changed compared to the previous version (see attached printscreen). QGIS 2.2 on the left and 2.4 on the right



Answer




Right click on the toolbar and select Identify Results


The options for the different select modes has been moved to the dock widget.


r - ggmap: create circle symbol where radius represents distance (miles or km)



I am trying to add a circle around a given latitude and longitude point where the circle size (radius) reflects distance in miles or km.


I managed to come up with this example using ggmap, where I can add a circle. but what do I need to do to have it represent distance. I would ultimately like the bounds of the circle to represent 20 miles from the center.


Code below:


library(ggmap)
library(ggplot2)

d <- data.frame(lat=c(33.79245),
lon=c(-84.32130))

emory <- get_map("Atlanta,Georgia", zoom=12)


p <- ggmap(emory)
p <- p + geom_point(data=d, aes(x=lon, y=lat),
color="red", size=20, alpha=0.5)
p


distance - Clustering spatial data in R?




I have bunch of data points with latitude and longitude. I want to use R to cluster them based on their distance.


I have already taken a look at this page and tried clustTool package. But I am not sure if clust function in clustTool considers data points (lat,lon) as spatial data and uses the appropriate formula to calculate distance between them.


I mean I cannot see how they differentiate between spatial data and ordinal data. I believe the distance calculation between two points on map (spatial) and two normal numbers is different. (Is it not?)


Also what happens if I want to consider a third parameter in my clustering?


Like say if I have (lat,lon) and one other parameter.


How is the distance calculated?


The other problem I have with clustTool is that it is designed with a GUI in mind. I don't know how I can skip the GUI overhead in the library because I don't need it.


What options do I have in R for cluster analysis of spatial data?




Thursday 30 January 2020

arcgis desktop - How to get size of file geodatabase feature class on disk?



Is there a simple way to determine the amount of hard drive space a feature class takes up?


I feel like I'm missing something simple, but I don't recall a method that does so.


You would think that right-clicking a dataset in ArcCatalog would do it for you.


The best I've ever been able to do was get the size of the entire gdb by looking in Windows Explorer.


If there isn't a way to do this in the ArcCatalog UI, I'd still be interested in a method to do so in code.


Any ideas?



Answer



If you're looking for an ArcObjects way to get it, then How to programmatically determine the size of a feature class in a file geodatabase? will provide that -- otherwise, you can enable the Size column in the Customize menu -> ArcCatalog Options -> Contents tab:


Enabling size column in ArcCatalog Options -> Contents tab Size column enabled reveals feature class is 50.69 MB


This works on file geodatabases but not on SDE geodatabases (in that case you could use some DBMS-specific queries to determine it though). It does not work on personal geodatabases.



qgis - Calculate population within radius


I have population data matched with polygons and a set of points that I have drawn buffers around. Assuming that the population within each geographical area, i.e. the polygons, is evenly distributed, how do I calculate the population within each buffer? I need to do this in QGIS.


(1KM is the buffer and befolkning_i grkrets is the polygon.)



enter image description here




pyqgis - python script in QGis - help using gdal.Polygonize


I want to create a very simple script in python that polygonizes a raster. Based on python cookbook + forums search I wrote the following code


#import des modules
from osgeo import gdal, ogr
from qgis.core import *
from qgis.gui import *
import sys

from PyQt4.QtCore import *

##srcRaster=Raster
##Result=output vector

# 1. Open Raster and get band 1
src_ds = gdal.Open(srcRaster)
srcband = src_ds.GetRasterBand(1)

# 2. Create output datasource

Result = QgsVectorLayer("Polygon", "temporary_points", "memory")
Result.startEditing()
Result.dataProvider().addAttributes([QgsField("ClassNo", QVariant.String)])
Result.updateFields()
gdal.Polygonize( srcband, None, Result, 0, [], callback=None )

which ends up with the error "in method 'Polygonize', argument 3 of type 'OGRLayerShadow *". The truth is, I don't really understand what should be the syntax of gdal_polygonize in a qgis python script (command line can be found here).


Could anybody please help ?




arcgis 10.3 - Connecting to Oracle User Schema Geodatabase using ArcObjects?



I've set up an Oracle user Schema geodatabase and loaded some data into it. I now want to connect to it programmatically using ArcObjects, the problem is I don't see how. Where can I set the newly created schema name in my connection property set, or how do I change the schema value after I have connected? In ArcCatalog it is done in the Geodatabase connection properties window : enter image description here


The code below works fine connecting to the default SDE schema. What do I have to change to connect to a user geodatabase? I have tried examining IDatabaseConnectionInfo2 and IWorkspaceProperties but I don't see where I can set or change the schema that is connected to.


 public static IWorkspace ConnectToTransactionalVersion(String dbclient, String dbConnProp, String user, String password, String database, String version, String authentication)
{
IPropertySet propertySet = new PropertySet();
propertySet.SetProperty("DBCLIENT", dbclient);
propertySet.SetProperty("DB_CONNECTION_PROPERTIES", dbConnProp);
propertySet.SetProperty("DATABASE", database);
propertySet.SetProperty("USER", user);
propertySet.SetProperty("PASSWORD", password);

propertySet.SetProperty("VERSION", version);
propertySet.SetProperty("AUTHENTICATION_MODE", authentication);
Type factoryType = Type.GetTypeFromProgID(
"esriDataSourcesGDB.SdeWorkspaceFactory");
IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance
(factoryType);
return workspaceFactory.Open(propertySet, 0);
}

I tried adding "SCHEMA" as a property in propertySet and it gets ignored, I also tried setting the version property to the new user schema default version but it errors with "version not found"





qgis - Line to polygon conversion


Currently we are working in a geological project..in which all the Contacts, Linear & some other geological features have to be captured.


The thing is most of the features have to be Captured as "Line" and then it has to be converted to "Polygon". So that we don't have to repeat the digitization for features like Contact....


1) Is it not at all possible in Qgis to Snap a line with the same line ? like that of in AutoCAD, Mapinfo..?? (Image 1)


2) I've captured a polygon as two separate lines...and when i use "Line to Polygon" (Enclose in the case of Mapinfo)....it actually creates two polygons instead of one single entity. Is it a bug ? (Image 2)



I am trying to convince some of colleagues (and the management as well) to switch over to QGIS..atleast in few systems...So i am trying to figure out the available functionalities in QGIS as that of MapInfo. So any sort of feedback will be of great help.


enter image description here


enter image description here


Update: Polygonize command does the trick of converting the two separate lines into one single polygon feature...(Image 2).



Answer



You can snap a line to itself, but only after saving.


My workaround is to digitize the line, but set the last point a bit offset to the first point, then save, and then move the last point to the first.


As a consequnce, I start digitizing a new closed line with a point that is already created by another adjacent polygon if possible.


Wednesday 29 January 2020

Splitting line at point positions using QGIS?


I am a rooky with QGIS and i am stuck from a basic problems.



I have a polyline and i have some points (2) on the polyline.


I would like to know if there is a tool or script that i could use to split my original polyline into 3 polylines? (Starting point to point A, Point A to B and B to ending point of the polyline)?


I am looking for an automatic process rather splitting the line manually.




united states - Geocoding USA addresses that cannot be sent over internet?


For a one-time project, I need to geocode a few thousand addresses. In the past I've used various online resources for this kind of thing (e.g., Google Maps API), but the addresses I'm working with have to be kept confidential - which means no sending it over the Internet, unless there's some iron-clad guarantee of privacy. What other options do I have?




dem - PostGIS interpolation / triangulation options



I'm looking for an open source workflow to automate DEM construction. We have a series of sites that fall within a LIDAR dataset. We want to create site specific DEM's for each site, and we're looking to automate the process.


So far, we have automated:



  • data loading into PostGIS (LIDAR and site points)

  • creation of site boundaries (combination of st_buffer and st_envelope)


Now we're looking for options to interpolate the data subsets at each site and export them to Surfer7 grid files.


Currently, we're querying the data in the PostGIS database using QuantumGIS, exporting to csv files, then manually importing and gridding the data in Surfer7. Hopefully we can automate this as well.


So for this use case, we'd like to triangulate our lidar data to a DEM. Other potential cases we can think of involve different interpolation methods - so if there are options for inverse distance and kriging - we're interested!


This is very much a learning exercise for us - we're working in baby steps!





GDAL - Perform Simple Least Cost Path Analysis


I am investigating methods to perform a simple least cost path analysis with gdal. By simple, I mean using the slope of a dem as the only cost factor.


I would prefer to do using the python or .net bindings, but will take anything. Can anyone suggest any good tutorials or the like?




python - OGR Layer Intersection


I am using what I would describe as the most primitive method of doing intersection with OGR. The short script below describes how I do it. Is there a best way of doing this?


   from osgeo import ogr, osr

shp1 = ogr.Open(file1)
shp2 = ogr.Open(file2)

SpatialRef = osr.SpatialReference()
SpatialRef.SetWellKnownGeogCS('WGS84')
# Create dst file here
dstshp = ogr.CreateDataSource('SomeFilename.shp')
dstlayer = dstshp.CreateLayer()
# define its attribute fields for dstlayer and create them


layer1 = shp1.GetLayer(0)
layer2 = shp2.GetLayer(0)

for feature1 in layer1:
geom1 = feature1.GetGeometryRef()
attribute1 = feature1.GetField('FieldName1')
for feature2 in layer2:
geom2 = feature2.GetGeometryRef()
attribute2 = feature2.GetField('FieldName2')
intersection = geom2.intersection(geom1)

dstfeature = ogr.Feature(dstlayer.GetLayerDefn())
dstfeature.SetGeometry(intersection)
dstfeature.setField(attribute1)
dstfeature.setField(attribute2)
dstfeature.Destroy() # and other features must be destroyed too

Answer



There are some errors in your script but it is not the most important problem:


You cannot create a valid shapefile without specifying the geometry of the layer:


driver = ogr.GetDriverByName('ESRI Shapefile')
dstshp = driver.CreateDataSource('SomeFilename.shp')

dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon)

And you don't know a priori (upfront) the geometry of the resulting intersection layer. The intersection of two polygon layers is different from the intersection of a polygon layer and a polyline layer for example.


For that, you can get the geometry of the intersection by:


For example (with two polygons shapefiles):


layer1.GetGeomType()
3 # -> polygon
# create an empty geometry of the same type
union1=ogr.Geometry(3)
# union all the geometrical features of layer 1

for feat in layer1:
geom =feat.GetGeometryRef()
union1 = union1.Union(geom)
# same for layer2
union2=ogr.Geometry(layer2.GetGeomType())
for feat in layer2:
geom =feat.GetGeometryRef()
union2 = union2.Union(geom)
# intersection
intersection = union1.Intersection(union2)

print intersection.GetGeometryName()
'MultiPolygon'

At this stage, you can save the resulting geometry to a shapefile (without the fields of the original layers):


dstshp = driver.CreateDataSource('SomeotherFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbMultiPolygon)

But if you want to use your script (a MultiPolygon is a collection of Polygons):


driver = ogr.GetDriverByName('ESRI Shapefile')
dstshp = driver.CreateDataSource('SomeFilename.shp')

dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon)
for feature1 in layer1:
geom1 = feature1.GetGeometryRef()
attribute1 = feature1.GetField('FieldName1')
for feature2 in layer2:
geom2 = feature2.GetGeometryRef()
attribute2 = feature2.GetField('FieldName2')
# select only the intersections
if geom2.Intersects(geom1):
intersection = geom2.Intersection(geom1)

dstfeature = ogr.Feature(dstlayer.GetLayerDefn())
dstfeature.SetGeometry(intersection)
dstfeature.setField(attribute1)
dstfeature.setField(attribute2)
dstfeature.Destroy()

Don't forget to define the fields before (look at Python GDAL/OGR Cookbook:Vector Layers). And it is much easier with the module Fiona


leaflet - Cartodb Layer source object not loading


I have tried to figure out why the layer is not showing up on the map. My code is as follows.


var markersLayerSource = {
user_name: 'someaccount',
type: 'cartodb',
sublayers: [{
sql: baseSql,

https: 'force https',
cartocss: cartoCss,
interactivity: 'cartodb_id'
}]
}

cartodb.createLayer(map, markersLayerSource, {https: true})
.addTo(map)
.on('done', function (layer) {
console.log(layer);

layer.getSubLayer(0).setInteraction(true);
layer.getSubLayer(0).on('featureOver', function(e, pos, pixel, data) {
// print data to console log
console.log("Event #" + data.cartodb_id + ", name " + data.name + ", max population: " + data.pop_max);
});
//layer.getSubLayer(0).show();
//cdb.vis.Vis.addInfowindow(map, layer.getSubLayer(0), ['cartodb_id']);
});

Any one with pointers on why i specifically need to have cdb.viv.Vis.addInfowindow inorder for the layer to show up on the map? Am I doing something wrong if one needs to render without the addinfowindow?




Answer



Please check all the variables you have introduced within your markerLayerSource (the console will tell you what is wrong anyway). I have replicated your code in this working example. But I have used a different methodology to add the tooltip:


var tooltip = layer.leafletMap.viz.addOverlay({
type: 'tooltip',
layer: layer,
template: '

{{name}}

{{pop_max}} inhabitants

',
width: 200,
position: 'bottom|right',
fields: [{ name: 'name' }, { pop_max: 'pop_max' }]
});


$('body').append(tooltip.render().el);

So first, check that all your variables are well declared (baseSql and cartoCss) and secondly, try to use this methodology to add the tooltip. This is a screenshot of the result:


enter image description here


raster - Modifying DEM elevation data in QGIS using a polygon


I've tried searching through the existing questions but am struggling to find a succinct answer to what I am trying to achieve.


For background, I am running a flood model in TUFLOW to help support a residential development on greenfield land. Part of the site has been shown to flood in a given event and I now want to reprofile areas of the land to allow additional flood storage. The hope is that with the reprofiled land, the flood extents will be greatly diminished and so there will be more land available to the development.


I have drawn polygons in QGIS over the exact areas that will need to be lowered by specific amounts i.e. the land in one area varies from 66.42-66.73m AOD. This will now be a uniform 66.10m AOD.


Is there a way of simply modifying the existing raster DEM file using the polygons?




python - -9999 (no data value) becomes 0 when writing array to GDAL memory file



I am trying to reproject and resample (1000m ->> 30m) a raster image (shape 238X406) using gdal.ReprojectImage() function in Python. The input file has -9999 cells as no data value. The result is an array (shape 10834x15055). The data type is float32.


When I write the result to a geotiff file, everything is expected.No data value is set to -9999 and output array has -9999 cells. However, when I write the result to a gdal memory file to save some time (25 seconds shorter processing time, down from 105 seconds to 80 seconds), this time all -9999s (no data value) become 0 (zero) in the output array of memory file. Both results are exactly same, but GeoTIFF file has -9999, memory file has 0. Even though memory file has no data value of -9999, but the output array is initialized to 0.0 instead of -9999.


I am using the same code to produce the results and the only difference is that I call memory driver when I want to write result to memory file (driver=gdal.GetDriverByName("MEM")), and GeoTIFF driver when I want to write result to GeoTIFF file (driver=gdal.GetDriverByName("GTiff"))


My gdal version is '1100000' from gdal.VersionInfo(). My operating system us ubuntu 12.04 LTS "precise".


# Create a file
outFileRead=driver.Create(outFilePath,X,Y,1,dataType,options)
print inFileRead.GetRasterBand(1).GetNoDataValue()
print inFileRead.GetRasterBand(1).ReadAsArray()
print inFileRead.GetRasterBand(1).ReadAsArray().shape
# Reproject

gdal.ReprojectImage(inFileRead,outFileRead,inProjection,outProjection,reSamplingType)
print outFileRead.GetRasterBand(1).GetNoDataValue()
print outFileRead.GetRasterBand(1).ReadAsArray()
print outFileRead.GetRasterBand(1).ReadAsArray().shape

the result of prints


-9999.0


[[-9999. -9999. -9999. ..., -9999. -9999. -9999.]


[-9999. -9999. -9999. ..., -9999. -9999. -9999.]


[-9999. -9999. -9999. ..., -9999. -9999. -9999.]



...,


[-9999. -9999. -9999. ..., -9999. -9999. -9999.]


[-9999. -9999. -9999. ..., -9999. -9999. -9999.]


[-9999. -9999. -9999. ..., -9999. -9999. -9999.]]


(406, 238)


-9999.0


[[ 0. 0. 0. ..., 0. 0. 0.]


[ 0. 0. 0. ..., 0. 0. 0.]


[ 0. 0. 0. ..., 0. 0. 0.]


...,



[ 0. 0. 0. ..., 0. 0. 0.]


[ 0. 0. 0. ..., 0. 0. 0.]


[ 0. 0. 0. ..., 0. 0. 0.]]


(15055, 10834)


This shows even if no data value is set to -9999 in memory file, the array is initialized to 0.0 instead of no data value. Array in GeoTIFF is initialized to no data value, -9999, as expected. I think this is a bug in the memory file.




qgis - OpenStreetMap: CRS not matching


I have different layers with different CRS, I think therefore the layers don't match. But I don't know how to fix this. Here is what I did (QGIS 3.0.1 on Debian Buster).


I connected to OpenStreetMap Tile Server https://a.tile.openstreetmap.org/{z}/{x}/{y}.png to request XYZ Tiles. This is the first layer.


The tile server uses the CRS WGS 84, Pseudo Mercator Projection EPSG 3857.



Now, I tried to built my own map, by doing a quick query with QuickOSM. Key is building, querying Extent of the map canvas (coordinates are 300.000, 6.100.000, Scale is 1:60.000, which is in France, somewhere south of Paris).


This gives me three additional layers, all called OsmQuery, see image.


enter image description here


The additional layers use CRS WGS 84 EPSG:4326. If you look at the map, you see the buildings belonging to Chateau-Landon (south) near La Madeleine-sur-Loing (north). The destinations are about 7km apart.


My thinking is: "The error might result from different projections. However, 7km is a lot. I find it kind of strange that OpenStreetMap uses two different CRS." But I'm quite a newbie...


What can I do to get an exact match of all the layers?


My QGIS version is as follows:


enter image description here




arcgis server - Split polyline at specific distance


I am working on a project that uses SQL server for geometry datatype and ArcGIS Server and asp.net website and entity framework for data access and ArcGIS JavaScript API for map viewer.


I want to split polyline object in equal distance for example if I have a polyline with 60 km I want to split it in 10 polyline with 6 km distance and send polylines on ArcGIS JavaScript map api.


I can use T-Sql and C# for splitting polyline. How can I do it?




qgis - Getting polygon shapefile node coordinates and point order



I am trying to get the coordinates of each node of a polygon in QGIS. What I basically did is use the following tools:


Vector -> Geometry tools -> Extract nodes


Vector -> Geometry tools -> Export Geometry column


However only around the first 300 points had successfully obtained coordinate values in the new xcoord and ycoord fields. The other three thousand plus points had null values. What can I do to get the other points' coordinates? Also is there an easy way to get the point order of each node?




Tuesday 28 January 2020

Points layer distance from the start of line layer in QGIS


I tried to find out and read multiple similar questions but not able to get answer to my question. I have two layers; points (various location data) and lines (roads). My aim is to calculate all points distance from the start of the road layer and populate this in points layer attribute (snapshot attached). Points are little away from the road line. I checked LRS but I did not find there. I am using QGIS 2.18.


Points distance from start


points distance of multiple roads




gdal - gdalwarp leaves horizontal artifacts regridding from EASE-Grid (laea) to Polarstero (stere)


I'm trying to understand how to use gdalwarp to warp an image. I believe I'm doing the basics correctly, but I might be missing some gdalwarp options.


The basic problem is that I see horizontal artifiacts in my output image.


Here's basic steps to reproduce.


Start with a simple 722x722 image. Here's one on imagur (Direct link to png)



Now use gdal_translate to apply metadata to make this a geotiff (this is a EASE-grid Lambert Azimuthal Equal Area projection on a 1924 authallic sphere.)


gdal_translate -a_srs '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs' \
-a_ullr -4524688.262500000 4524688.262500000 4524688.262500000 -4524688.262500000 \
-co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 \
-of GTiff ./UiMbqSd.png ./UiMbqSd.withmetadata.tif

This seems to have the expected geographic information (~12.5km grid) confirmed with gdalinfo


 gdalinfo UiMBqSd.withmetadata.tif



Driver: GTiff/GeoTIFF
Files: UiMbqSd.withmetadata.tif
Size is 722, 722
Coordinate System is:
PROJCS["unnamed",
GEOGCS["unnamed ellipse",
DATUM["unknown",
SPHEROID["unnamed",6371228,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],

PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["latitude_of_center",90],
PARAMETER["longitude_of_center",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-4524688.262500000186265,4524688.262500000186265)
Pixel Size = (12533.762500000000728,-12533.762500000000728)
Metadata:

AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=DEFLATE
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-4524688.263, 4524688.263) (135d 0' 0.00"W, 29d42'45.71"N)
Lower Left (-4524688.263,-4524688.263) ( 45d 0' 0.00"W, 29d42'45.71"N)
Upper Right ( 4524688.263, 4524688.263) (135d 0' 0.00"E, 29d42'45.71"N)
Lower Right ( 4524688.263,-4524688.263) ( 45d 0' 0.00"E, 29d42'45.71"N)
Center ( 0.0000000, 0.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)

Band 1 Block=722x2 Type=Byte, ColorInterp=Red
Mask Flags: PER_DATASET ALPHA
Band 2 Block=722x2 Type=Byte, ColorInterp=Green
Mask Flags: PER_DATASET ALPHA
Band 3 Block=722x2 Type=Byte, ColorInterp=Blue
Mask Flags: PER_DATASET ALPHA
Band 4 Block=722x2 Type=Byte, ColorInterp=Alpha

Next I regrid it to a polarstero projection.


gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \

./UiMbqSd.withmetadata.tif ./regridded_lon0_-45.tif

The problem I'm seeing is horizontal artifacts that appear to protrude around the middle of the image. Here's the png representation of the problem (http://imgur.com/Nrx4ZoS)


I thought it might be something weird about my corners, but I see these artifacts with different lat_0 values.


gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \
./UiMbqSd.withmetadata.tif ./regridded_lon0_0.tif

Here's an example with no rotation (http://imgur.com/uiiF9Ir)


I'm currently running this on a mac:


> gdalwarp --version

GDAL 1.11.1, released 2014/09/24

But I've tested and seen the same behavior on UBUNTU 12.04:


GDAL 1.10.1, released 2013/08/26

Answer



Try adding the -et (error threshold) option with lower thresholds than the default (0.125). When I use "-et 0.01", the horizontal artifacts disappear:


gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \
-et 0.01 \
./UiMbqSd.withmetadata.tif ./regridded_lon0_-45.tif

Monday 27 January 2020

qgis - How to create lines to visualize differences between polygon features in PostGIS?


I've a PostGIS table polygon_b with some polygon features. There is also a table polygon_a which contains the same polygons as polygon_b but with minor changes. Now I want to create lines to visualize the differences between the polygon features.


enter image description here enter image description here enter image description here


I suppose that ST_ExteriorRing and ST_Difference will do the job but the WHERE clause seems to be quite tricky.



CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, yourSRID) AS geom
FROM
(SELECT
(ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
FROM polygon_a, polygon_b
WHERE
-- ?
) AS g;


Can anyone help me?


EDIT 1


As posted by 'tilt' I've tried ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom) but the result is not as expected.


CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, your_SRID) AS geom
FROM
(SELECT
(ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom

FROM polygon_a, polygon_b
WHERE
ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom))
AS g;

enter image description here


EDIT 2


workupload.com/file/J0WBvRBb (example dataset)




I've tried to turn the polygons into multilines before using ST_Difference, but the results are still strange.



CREATE VIEW multiline_a AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_a.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_a;

CREATE VIEW multiline_b AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_b.geom))::geometry(multilinestring, 4326) AS geom
FROM

polygon_b;

CREATE VIEW line_difference AS SELECT
row_number() over() as gid,
g.geom
FROM
(SELECT
(ST_Dump(COALESCE(ST_Difference(multiline_a.geom, multiline_b.geom)))).geom::geometry(linestring, 4326) AS geom
FROM
multiline_a, multiline_b)

As g;

enter image description here



Answer



Here are a few new tricks, using:



  • EXCEPT to remove geometries from either table that are the same, so we can focus only on geometries that are unique to each table (A_only and B_only).

  • ST_Snap to get exact noding for overlay operators.

  • Use the ST_SymDifference overlay operator to find the symmetric difference between the two geometry sets to show the differences. Update: ST_Difference shows the same result for this example. You can try either function to see what they get.



This should get what you expect:


-- CREATE OR REPLACE VIEW polygon_SymDifference AS
SELECT row_number() OVER () rn, *
FROM (
SELECT (ST_Dump(ST_SymDifference(ST_Snap(A, B, tol), ST_Snap(B, A, tol)))).*
FROM (
SELECT ST_Union(DISTINCT A_only.geom) A, ST_Union(DISTINCT B_only.geom) B, 1e-5 tol
FROM (
SELECT ST_Boundary(geom) geom FROM polygon_a
EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b

) A_only,
(
SELECT ST_Boundary(geom) geom FROM polygon_b
EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_a
) B_only
) s
) s;

rn | geom
----+-------------------------------------------------------------------------------------

1 | LINESTRING(206.234028204842 -92.0360704110685,219.846021625456 -92.5340701703592)
2 | LINESTRING(18.556700448873 -36.4496098325257,44.44438533894 -40.5104231486146)
3 | LINESTRING(-131.974995802602 -38.6145334122719,-114.067738329597 -39.0215165366584)
(3 rows)

three lines




To unpack this answer a bit more, the first step with ST_Boundary gets the boundary of each polygon, rather than just the exterior. For instance, if there were holes, these would be traced by the boundary.


The EXCEPT clause is used to remove geometries from A that are part of B, and rows from B that are part of A. This reduces the number of rows that are part of A only, and part of B only. E.g., to get A_only:


SELECT ST_Boundary(geom) geom FROM polygon_a

EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b

Here are the 6 rows of A_only, and 3 rows of B_only: A_only B_only


Next, ST_Union(DISTINCT A_only.geom) is used to combine the linework into a single geometry, typically a MultiLineString.


ST_Snap is used to snap nodes from one geometry to another. For instance ST_Snap(A, B, tol) will take the A geometry, and add more nodes from the B geometry, or move them to the B geometry, if they are within tol distance. There are probably several ways to use these functions, but the idea is to get coordinates from each geometry that are exact to each other. So the two geometries after snapping look like this:


A snapped B snapped


And to show differences, you can choose to use either ST_SymDifference or ST_Difference. They both show the same result for this example.


How to get the geometry type of an empty PostGIS table?


I have a web application that needs to know the geometry type of a PostGIS table's geom field before inserting rows into it. I use the following query to determine the geometry type of the geom column:


SELECT GeometryType(geom) FROM my_schema.building LIMIT 1


This returns the geometry type of an actual row, so it does not work when my table is empty. How do I determine the geometry type of the geometry column itself?



Answer



The query could be run against the geometry_columns table in this way


SELECT type 
FROM geometry_columns
WHERE f_table_schema = 'my_schema'
AND f_table_name = 'building'
and f_geometry_column = 'geom';

(or, if you are using a geography type, subsititute geometry_columns with geography_columns and f_geometry_column with f_geography_column)



Spatial clustering with PostGIS including attributes



Spatial clustering with PostGIS including attributes


I have a points table with some attributes like : year, type, etc. And I'd like to do a points clustering in PostGIS. I followed the subject in Spatial clustering with PostGIS


But I have the following problems:


first problem without considering the attributes: I got a table with points and I can visualise it in Qgis but not in ArcMap and the result has no SRID (the original table has srid :4326) and I'd like the output with geom as point or polygon (not geometry collection); second problem: I'd like to do a clustering for every year (2014,2015, 2016) and every type (a, b, c) and get the result in one table but I do not know if it is possible and how to do it?


CREATE TABLE table_clusterd AS (
SELECT row_number() over () AS id,
ST_NumGeometries(gc) as radius,
gc AS geom_collection,
ST_Centroid(gc) AS centroid,
ST_MinimumBoundingCircle(gc) AS circle,

FROM (
SELECT unnest(ST_ClusterWithin(geom, 100)) gc
FROM myTable
)
f);

Answer



To answer your main question, to cluster by attributes, you simply use GROUP BY attribute. AS ST_ClusterWithin will operate on any geometry type, it returns a GeometryCollection. To only pull out points, use ST_CollectionExtract(unnest(ST_ClusterWithin(geom, dist)),1). Finally, you can call ST_SetSRID(geom, 4326) to convert back to 4326. Here is an example that clusters on attributes, a, b, c and years, 2014, 2015 and 2016.


WITH 
testvals (att, year, geom) AS
(SELECT

('[0:2]={a,b,c}'::text[])[trunc(random()*3)],
('[0:2]={2014, 2015, 2016}'::int[])[trunc(random()*3)],
ST_MakePoint(random()*100, random()*100)
FROM generate_series(1, 1000)),
clusters(att, year, geom) AS
(SELECT
att,
year,
ST_CollectionExtract(unnest(ST_ClusterWithin(geom, 50)), 1)
FROM testvals

GROUP BY att, year
ORDER BY att, year)
SELECT
row_number() over() AS id,
att,
year,
ST_Numgeometries(geom) as num_geoms,
ST_AsText(
ST_SetSRID(
ST_Centroid(

ST_MinimumBoundingCircle(geom))
, 4326)
)
FROM clusters;

To convert this in to a create table, simply remove the ST_AsText (that is for debugging) and the WITH clause and replace testvals with your own table name and SELECT from clusters.


This produces something like the following:


id | att | year | num_geoms | geom
----+-----+------+------------------+-----------------------------


1 | a | 2014 | 104 | POINT(47.8713526488699 47.5714471330078)



2 | a | 2015 | 133 | POINT(50.5456980698809 48.7293153893785)


3 | a | 2016 | 92 | POINT(51.6473020685052 48.6573579298721)


4 | b | 2014 | 95 | POINT(53.6228174459529 45.5578126914714)


etc. As I have used 1000 points spread over a 100 x 100 grid with the distance parameter of 50, there are 9 clusters produced, as you would expect, for each attribute, year combo. For real world data, you will likely get more clusters. Note, if there is only one point in the cluster, then ST_Centroid(ST_MinimumBoundingCircle(... will return EMTPY POINT, which you might need to handle.


This crazy looking construct,


('[0:2]={a,b,c}'::text[])[trunc(random()*3)]

is a quick way to randomly generate a series of a, b, and c. For more see this answer and generate_series is used like a loop to create x number of samples.


Depending on the distance parameter to ST_ClusterWithin, and the number of rows created by generate_series, the above demo will return anything from 9 rows (distance 100), ie, all combinations of a, b, c and 2014, 2015, 2016 to 100 rows (distance 1). Usage may vary.


spatialite - Is it possible to add a new line in a non-spatial table using Qgis?



While using QGIS, and having loaded a non spatial table (from a dbf or a Spatialite table), is it possible to add a new row to its table attributes?


Normally in a spatial layer you have to use "add features button", that obviouslly won't work in this case.


My idea was to use Qgis as somekind of frontend to a Spatialite database, with the ability to edit\add features in spatial tables as well in the non-spatial ones.



Answer



In 1.8 there should be an add button in the attribute table


enter image description here


arcgis desktop - Creating Drop Down List Variable with ModelBuilder?


I am building a model and am trying to create a variable (that will be turned into a parameter and then added as a field) that is a drop down list where only specific values can be chosen.


Any help in how to do this. I have:




  1. right clicked the white space in ModelBuilder

  2. clicked on create variable

  3. selected string

  4. right clicked on the oval and chose 'list of values'

  5. Entered my list items


This doesn't seem to act as a drop down when running the model. Am I doing something wrong? These are text values I'd like in the drop down.



Answer



Often a particular tool that will be accepting this variable may have the option of a value list filter for the parameter.



python - Understanding UnicodeDecodeError: 'utf8' codec in ArcPy script?



I work with arcmap 10.3 and python 2.7.8. I have more than 500 shapefiles that located in many folders and subFolders. All Sub Folders are located in one large directory. I try with arcpy to detect all shapefiles that have in their attribute table ,in field name "YEUD", the value 20. I search all shape files that begin with letters "mig". Finally i tried to print all the shapefiles that had been found with value 20 in it. When i run this code:


import arcpy,os,fnmatch,unicodedata,codecs

rootPath = r"C:\Project\layers"
pattern = 'mig*.shp'
for root, dirs, files in os.walk(rootPath):
for filename in fnmatch.filter(files, pattern):
shp = os.path.join(root, filename)
if arcpy.ListFields(shp, "YEUD"):
print("{} has YEUD field".format(shp))

with arcpy.da.SearchCursor(shp, ["YEUD"]) as rows:
for row in rows:
if row[0] == 52:
print("{} has a record with YEUD = wanted row".format(shp))
break

i get an error when the python meet files and folders with right to left font:


UnicodeDecodeError: 'utf8' codec can't decode byte 0xe7 in position 23: invalid continuation byte

For completeness, i asked this question in https://geonet.esri.com/message/519769#519769 and marked it as correct answer for files and folders names that written in left to right fonts ,but when i run this code i get an error when the python meet files and folders names with right to left fonts.



In GeoNet i didn't receive helpful answer. I also searched answers in stackOverflow but didn't understand how to unicode the script.




gdal - Does gdal_calc only support 26 input raster files at a time?


In gdal_calc.py, input files are specified with letters A-Z.


Here are two files, input.tif and input2.tif, being averaged into result.tif:


gdal_calc.py -A input.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"


I have more than 26 files. What can I do about this?


One possibility: keep adding up files, 25 at a time, and then divide the values by the total number of input files. I think my numbers will become too large.


It looks like the source code for gdal_calc.py could be adapted to process as many images as needed, but the script is a mammoth nested looping function called doit(), so that's a bit scary. Should I really be hacking source code, or is there a better way?


My images are all single-band GeoTIFFs, by the way.


EDIT: After posting, I realized this question is somewhat of a duplicate of this rather incoherent one. Fortunately, the other question does have a good answer showing how to average more than 26 images together from the command line. The answer only works for a predetermined number of images, though. It would be nice to have a general solution that supports a theoretically unlimited number of images.




KML in QGis with additional Data




i´ve developed a android app that extracts the exif-information of a photo and creates a KML-file (Placemark) - some additional information like camera-type, focal length,... are stored in the tag "ExtendedData". When i try to import the KML-file in qgis the attribute table only shows the name and description field.


is it possible to import a the KML-file in qgis without loosing any information stored in "ExtendedData"?


The KML-File is like this:






20130404_132050


Orientierung
1



Brennweite
279/100



15.448561,47.11721,0





Attributes in Qgis: enter image description here


Hope you can help me! thanks, michael


EDIT: ok works fine now - i did not had the newest version of QGIS (2.0.1), the KML-File gets importet now correctly with all the attributes in "ExtendedData"


enter image description here




Sunday 26 January 2020

arcpy - Remove duplicates from a field


I need to remove duplicates from a field named "Intersecti"


enter image description here


Here is my code. I used list..set. It's not removing any duplicates.


import arcpy
duplicates = "G:\\xStreetNew\\Duplicates.shp"
duplicates_List = []

with arcpy.da.SearchCursor(duplicates, ['Intersecti']) as Cur:
for row in Cur:
duplicates_List.insert(0,row[0])
duplicates_List = list(set(duplicates_List))
print duplicates_List

Answer



You need to:



  • Use an update cursor

  • Split the string into a list


  • Remove duplicates

  • Put the list back into a string

  • Assign the string to that row/field

  • Apply the update


Try this:


import arcpy
duplicates = "G:\\xStreetNew\\Duplicates.shp"

with arcpy.da.UpdateCursor(duplicates, ['Intersecti']) as cur:

for row in cur:
row[0] = ','.join(list(set(row[0].split(','))))
cur.updateRow(row)

Or--for clarity-- you could replace the jumbled, multi-operation line, like this:


import arcpy
duplicates = "G:\\xStreetNew\\Duplicates.shp"

with arcpy.da.UpdateCursor(duplicates, ['Intersecti']) as cur:
for row in cur:

v = row[0]
v = v.split(',')
v = list(set(v))
v = ','.join(v)
row[0] = v
cur.updateRow(row)

duplication - Creating duplicate features - many to one - in qgis


Essentially I have a shapefile of polylines and wish to create duplicates and based on attribute data.


I have route data for individuals as polylines. Then I have attribute data which tells me the number of people and what are the members of the group. The number of individuals is variable ranging from 1 -10 and I would want to duplicate to create a single line for each individual.


This is to examine footfall. I realise that with a grid I could simply sum the number of members for each polyline togive me footfall. However I wish to examine different members of the group differently and treat their impacts differently, therefore I need single polylines for each member of the group.


A similar question has been asked (below), but this is in ArcGIS and I have only QGIS available.


Creating duplicate features based on many-to-one conversion of related table?





arcobjects - Looking for a region growing algorithm


I work with a raster set representing land values and I would like to automatically create polygons with random shape but with certain size and/or land value given a centroid point of each polygon. The centroid point will act as the initial raster cell (a pixel called as a seed cell) that will belong in that certain polygon and then a region growing algorithm will begin adding pixels around a centroid pixel until the constraint of size and/or land value of a polygon will be satisfied.


Is anybody who has any idea about an existing region growing algorithm? or at least which is the class of ArcObjects or Interfaces that I may work with, so as to read the value of each pixel, create a new polygon by adding new pixel to it?


I hope that the above make sense and somebody may help me.


Edit



Basically I intend to automatically create polygons with regular shapes and a series of constraints such as size, land value etc. using a genetic algorithm (GA). But first I need to create a set of random solutions in order to feed the GA.So, for instance on land block which is an area enclosed by roads I have already saying 4 centroid points. Each centroid point represents the approximate location of each new land parcel. In addition, each centroid is associated with attributes of each parcel i.e. size and land value. So, I want to begin by each centroid as a seed point of the region growing algorithm to starting creating a random shape for each parcel based on each centroid.


I hope that the above make sense. I look forward how can I create this region growing algorithm or if there is already one in VBA and ArcObjects.




How can I import georeferenced photos (jpg, kml, kmz) using QGIS or ArcGIS for Desktop?


I've got several georeferenced digital photos as jpg, kml, and kmz. How can I import them into my GIS projects (QGIS or ArcGIS for Desktop) so I can open them as a popup?




Saturday 25 January 2020

Reprojecting coordinates without building geometry object in ArcPy?


Is it possible with arcpy to reproject single x,y coordinates without building a complete geometry object?


For instance, something like this:


input = [[-8081441,5685214], [-8081446,5685216], [-8081442,5685219], [-8081440,5685211], [-8081441,5685214]]
output = [arcpy.A_Reproj_Fn((coord[0], coord[1]), source_SR, dest_SR) for coord in input]


I have huge arrays of coordinates to reproject, but for performance issue I do not want to build a geometry object before reprojecting.



Answer



There are some options available. Here is a very simple benchmarking. I am on a decent i7 machine with 16GB RAM and SSD disk, Windows 10, Python 2.7 (the built-in Python in ArcGIS 10.4.1).


Ideally you would be using pyproj from Anaconda (it is easy to connect Anaconda with arcpy, check this guide). Excellent in terms of performance: 8secs for transforming 1mln XY pairs.



import time
from functools import wraps
import numpy as np
from pyproj import Proj, transform


src_sr = Proj(init='epsg:4326')
dest_sr = Proj(init='epsg:3857')

def timethis(func):
'''
Decorator that reports the execution time.
'''
@wraps(func)
def wrapper(*args, **kwargs):
start = time.time()

result = func(*args, **kwargs)
end = time.time()
print(func.__name__, round(end-start))
return result
return wrapper

#----------------------------------------------------------------------
@timethis
def prepare_data(size):
"""create numpy array with coords and fields"""

ys = np.random.uniform(0,89,size).tolist()
xs = np.random.uniform(-179,179,size).tolist()
return zip(xs,ys)

@timethis
def pyproj(data):
res = [list(transform(inProj,outProj,pair[0],pair[1])) for pair in data]
return res

data = prepare_data(1000000)

print data[0:2]
print pyproj(data)[0:2]

The output:


('prepare_data', 0.0)
[(-68.23467501281871, 15.374738820652137), (3.4997809786621303, 22.93714083606868)]
('pyproj', 8.0)
[[-7595849.276871487, 1732425.641861154], [389593.8364326526, 2624418.652993309]]



Next in terms of performance is using arcpy.da.NumPyArrayToFeatureClass and then back with arcpy.da.FeatureClassToNumPyArray. The worst one is arcpy.PointGeometry:


import arcpy
import numpy as np
import time
from functools import wraps

def timethis(func):
'''
Decorator that reports the execution time.
'''

@wraps(func)
def wrapper(*args, **kwargs):
start = time.time()
result = func(*args, **kwargs)
end = time.time()
print(func.__name__, round(end-start))
return result
return wrapper

temp_fc = r'in_memory\temp_fc'


src_sr = arcpy.SpatialReference(4326)
dest_sr = arcpy.SpatialReference(3857)

#----------------------------------------------------------------------
@timethis
def prepare_data(size):
"""create numpy array with coords and fields"""
ys = np.random.uniform(0,89,size).tolist()
xs = np.random.uniform(-179,179,size).tolist()

if arcpy.Exists(temp_fc):
arcpy.Delete_management(temp_fc)
return zip(xs,ys)

#----------------------------------------------------------------------
@timethis
def reproject_numpy(data):
"""reproject coords with arcpy.da module"""
arr = np.array(data,np.dtype([('x', np.float64), ('y', np.float64)]))
arcpy.da.NumPyArrayToFeatureClass(arr,temp_fc,['x', 'y'],src_sr)

res = arcpy.da.FeatureClassToNumPyArray(temp_fc,'SHAPE@XY',spatial_reference=dest_sr)
return [list(i[0]) for i in res]

#----------------------------------------------------------------------
@timethis
def reproject_geometry_objects(data):
"""reproject coords with arcpy.PointGeometry.projectAs()"""
geometries = (arcpy.PointGeometry(arcpy.Point(pair[0],pair[1]),src_sr).projectAs(dest_sr) for pair in data)
res = [[geom.firstPoint.X,geom.firstPoint.Y] for geom in geometries]
return res


sizes = [10**i for i in xrange(8)]
for size in sizes:
print
print '{size}'.format(size=size).center(50,'-')
data = prepare_data(size)
print data[0:2]
print reproject_numpy(data)[0:2]
print reproject_geometry_objects(data)[0:2]


The output:


------------------------1-------------------------
('prepare_data', 2.0)
[(-167.4990262605815, 20.72648057680592)]
('reproject_numpy', 0.0)
[[-18645906.311743677, 2359294.0419395217]]
('reproject_geometry_objects', 0.0)
[[-18645906.311697092, 2359294.0419164128]]

------------------------10------------------------

('prepare_data', 2.0)
[(-84.10710438738408, 54.32017524422917), (149.53056786201995, 60.5826412111139)]
('reproject_numpy', 0.0)
[[-9362760.0324575398, 7231028.3662902983], [16645666.672426889, 8530614.7980484683]]
('reproject_geometry_objects', 0.0)
[[-9362760.032500302, 7231028.3663340295], [16645666.672429098, 8530614.79807427]]

-----------------------100------------------------
('prepare_data', 0.0)
[(99.45852992001386, 10.777413690256289), (-128.66525647722594, 22.3855564491687)]

('reproject_numpy', 0.0)
[[11071672.905741969, 1206874.3036907075], [-14322950.83380558, 2557879.2814767426]]
('reproject_geometry_objects', 0.0)
[[11071672.905743506, 1206874.3037197446], [-14322950.833830737, 2557879.281497049]]

-----------------------1000-----------------------
('prepare_data', 0.0)
[(-173.27374656367525, 8.223262860596597), (-139.61519355591957, 10.62062833548439)]
('reproject_numpy', 0.0)
[[-19288745.235347211, 918568.44701982441], [-15541892.253658244, 1189112.2559826041]]

('reproject_geometry_objects', 1.0)
[[-19288745.235311065, 918568.4469744475], [-15541892.253649296, 1189112.25603746]]

----------------------10000-----------------------
('prepare_data', 0.0)
[(11.194744968264729, 7.877971732593098), (54.58669010557381, 7.112696412558614)]
('reproject_numpy', 0.0)
[[1246193.3093983289, 879748.16972375475], [6076562.5466901502, 793823.26923564379]]
('reproject_geometry_objects', 9.0)
[[1246193.309427791, 879748.1696780229], [6076562.546642701, 793823.2691861242]]


----------------------100000----------------------
('prepare_data', 0.0)
[(-119.9078743699582, 9.317778132275732), (-46.418166430278006, 4.464482461184021)]
('reproject_numpy', 1.0)
[[-13348083.516972212, 1041852.8393190451], [-5167246.6505450243, 497487.58635374776]]
('reproject_geometry_objects', 90.0)
[[-13348083.516967565, 1041852.8393501436], [-5167246.650575973, 497487.5863742894]]

---------------------1000000----------------------

('prepare_data', 0.0)
[(-69.46181556994199, 5.285715860826253), (-150.64679833442418, 17.05875285677861)]
('reproject_numpy', 7.0)
[[-7732453.938828676, 589239.59320861078], [-16769924.88017785, 1927665.2915218703]]

---------------------10000000---------------------
('prepare_data', 5.0)
[(110.4244111629609, 30.50859019906087), (-89.59380469024654, 30.291453068724223)]
('reproject_numpy', 84.0)
[[12292389.221812241, 3569093.3287834972], [-9973536.7163228001, 3541068.701018624]]


1mln XY pairs can be transformed in ~7secs with the numpy array conversion method. 10mln XY pairs can be transformed in ~85secs with the numpy array conversion method.


Not too bad! Pretty much the same as using pyproj. Unless you can't use pyproj, stick to this one.


qgis - "gdalbuildvrt command not found" in Build Virtual Raster tool


I am trying to use the "Build Virtual Raster" tool in QGIS 3 to build a mosaic. When I run the tool, I get the following message:


GDAL command output:

/bin/sh: gdalbuildvrt: command not found

gdalbuildvrt is available in my terminal, but it appears that QGIS does not know where to find it. How do I fix this?


I am running:



  • QGIS 3.0.3

  • macOS 10.13.4



Answer



This is a know issue, however there is a workaround described in the readme file from the installer:



QGIS 3 does not find external tools needed for Processing. Some configuration options are missing. A workaround is to use QGIS custom variables. Open QGIS Preferences → System.



  • Under the Environment section, turn on Use Custom Variables.

  • Add (+) a variable

  • Select Prepend from the Apply popup

  • Enter PATH for the variable

  • For the Value enter: /Library/Frameworks/Python.framework/Versions/3.6/bin:/Library/Frameworks/GDAL.framework/Versions/2.2/Programs:

  • Quit QGIS 3 and start it again


carto - Sunlight layer into CartoDB


I have a list of tweets on a torque map made with CartoDB. I want to add an additional layer that shows the sunlight around the world at certain times, similar to the Twitter #sunrise map found here.


Does anyone know how to add a layer like this?



Answer



It is a Leaflet plugin you can implement in any CartoDB map, you can find it here: https://github.com/joergdietrich/Leaflet.Terminator



python - how to overlay shapefile and raster?


I have a shapefile with polygons. And I have a global raster file. I want to overlay the shapefile's polygons onto the raster grid and calculate the mean raster value for each polygon.


How can I do this using GDAL, writing the results to the shapefile?



Answer



Following advice I got on the gdal-dev mailing list, I used StarSpan:


starspan --vector V --raster R1 R2 ... --stats mystats.csv avg mode


The results are saved to CSV format. At that time, that was already enough for me, but it should be possible to somehow forge a Shapefile from that info.


dem - QGIS: Edit elevation in DTM to get correct slope of a river


I am new to QGIS and just need it to do a hydrolocic analysis of an area of which I later make a flow model with HEC-RAS 2D.


I made it according to this tutorial: https://www.youtube.com/watch?v=Ro-RRzMMw-c&t=1774s


Unfortunately there is a railway with lower terrain nearly at the end of my river which is misunderstood as the stream so I can't make the catchment area to the very end of the river.


How can I edit the elevation of the DTM to get the correct stream (or is there another way)?


My river "Kraimser Bach" with the wrong end of the stream. It would go along the railway to north direction and not into the bigger river:


wrong catchment are



I tried to chagne the z-value of the DTM along this shapefiles to get the correct stream but it didnt worked the way I tried it. DTM with shapefiles i would like to edit the z-value in the DTM


I tried to do it with rasterize from polygons but that didn't work. (According to this: Modifying DEM elevation data in QGIS using a polygon) I also tried it with the SAGA tool "Burn streams into dem" but therefore I would need the real slope and so it didn't worked as well.




postgis - Spatial joining syntax



A standard PostGIS point-in-polygon query might look like this. (example tutorial)


SELECT p.*,b.polyname
FROM points p, shapes s
WHERE ST_Within(p.the_geom,s.the_geom)

Imagine this was not a spatial join. Then it would look like this:


SELECT i.*,c.categoryname
FROM items i, categories c
WHERE i.zipcode = c.zipcode


Most SQL folks would advise rewriting that query using a variety of the JOIN keyword syntax:


SELECT i.*,c.categoryname
FROM items i
JOIN categories c ON i.zipcode = c.zipcode

Does PostGIS offer this sort of syntax for spatial joins? I can't find it. All I see are a lot of functions that return true or false, leading to the workaround syntax. Ideally, I would think a spatial join would look like:


SELECT p.*,b.polyname
FROM points p
JOIN shapes s
ON

p.the_geom IS WITHIN s.the_geom

Answer



To use [FULL|INNER|LEFT|RIGHT|OUTER|CROSS] JOIN syntax supported by PostgreSQL, the query would still need ST_Within(). I've never seen any spatial database use IS WITHIN, since there could be parameters for the spatial relation, such as ST_DWithin. It should look like this:


SELECT p.*
FROM points p
INNER JOIN shapes s ON ST_Within(p.the_geom, s.the_geom)

remote sensing - What is NIR LiDAR?


Can LiDAR sensors collect near-infrared (NIR) data at the same time or does there need to be a separate NIR camera that is combined with the LiDAR data after collection?


Blue Marble makes reference to it in their latest release of global mapper: http://www.directionsmag.com/pressreleases/updated-global-mapper-lidar-module-now-supports-nir-lidar-data/444480



Answer



The simple answer is that lidar sensors coupled with NIR cameras can collect point cloud data that can then have the NIR values "embedded" with them, the same way RGB values can be assigned to point cloud data collected with high res photos.


Friday 24 January 2020

Performing gdal translate to png/jpeg?



I tried gdal translate to png:


gdal_translate image.bsq image.jpg -of JPEG -outsize 10% 10% -scale

enter image description here


and jpeg:


gdal_translate image.bsq image.png -of PNG -outsize 10% 10% -scale

enter image description here


jpeg looks better and file size is smaller. In general should I expect poor result with png?


Using the -r cubic and bilinear are almost the same to my eyes:



enter image description here




Why can PostGIS create LINESTRING from the same values of coordinates of the POINT type?


So, I needed to create linear objects from point objects.





  1. I split the objects into points with the following query:

    CREATE TABLE roads_arc_var3 AS SELECT (ST_DumpPoints(geom)).geom FROM roads_arc_var1;




  2. I created the following table, so that it created two identical fields with the same values of coordinates:


    CREATE TABLE roads_arc_var4 AS SELECT a.geom AS a, b.geom AS b FROM roads_arc_var3 AS a, roads_arc_var3 AS b WHERE a.geom = b.geom




  3. I added the field "geoline" to the table:


    SELECT AddGeometryColumn('roads_arc_var4', 'geoline', 4326, 'LINESTRING', 2);





  4. I changed the values of the initial and final coordinates in some records in order to create lines from them, and I expected that in the 'geoline' field line values will appear only in those records in which the values of the coordinates in the starting point are not equal to the values of the coordinates end point, the rest will remain empty


    UPDATE roads_arc_var4 SET geoline = ST_MakeLine(a, b);



  5. But, the picture was different, see figure


enter image description here


What is the explanation for the PostGIS developers of this situation, because it's about 2D?


Or it is a problem with rounding off the values of coordinates, but in this case the user needs to understand what is its length, etc.




Answer



Many of the PostGIS functions may return topologically invalid geometries. It is often more practical to do so than to throw an error or stop the process.


SELECT ST_AsText(
ST_MakeLine(
ST_GeomFromText('POINT (1 1)'),
ST_GeomFromText('POINT (1 1)')));

Result: "LINESTRING(1 1,1 1)"


SELECT ST_IsValid(
ST_MakeLine(

ST_GeomFromText('POINT (1 1)'),
ST_GeomFromText('POINT (1 1)')));

Result: "f"


SELECT ST_IsValidReason(
ST_MakeLine(
ST_GeomFromText('POINT (1 1)'),
ST_GeomFromText('POINT (1 1)')));

Result: "Too few points in geometry component[1 1]"



Fortunately it is not hard to add some logic into your SQL. You can for example add WHERE NOT ST_Equals(a.geom,b.geom) or then check the result with ST_IsValid and save only the valid values.


Arcpy mapping text element reset


I am writing a script to create thousands of maps and I am trying to create a table that updates for every map. Basically an updating table of the information specific to each map. I have been trying to link a an excel table to the layout view but arc does not seem to have that functionality. So what I am doing instead is creating a manual table in the layout view. I have 27 values for each map that I have been storing in a dictionary (using searchcursors to extract the info and store into dictionary).


I cannot copy and paste my table into here but it essentially looks this:


Receptors               AEGbuffer    CONbuffer    ERGbuffer
buffer 1 10 19
schools 2 11 20
childcares 3 12 21
hospitals 4 13 22
nursinghomes 5 14 23
critical 6 15 24

sw 7 16 25
streams 8 17 26
respop 9 18 27



for k,v in sorted(table.items()):
for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
if elm.text == "1":
if k == 'AEGbuffer': elm.text= v
if elm.text == "2":

if k == 'AEGschools': elm.text = v
if elm.text == "3":
if k == 'AEGchildcares': elm.text = v
if elm.text == "4":
if k == 'AEGhospitals': elm.text = v
if elm.text == "5":
if k == 'AEGnursinghomes': elm.text = v
if elm.text == "6":
if k == 'AEGcritical': elm.text = v
if elm.text == "7":

if k == 'AEGsw': elm.text = v
if elm.text == "8":
if k == 'AEGstreams': elm.text = v
if elm.text == "9":
if k == 'AEGrespop': elm.text = v

this is just for the first row, the rest of the code is the same process for the pasting of the dictionary values.


where I am stuck: I do not know how to properly reset the newly changed values to their original 1-27 values so the next map can be easily produced. is there a tool or trick to this?


what I have tried: after I export the map I reopen the text elements and reverse the process like..


for k,v in sorted(table.items()):

for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
if k == 'AEGbuffer':
if elm.text == v: elm.text = "1"
and so on......

but this approach hasn't worked great. a few of the values are wrong when there are the same values a few times on the table.



Answer



Instead of setting the text contents of each element set their name to something like "Text1" or even better you could make the element names match up with the key names in your table like 'AEGbuffer'. If it were me I would also add a suffix to the fields so i can make sure I am operating only on the table elements. You might come back later and add text elements that you don't want to include in your table loop.


enter image description here


Then your giant loop could look like this



for table_elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "*_table"):
if text_elm.name == k:
text_elm.text = v

To reset you just


for table_elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "*_table"):
table_elm.text = ""

EDIT: Based on your comments I think your real problem is your data format.


{'hospitals':{'AEGBuffer': 'None','ConBuffer': '65,200', 'ERGBuffer': 'None'}, 'schools':{'AEGBuffer': '600','ConBuffer': '7,300', 'ERGBuffer': '550'},'respop':{}}



lets you do:


#your k is now going to be the row name, i.e. hospitals and v is going to be another dictionary
for k,v in sorted(table.items()):
#this relies on elements being named hospitalERGBuffer, schoolCONBuffer etc.
for bufferType, bufferValue in v.items():
fieldName = str(k + bufferType) # i.e. 'hospitalsERGBuffer'
table_elm = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", fieldName + "*")[0] #there's only going to be one if you've done this right.
table_elm.text = bufferValue

Is there a 3rd party ArcGIS viewer (builder) for Javascript?


Esri provides an ArcGIS Viewer for Flex as well as an ArcGIS Viewer for Silverlight.



Strangely though, there are no plans for a similar viewer for Javascript.


Are there any 3rd parties that have created a "viewer" for javascript?


Though they are referred to as "viewers", they are actually viewer-builders: they allow administrators to create web applications without writing any code. The javascript viewer blogged about here does not allow user to choose tools to include in the app, as can be done with the flex and silverlight versions.




Update


I have developed some capabilities that have been packaged as a combination of an SOE and an Add-in for the silverlight viewer. This works great - when the admin builds a viewer application, he can pick my add-in tool from the list and configure it to use the URL of a server where my SOE is installed. However, from what I can tell there are no add-ins for the javascript viewer. Or am I missing something?



Answer



Here is ESRI's Web App Builder http://doc.arcgis.com/en/web-appbuilder/


It does require an AGS account.


Thursday 23 January 2020

Getting distance between two points over linestring using PostGIS?


I have a linestring: LINESTRING(-1.3326397 50.9174932,-1.3319842 50.9166939)


and I select two points that lay on that linestring.


How do I get the distance between these two points (in meters)?




python - Splitting line every x meters using QGIS?




I have a line that I ultimately want to split up into points. The points should be every 100 meter along the line. So I don't want to extract the nodes.


Are there any Open Source (QGIS, Python) tools around for that?


The use case is that I have a bus line without bus stops. Though I know every 100 meter the bus stops. This way I want to generate bus stops to use as a GTFS feed.



Answer



You can use the Convert lines to points tool (you need SAGA GIS installed and the processing toolbox plugin enabled) and set your distance:


Convert lines to points


This is what I received for my line layer:


Line and points


I used the Measure Line tool from the toolbar to do a quick check between points:


Measured line



Hope this helps!


installation - Multiple versions of GDAL (32bit and 64bit) on same machine


I currently have GDAL installed on my Windows 7 64bit workstation. It is the 32bit version. I chose the 32 bit version because I wanted to use the GDAL Python bindings with a 32bit installation of Python (version 2.6 installed with ArcGIS).


Would it be possible to also install the 64bit version of GDAL on the same computer so that I can take advantage of more memory and processors when running the GDAL command line utilities?




Answer



It shouldn't be a problem to have several Python/GDAL installations on a single Windows computer.


For example, let me count how many GDAL installations are on my 64-bit work PC:



  1. ArcGIS 10.0 with bundled Python 2.6 (32-bit). I only use this Python installation for Esri things

  2. Python 2.7 (64-bit) with GDAL, etc. Good for using lots of RAM (and I really do use it)

  3. Python 3.2 (64-bit) with GDAL, etc. Mostly for future-proof testing of stuff I make

  4. OSGeo4w has GDAL, although I needed to modify the Osgeo4W.bat file to make it work

  5. Some groundwater software that I use includes GDAL

  6. Some geological modelling tool I've been demoing includes GDAL



Occasionally, there could be a clash between different installations, but there is usually a way to fix it.


qgis - How to avoid the CRS Dialog when creating a vector layer memory?


I am trying to create a vector memory in a plugin. Basing on some examples,I was able to create a memory vector layer and to add it in QGis with : QgsMapLayer.instance().addMapLayer(myLayer,True)


But when the layer is added , QGis always displays the crs dialog to specify the crs of the layers


:crs dialog


I tried to specify myself the crs with : layer.setCrs(QgsCoordinateReference(4326)) but it still shows me the crs dialog.


Do you have any ideas to how avoid this dialog ?




arcobjects - How to create a transparent buffer of a polyline?



I use vb.net arcobjects create a buffer around a polyline feature once it is selected. The problem is the buffer has a obvious colour. Since I have other functions on the selected feature, I cannot clear the buffer because user may still need the feature selected.


When I select another feature, the new buffer will created and old buffer cleared. But there always a buffer which makes the polyline there looks different from other features in the same layer.


So I am thinking whether can make the buffer transparent so that it wouldn't affect the overview of the map.


My code is from an ESRI web sample:


 Public Sub CreateGraphicBuffersAroundSelectedFeatures(ByVal activeView As ESRI.ArcGIS.Carto.IActiveView, ByVal distance As System.Double)

'parameter check
If activeView Is Nothing OrElse distance < 0 Then
Return
End If


Dim map As ESRI.ArcGIS.Carto.IMap = activeView.FocusMap
' Clear any previous buffers from the screen
Dim graphicsContainer As ESRI.ArcGIS.Carto.IGraphicsContainer = CType(map, ESRI.ArcGIS.Carto.IGraphicsContainer) ' Explicit Cast
graphicsContainer.DeleteAllElements()

' Verify there is a feature(s) selected
If map.SelectionCount = 0 Then
Return
End If


' Reset to the first selected feature
Dim enumFeature As ESRI.ArcGIS.Geodatabase.IEnumFeature = CType(map.FeatureSelection, ESRI.ArcGIS.Geodatabase.IEnumFeature) ' Explicit Cast
enumFeature.Reset()
Dim feature As ESRI.ArcGIS.Geodatabase.IFeature = enumFeature.Next()

' Buffer all the selected features by the buffer distance and create a new polygon element from each result
Dim topologicalOperator As ESRI.ArcGIS.Geometry.ITopologicalOperator
Dim element As ESRI.ArcGIS.Carto.IElement
Do While Not (feature Is Nothing)

topologicalOperator = CType(feature.Shape, ESRI.ArcGIS.Geometry.ITopologicalOperator) ' Explicit Cast
element = New ESRI.ArcGIS.Carto.PolygonElementClass()
element.Geometry = topologicalOperator.Buffer(distance)
graphicsContainer.AddElement(element, 0)
feature = enumFeature.Next()
Loop

activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGraphics, Nothing, Nothing)

End Sub


I am using vb.net arcobjects and arcmap 10.2. Any advice on it?



Answer



The following VBA code grabs the first polygon graphic on the Map and sets the transparency. The problem is that graphics don't appear to support true transparency. In the Help for IColor.Transparency is states:



...For graphic elements, 0 for transparent and 255 for opaque are the only supported values...



Public Sub MakeTransparent()
' Get map document
Dim pMXDoc As IMxDocument

Set pMXDoc = ThisDocument

' Get Map
Dim pMap As IMap
Set pMap = pMXDoc.FocusMap

' Get Graphics Container for map
Dim pGC As IGraphicsContainer
Set pGC = pMap
pGC.Reset


' Get first graphic, assumed to be polygon
Dim pElement As IElement
Set pElement = pGC.Next

' Get symbol of graphic and its colour
Dim pFillShapeElement As IFillShapeElement
Set pFillShapeElement = pElement
Dim pFillSymbol As IFillSymbol
Set pFillSymbol = pFillShapeElement.Symbol

Dim pColour As IColor
Set pColour = pFillSymbol.Color

' Set transparency
pColour.Transparency = 0 ' Only 0 or 255 is support
pFillSymbol.Color = pColour
pFillShapeElement.Symbol = pFillSymbol

' Refresh map
pMXDoc.ActiveView.Refresh

End Sub

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...