Monday 29 February 2016

Generating Contours from csv Concentrations


I have a shp file showing locations of various wells on a site. I have linked (joined) a csv file of analytical results to this file. I am now looking for a method to generate contours based upon the analytical results. Idea been to identify by contours areas of similar analytical concentrations. Any ideas?



Answer



First you will want to generate a raster surface. To do that you want to start with points with values. You should be able to create a point file with the values from the analytical results that you have (assuming you have coordinate values in your table). Then the size of the raster surface that you create will depend on the spatial density of points that you have. It is best if you have at least one point (better with a few points) for every cell in the raster.


To create the surface you must use some method of interpolation. You can use Kriging and get into the whole art of creating a surface or use an easier method like IDW (I recommend IDW). This is available with Spatial Analyst (in that toolbox). If you are using QGIS see Spatial Analysis (Interpolation) and Inverse Distance Weighted.


Next typically to make smooth contours, you want to smooth your raster surface. NOTE: This degrades the data, so you might want to try using unsmoothed data to see how it looks.


Optional Smoothing: Use FocalStatistics and a Mean (I've used a 9x9 neighborhood to good effect). Focal Statistics is in the Spatial Analyst toolbox and from what I understand there is a method for doing this with GRASS see ... "Is there a Focal Statistics (spatial analyst) tool in ArcGIS equivalent in QGIS?"


Next whether you smoothed the raster or not, you are ready to generate contours. In ArcMap the Contour Tool is again under spatial analyst tools. but you can generate contours at the OSGEO4W shell with the command gdal_contour. Here is a sample command that I've used.



gdal_contour -a contour inputraster output.shp -i 10

And that will create contours in intervals of 10 (most commonly feet or meters but in your case it will be in intervals of whatever units your analytical results are measured in).


arcgis desktop - Batch export all layers in ArcMap document that have selected features?


I have several layers in an ArcMap Document. All of these layers are stored as shape files.



Specific features of the said layers are selected. I want to export the data as shapefile programatically. How can I do this using arcpy?


Sample view of my arcmap workspace



Answer



Try..


    ## This script will look for all the layers that have feature selected in them in the TOC and export them in seperated shapefile.
##output layer name will be the original name+ _selected e.g. luzon_loss_selected
import arcpy
arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames( mxd, "Layers")[0]

layers = arcpy.mapping.ListLayers(mxd, "*", df)
output_folder = r'C:\Users\YOUR_USER_NAME\Desktop\test' ##folder path of output
for layer in layers:
Sel=arcpy.Describe(layer)
if Sel.FIDset:
arcpy.FeatureClassToFeatureClass_conversion(layer,output_folder,layer.name + "_selected","","","")

arcgis desktop - Unable to rotate Marker Line Symbols


I have a line feature class that I am trying to symbolize with a Marker Line Symbol type. I have my symbology set up but I am having a problem with the rotation angle.



When I enter an angle in the box and click ok, the symbol does not rotate. If I go back into the Symbol Selector dialog box, the angle is still shown as 0. In this example, I am trying to enter 180 in the Angle box to make the arrow turn 180 degrees. Is this a bug, or is there a workaround?


enter image description here



Answer



I think it's a bug, or more perhaps more worded more precisely: it's a result of inadequate user interface portrayal of the symbol properties and their effects. My deduction is that the same dialog box for point symbol properties is re-used, but not all the controls there are applicable for Simple Line Symbols.


It would help to also add an additional graphic showing what the desired final complete line will look like, even a pencil sketch. For example, are you really after a marker line, where every x distance has a symbol, or is it just arrows at the middle and/or ends?


For a true symbol marker line, you need to use Cartographic Representation, which requires Desktop Advanced license. If you're fortunate enough to have that, see how to get started with fix symbol rotation to page.


Rich man's (CartoRep) example: equally spaced fixed rotation marker symbol.


[example: fixed marker along outline


Poor man's repeated marker line symbol: it works, angle is fixed at 45deg to page according to symbol properties, but the distance between each point varies.


example: simple repeated marker with rotation



click path for simple repeated marker with rotation


To round out the examples, for arrows at ends like you'd use to depict stream flow, see Line Properties >> Line Decorations and rotate to follow line. Here we've created a complex symbol, two arrow types on one line, by stacking lines (see "Layers" panel).


Stream flow with different arrows at beginning and end


qgis - PyQGIS make screenshot of mapCanvas after setExtent is called


I have the following function:


def createScreenshots(layer):
counter = 0
for feature in layer.getFeatures():
counter += 1

point = feature.geometry().asPoint()
qgis.utils.iface.mapCanvas().setCenter(point)
qgis.utils.iface.mapCanvas().refreshAllLayers()
qgis.utils.iface.mapCanvas().saveAsImage("D:\\\\m\\testing\\" + str(counter) + ".png")

Screenshots are not in real coordinates, .pgw header is not in right coordinates, therefore resulting png images are shifted.


Question: How to wait for rendering of canvas after setCenter is called?



Answer



ids = None


def createScreenshot4(layer):
global ids
ids = layer.allFeatureIds()

def exportMap():
global ids
qgis.utils.iface.mapCanvas().saveAsImage( "D:\\\\ma\\bo\\{}.png".format( ids.pop() ) )
if ids:
setNextFeatureExtent()
else: # We're done

qgis.utils.iface.mapCanvas().mapCanvasRefreshed.disconnect( exportMap )

def setNextFeatureExtent():
reqq = QgsFeatureRequest()
reqq.setFilterFid(ids[-1])
for feature in layer.getFeatures(reqq):
point = feature.geometry().asPoint()
qgis.utils.iface.mapCanvas().setCenter(point)
qgis.utils.iface.mapCanvas().refreshAllLayers()



qgis.utils.iface.mapCanvas().mapCanvasRefreshed.connect( exportMap )
setNextFeatureExtent() # Let's start

Resulted code with setCenter and refreshAllLayers.


coordinate system - Assigning CRS to shapefile when it doesn't have one, in R?


I think that by using readOGR to import a shapefile renders worthless to put proj4string to the script since if any coords are available in the shapefile, it will be transported through the readOGR.


Therefore, proj4string is expedient only if we use another import function for example readShapeSpatial (for polygons, lines or points) or by using readOGR, but in the case that the shapefile does not contain the coordinates/projection:


When checking if the shapefile has projection information and if it has not, should I add proj4string although I have put readOGR?



Answer



As JeffreyEvans states, readOGR from the rgdal library imports a CRS if there is one embedded in the shapefile. You can check by (example using a shp I've been playing with):


proj4string(india)  # from the sp package
# [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

If this returns NA you can specify a CRS with:



proj4string(india) <- CRS("+init=epsg:24383")

You can obtain the EPSG CRS strings from http://spatialreference.org/


To answer your specific question, yes, you should define a projection if one is not specified.




If one is not certain which CRS to assign/define, refer to the following post:



And if one wants to reproject a spatial object in R, i.e., change its projection, then refer to:



Google Earth Engine: Export an entire collection


I would like to export all the images that the collection (code) contains. How can I do that?


var colection = ee.ImageCollection ("MODIS/NTSG/MOD16A2/105")


Answer



If you want an automatized process you'll have to use the Python API. If you can do that, follow:



  1. Install geetools (https://github.com/gee-community/gee_tools)



pip install geetools




  1. Export Collection using a python script. As you want to export a MODIS collection, I highly recommend you use a region parameter.



from geetools import batch
import ee
ee.Initialize()

col = ee.ImageCollection("MODIS/NTSG/MOD16A2/105")
region = ee.Geometry.Polygon([XXXX])

help(batch.ImageCollection.toDrive) # See help for function


tasks = batch.ImageCollection.toDrive(col, 'MyFolder', region=region, scale=1000)

Otherwise, if you don't mind clicking "Run" for every image, you can use geetools in the code editor (https://github.com/fitoprincipe/geetools-code-editor), follow:


var tools = require('users/fitoprincipe/geetools:batch')
var col = ee.ImageCollection("MODIS/NTSG/MOD16A2/105")
var region = ee.Geometry.Polygon([XXXX])

batch.ImageCollection.toDrive(col, 'MyFolder',
{region: region,
scale: 1000})

r - Cluster geometries that touch each other


I'm trying to cluster the geometries in a sf object that touch each other. For the moment I'm using the same approach that it's explained in the answer here, but the problem is that the same idea doesn't work if I consider a bigger number of geometries because the touching matrix becomes too memory-heavy and I'm not using its sparse properties.


A small and (I hope) reproducible example


# packages
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright

# download data
milan_primary_highways <- opq("Milan, Italy") %>%
add_osm_feature(key = "highway", value = "primary") %>%
osmdata_sf() %>%
trim_osmdata(bb_poly = getbb("Milan, Italy", format_out = "polygon")[[1]]) %>%
magrittr::use_series(osm_lines)


# Determine the non-sparse matrix of touching geometries
touching_matrix <- st_touches(milan_primary_highways, sparse = FALSE)
#> although coordinates are longitude/latitude, st_touches assumes that they are planar

# Cluster the geometries
hc <- hclust(as.dist(!touching_matrix), method="single")

# Cut the dendrogram
groups <- cutree(hc, h=0.5)


# Results
table(groups)
#> groups
#> 1 2 3 4 5 6 7 8 9 10 11 12
#> 1444 21 18 8 4 5 8 2 3 5 2 2
plot(st_geometry(milan_primary_highways), col = groups)


Is there an alternative approach that takes advantage of the sparse property of the touching matrix (that I cannot use because as.dist just accepts a numeric matrix). Do you want me to provide an example when this approach does not work? The error is simply "R cannot allocate a vector of size n GB" when I use the function hclust(as.dist(x)).


I don't know if it's important but I tried to explain the reason why I need this clustering here




Answer



Make an adjacency list:


touching_list = st_touches(milan_primary_highways)

this isn't NxN so won't grow - its a list of length N where each element is only the adjacent elements. This might help it scale up to larger problems.


Now use igraph to create a graph objects and get the connectivity:


library(igraph)
g = graph.adjlist(touching_list)
c = components(g)


The membership element of c is of length N and says which cluster it is in. It summarises the same as your example:


> table(c$membership)

1 2 3 4 5 6 7 8 9 10 11 12
1444 21 18 8 4 5 8 2 3 5 2 2

So if you want to assign it to the data, do:


 milan_primary_highways$groups = c$membership

arcgis pro - For-Looping through Bookmarks in Python Script tool?


I'm using this script tool, in ModelBuilder for ArcGIS Pro, to set up a map layout, reference a bookmark, and export the layout as a PDF. This model iterates through Urgent Care centers and its corresponding data one at a time (in alphabetical order by center name). Since my model loops through each Urgent Care center (via the Iterate Field Values), it repeats the same bookmark through every loop/iteration. How do I tell it to go from one bookmark to the other?? (I have about 90 centers, so that means 90 bookmarks also)


Here is my script tool:


import arcpy

#input layer

lyr = arcpy.GetParameterAsText(0)

# input name of layout
p = arcpy.mp.ArcGISProject("CURRENT")
lyt = p.listLayouts("Layout_King")[0]

# Reposition the scale bar
scaleBar = lyt.listElements("MAPSURROUND_ELEMENT", "Scale Bar")[0]
mf = scaleBar.mapFrame
scaleBar.elementPositionX = mf.elementPositionX + 0.0

scaleBar.elementPositionY = mf.elementPositionY - 0.5

# Reposition the north arrow
northArrow = lyt.listElements("MAPSURROUND_ELEMENT", "North Arrow")[0]
mf = northArrow.mapFrame
northArrow.elementPositionX = mf.elementPositionX + 8.8
northArrow.elementPositionY = mf.elementPositionY + 0.7

# Align the title with the center of the map frame
title = lyt.listElements("TEXT_ELEMENT","Name of Map Text")[0]

mf = lyt.listElements('MAPFRAME_ELEMENT',"Map Frame")[0]
title.elementPositionX = mf.elementPositionX + (mf.elementWidth / 3.7)
title.elementPositionY = mf.elementPositionY + (mf.elementHeight / 0.98)

# Reposition the Legend and fix legend title
legend = lyt.listElements("LEGEND_ELEMENT", "Legend")[0]
legend.title = "Legend"
legend.elementPositionX = mf.elementPositionX + 7.7
legend.elementPositionY = mf.elementPositionY + 7.15


# setting layout to bookmark
aprx = arcpy.mp.ArcGISProject("Current")

# add name of layout
lyt = aprx.listLayouts("Layout_King")[0]
mf = lyt.listElements("MAPFRAME_ELEMENT")[0]

# add name of bookmark
bkmks = mf.map.listBookmarks()
bkmks.sort(key=lambda x: x.name, reverse=True)

for bkmk in bkmks:
mf.zoomToBookmark(bkmk)
lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\" + bkmk.name + ".pdf")

I have heard that this could be accomplished via a for-loop, but don't know how to write one.




Sunday 28 February 2016

arcgis desktop - Changing/adding keyboard shortcuts in ArcMap?


How can I change or add keyboard shortcuts in ArcMap?


I have tried Customize --> Customize mode --> Commands --> Keyboard. However, this does not seem to actually work. I can add shortcuts and link them, but nothing happens when I use them even though they seems to be saved (my created shortcut still shows up after restarting Arcmap).


I'm not using English keyboard layout so, for example, I can't use the [Ctrl] + [Shift] + [=] for Zoom to selected features, since I need Shift to even get a =. It would be really useful to be able to change or add alternative shortcuts for the tools I use regularly.


Does anyone know how? Preferably something within Arcmap, but if it comes down to writing an AutoHotkey script that'll work too if you can give me some hints.


There are a few questions relating to shortcuts already, but none that I found is giving me a useful answer. I have also looked at all the lists of shortcuts that's been implemented, both on ArcGIS Resources Online and this very useful pdf. No hints on how to change the shortcuts though.


COMMENT: If the shortcut includes [ALT GR] it doesn't seem to work properly.



Answer



You need to go to Customize> Customize Mode on the Menu bar, in ArcMap.


Once the Customize dialog appears, click on the keyboard button.



For more details, have a look at this help article: Assigning a shortcut key


Comparison of Open Source Desktop GIS Packages


I'm curious as to what comparisons are out there between different Open Source Desktop GIS software packages. I'm aware that there has been a QGIS and gvSIG comparison on this site but I am looking for comparisons between any OpenSource Desktop GIS.


A substantial list of Open Source GIS software can be found in response to the question: What are some Free and Open Source GIS Desktop packages


I don't expect a single comparison with everything in it, but comparative subsets. Also any comprehensive benchmark/testing that may have been performed between any grouping of these products.



Answer




From what I've seen, this is the most exhaustive comparison matrix of GIS software out there: https://docs.google.com/spreadsheet/ccc?key=0Albk_XRkhVkzdGxyYk8tNEZvLUp1UTUzTFN5bjlLX2c&hl=en#gid=0


There are quite a lot of variables in this matrix, some which may be out of scope for your current study, but the issue of application scalability can be drawn from such a topic.


Other links http://en.wikipedia.org/wiki/List_of_geographic_information_systems_software http://en.wikipedia.org/wiki/Comparison_of_geographic_information_systems_software


Leaflet draw: prevent drawing through a custom overlay


I have a leaflet draw implementation (in Angular) with a custom sidebar that overlaps the Map. It should not be possible to draw through that sidebar. For circlemarkers I could solve this problem easily by using


pointer-events: all;


and


event.stopPropagation();

for the clickevents. That way any clicks will not arrive on the map, so no circlemarker is drawn.


My problem: this solution does not work when drawing eg. rectangles. Even with the code above, any click on the sidebar starts drawing a new rectangle or finishes the drawing of a rectangle.


How can I prevent the Map from registering these events on top of the sidebar. Why does it work for circlemarkers an not for polygons?




qgis - Querying through Relations


Right now I am stumped on how to get two layers to interact. I have one layer that is just a map with town borders, with attributes for town codes. My other layer are data for businesses, wages for example. I setup a relate in the project properties menu between them using town codes, and I can see in the attribute table the the map layer knows what businesses exist in each town now.


However, I can't seem to make use of it. If, for example want to select by expression to just select towns where the average wage is over a certain level, I can't even begin because I can't write an expression using values from both layers, even though I set them up in a relation. How would I really go about doing this?


I am trying to figure out how to query, select, or anything using attributes from two related layers.



Answer



I'd suggest a two step process:




  1. If your businesses do not have the town codes associated with them, then you need to do a spatial join. In QGIS this is under Vector - Data Management Tools - Join Attributes by Location. The target layer will be your businesses and the join layer will be the towns (I'm assuming the towns are a polygon layer). Every business (assuming these are points) will be assigned the attributes of the town where it is located. Once this is done, you'll be able to query the new business layer using the attributes of the towns - and create summary tables (next step).





  2. If you want to join summarized data from the businesses back to the towns (so you can map by town or do town-based queries) you'll need to download a QGIS plugin that will allow you to sum by attributes to create new tables. Under Plugins - Manage and Install Plugins, take a look at Group Stats or Statist. You can use one of these to summarize the businesses by the town attributes you assigned in step 1, to create average wages for the town (for example). Once you have business data summarized by town, then you can do a regular one to one attribute table join, to join your summary table back to the towns layer to map the summarized data.




wms - GeoServer: Receive multiple tiles from layer


Using GeoServer you can do a WMS request to receive a single PNG representation of a certain layer using something like this.


http://localhost:8080/geoserver/wms?service=WMS&version=1.1.0&request=GetMap&layers=spearfish&styles=&bbox=589425.9342365642,4913959.224611808,609518.6719560538,4928082.949945881&width=512&height=359&srs=EPSG:26713&format=image%2Fpng

Is it possible to alter this request to directly request different "zoom levels" or "tiles" of the layer? I assume you need to do a Web Map Tile Service (WMTS) request but I can't seem to find any documentation on how to do this. Any guidance would be helpful.



Answer




If you want to get a 256*256 PNG tile of a certain map service from geoserver you can do it with including tiled=TRUE in your request:


http://localhost:8080/geoserver/wms?LAYERS=spearfish&STYLES=&FORMAT=image%2Fpng&TILED=true&TILESORIGIN=589425.93423656%2C4913959.2246118&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A26713&BBOX=589425.93423656,4913959.2246118,609518.67195605,4934051.9623313&WIDTH=256&HEIGHT=256

But there's no zoom parameter in the request, and the client side application will take care of the certain zoom level and the bounding box.


But if you insist on using WMTS (which capability is included in GeoServer with GeoWebCache integration) you can do it with the following request:


http://v2.suite.opengeo.org/geoserver/gwc/service/wmts/?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=medford%3Abuildings&STYLE=_null&TILEMATRIXSET=EPSG%3A900913&TILEMATRIX=EPSG%3A900913%3A13&TILEROW=3030&TILECOL=1300&FORMAT=image%2Fpng 

Here the TILEROW and TILECOL are the key parameters. But if you want to extend your knowledge about WMTS I would suggest you to check out the official OGC standard here:


http://www.opengeospatial.org/standards/wmts


And you can find useful information about WMTS in geoserver at the following link:



http://docs.geoserver.org/latest/en/user/webadmin/tilecache/defaults.html?highlight=wmts


I hope this helps!


Saturday 27 February 2016

gdal - Installing rasterio and Fiona together?


I work on Windows; I use Christoph Gohlke wheels; Python 3.7


After installing GDAL-3.0.1 and Fiona-1.8.6


an error occurs on import of fiona from fiona.ogrext import Iterator, ItemsIterator, KeysIterator ImportError: DLL load failed:


by downgrading to GDAL-2.4.1 I could surpass the error, yet rasterio starts rasing an error : from rasterio._base import gdal_version ImportError: DLL load failed:




fields attributes - Is there are way to import non-geo JSON files into QGIS?


Here is a json file containing various attributes relating to countries. However, it does not contain any geographic data.


http://bost.ocks.org/mike/nations/nations.json


I would like to join these attributes with another file that has the geographic data for the countries. But I can't figure out how to get it into QGIS.


Is there a way of loading non-geojson json files into QGIS?



Answer



Yes, you can import non-geo JSON files into QGIS but there are no native method for this in QGIS. It's because in your link, it's non tabular data. For each "column" (like income) per country, you have year and value. You need to preprocess the data to change the structure to flatten it to 2 dimensions.


One choice could be to use PyQGIS console to solve this issue. In first, you get the data with standard Python urllib2 and json libraries and process it with Python. Then, using a memory provider, you can integrate the data as a layer and make your join. Be aware that if you need to reuse the data later, you will need to execute again the operations. So, you will need to save your memory layer or use the "Memory Layer Saver" plugin (but this plugin keeps data persistent per project).



Another choice could be to avoid QGIS completely at the beginning and generate a CSV from the data using Python (outside QGIS). Then, you can load the CSV file in QGIS.


For PyQGIS, look at the PYQGIS cookbook.


Where can I find historic country border data?


I am looking for border polygons from World War I (from 1914-1918). Is there a resource out there that maintains this kind of information?




Creating square grid polygon shapefile with Python?


I have the following coordinates


minx, maxx, miny ,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073


I wish to create a square grid of size 1 m using python.


import math


minx,maxx,miny,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073

size = 1

def set_bbox(minx, maxx, miny, maxy, distx, disty):
nx = int(math.ceil(abs(maxx - minx)/distx))
ny = int(math.ceil(abs(maxy - miny)/disty))
new_maxx = minx + (nx*distx)
new_miny = maxy - (ny*disty)
return ((minx, new_maxx, new_miny, maxy),ny,nx)

# shift the bottom (right - down)

coord, ny, nx = set_bbox(minx,maxx,miny,maxy,size,size)
# left-up origin
origin = coord[0],coord[3]
# number of tiles
ncell = ny*nx


arcgis desktop - Solution to "Error 000918: Cannot retrieve feature class extent"?



I have an ArcGIS model where I merge two multipoint feature classes, Then I would like to use "Point to Raster" on it. Cell extent is specified through the input DEM which all data derives from. It worked many times, but now after changing the storage location of an intermediate result right before the "Point to Raster" tool, I get the error message "Error 000918: Cannot retrieve feature class extent". What does that mean? I didn't change the input data, but the model doesn't work anymore.


System: Win7, ArcGIS for Desktop Advanced 10.2.1, 3D and Spatial Analyst.



edit: When I run the model to this point (Zusammenführen/Merge (2)) as shown below, the problem does not occur anymore. Why is it not possible to run the whole model at once?


See the specified part of the model below: enter image description here




coordinate system - Is there projection issue with Natural Resources Canada (NRC) Canada Atlas hydrology data?


I have been using the Natural Resources Canada (NRC) data for various generalized maps but it is not until now when i overlaid higher resolution features obtained from the provincial government (Ministry of Natural Resources Land Information Ontario) that I notice a shift in the the NRC Hydrology (watercourse, water bodies) or approximately 800m.


What is strange is that all other data from the same source such as transportation, places, etc. match perfectly with when overlaid with the higher resolution provincial data.


The Provincial data is in NAD83 UTM16 and the Atlas Canada data in WGS84.


ArcMap is projecting on the fly which is usually not an issue (other than performance).


alt text



Answer



GSRUG http://www.geod.nrcan.gc.ca/tools-outils/gsrug_e.php


It is different with every province in Canada http://www.fmepedia.com/index.php/Canadian_Spatial_Reference_System_%28CSRS%29_High_Precision_Datum



lots of info here: http://www.geod.nrcan.gc.ca/faq_e.php#faq1


Friday 26 February 2016

qgis - How to change the edit widget for fields in a layer with PyQGIS?


Is possible change the edit widget with PyQGIS?


I have made an extension that loads a PostGIS layer and I want to edit the date fields with the date widget by default without any change in the settings of layer ( "layer properties dialog").



Answer



You can use the method setEditorWidgetV2() from QgsVectorLayer.


To set the edit widget of a field to DateTime, do this:


vLayer.setEditorWidgetV2( fieldIndex, 'DateTime' )


Parameters can be set via a dictionary of key-values, like this:


vLayer.setEditorWidgetV2Config( fieldIndex, {'display_format': 'yyyy-MM-dd', 'allow_null': False, 'field_format': 'yyyy-MM-dd', 'calendar_popup': True})

enter image description here


You can find a description of the available options at the QgsDateTimeEditWrapper documentation.


A list of available widgets can be found here.




On the other hand, if you want all your Date fields to have the DateTime edit widget, you could do something like this:


fields = vLayer.dataProvider().fields()

dateFields = []
# Get indices of date fields
for field in fields:
if field.typeName() == 'Date':
dateFields.append( vLayer.fieldNameIndex(field.name()) )

# Set the DateTime editor widget for date fields
for dateField in dateFields:
vLayer.setEditorWidgetV2( dateField, 'DateTime' )

Following symbolic links in QGIS to load data on Linux?



Can QGIS 2.4 follow Linux soft links?


I want to add a number of raster layers from a store of NED30m elevation data, and thought I'd do it by creating a collection of softlinks in one directory to the elevation files nested in the store.


However, when I try to add a layer from one of the links, QGIS 2.4 complains:



Raster layer: somepath/elev_ned30m_links/n46w105.adf is not a supported raster data source



However, n46w105.adf is the result of a linux ln -s command:



ln -s /Mapstore/ned30m/grid/n46w105_1/w001001.adf




QGIS will open w001001.adf just fine if approached directly, but that means digging down into each tile's directory structure, while the softlinks can be created programmatically with a script.



Answer



Yes QGIS can follow symbolic links.


However, an ArcGIS Grid coverage can't be linked like that, the raster data is not just the .adf file, it is a directory of related files stored in a parent directory (workspace) with an INFO directory and you can not change either the directory structure or the .adf filename (both of which you did with the ln -s command), see format description:


structure


However, you have another option...


Use your script to create VRT headers pointing to the Grid .adf files and store them in a single directory. Replace the bit in your script that does the ln -s ... with


gdal_translate -of VRT /absolute/path/to/wNNNN.adf /absolute/path/to/vrt/directory/someraster.vrt

Alternatively, if you wish to mosaic the GRIDS, create a textfile with the paths of all the GRIDs (either /path/to/GRID folder path or /path/to/GRID/w001001.adf file path) and then use



gdalbuildvrt -input_file_list /path/to/gridlist.txt your_new_mosaic.vrt

postgis - Create closest point on multilinestring to use in shortest_path()


I am working a routing based application where I need to create a route between two points of interest. The problem which I am facing occurs after calculating network topology using assign_vertex_id(). When I run shortest_path I get the path between vertexes and not from the point of interest (as shown in image).


enter image description here


Can you please help me out with this: how can I use this query or any other function to do it?



Thanks for any tips.


EDIT 1


I am using this function called multiline_locate_point() similar to line_locate_point() which helps me to find nearest LINESTRING from MULTILINESTRING to the point.


here it is:


-- Function: multiline_locate_point(geometry, geometry)

-- DROP FUNCTION multiline_locate_point(geometry, geometry);

CREATE OR REPLACE FUNCTION multiline_locate_point(amultils geometry, apoint geometry)
RETURNS geometry AS

$BODY$
DECLARE
mindistance float8;
nearestlinestring geometry;
nearestpoint geometry;
i integer;

BEGIN
mindistance := (distance(apoint,amultils)+100);
FOR i IN 1 .. NumGeometries(amultils) LOOP

if distance(apoint,GeometryN(amultils,i)) < mindistance THEN
mindistance:=distance(apoint,GeometryN(amultils,i));
nearestlinestring:=GeometryN(amultils,i);
END IF;
END LOOP;

nearestpoint:=line_interpolate_point(nearestlinestring,line_locate_point(nearestlinestring,apoint));
RETURN nearestpoint;
END;
$BODY$

LANGUAGE plpgsql IMMUTABLE STRICT
COST 100;
ALTER FUNCTION multiline_locate_point(geometry, geometry) OWNER TO postgres;

this gives me the_geom() of point on line. Now the problem is shortest_path() only takes vertex (integer number) as input to create route. So how can I pass this point to shortest_path() so that it creates route from the point and not from any vertex ??




Are polygons stored clockwise or counterclockwise in a shapefile?



I'm working with polygon coordinates in QGIS using PyQGIS using ESRI shapefiles as input data.


Are the polygon vertices (I am referring only to the outer boundary in the case of polygons with holes) stored as clockwise or as counterclockwise in a shapefile or can be in mixed order?


Is the order always the same or can be different depending on the source of the data: digitizing, converting from other sources as *.csv files, etc.?



Answer



According to the Shapefile white paper (page 8), polygons are stored in a clockwise fashion except for interior polygon hole parts (i.e. 'donut holes'), which are stored in a counter-clockwise fashion. In fact, that's how polygon holes are recognized.


Loading a raster using QGIS into a PostGIS2.0 enabled database



Following on from my attempt to install a postgresql 9.1 db with postgis 2.0 on windows 7, I am trying to figure out how to load a raster. I have successfully loaded a shapefile, and am trying to do the same with a .tif file using the Load Raster to PostGIS plugin (version 0.5.1) in QGIS (version 1.7). I have set up the connection to my db, and am using the following settings:


enter image description here


When I click on the OK button, I get the below error message. I've tried this with a .adf file and also a .tif file, both projected in Albert Equal Area Conic with an SRID of 102003.


Checking parameters...
Connecting to database...
Storing overview 1 on database...
Failed.
Finished storing overview 1.
Storing overview 2 on database...
Failed.

Finished storing overview 2.
Storing overview 3 on database...
Failed.
Finished storing overview 3.
Finished.

This process has not inserted anything into my database, and I don't understand the error message. Some previous related questions are here (asked by me, using a different process on loading rasters) and here (using the same plugin but much earlier in the year).



Answer



celenius,


It's quite possible that the QGIS raster loader hasn't been revised yet to fit with the new changes we have in place. I suggest trying the raster2pgsql.exe packaged with the latest windows experimental. It should load tifs just fine.



Instructions here.


http://www.postgis.org/documentation/manual-svn/using_raster.xml.html#RT_Loading_Rasters


Thursday 25 February 2016

python - Raster diff: how to check if images have identical values?


Is there a means to check to see if any 2 given raster layers have identical content?


We have a problem on our corporate shared storage volume: it's now so big that it takes over 3 days to conduct a full back up. Preliminary investigation reveals one of the biggest space consuming culprits are on/off rasters that really should be stored as 1-bit layers with CCITT compression.


a typical present/not-present raster


This sample image is currently 2bit (so 3 possible values) and saved as LZW compressed tiff, 11 MB in file system. After converting to 1bit (so 2 possible values) and applying CCITT Group 4 compression we get it down to 1.3 MB, almost a full order of magnitude of savings.



(This is actually a very well behaved citizen, there are others stored as 32bit float!)


This is fantastic news! However there are almost 7,000 images to apply this too. It would be straightforward to write a script to compress them:


for old_img in [list of images]:
convert_to_1bit_and_compress(old_img)
remove(old_img)
replace_with_new(old_img, new_img)

...but it's missing a vital test: is the newly compressed version content-identical?


  if raster_diff(old_img, new_img) == "Identical":
remove(old_img)

rename(new_img, old_img)

Is there a tool or method which can automatically (dis)prove the content of Image-A is value-identical to the content of Image-B?


I have access to ArcGIS 10.2 and QGIS, but am also open to most anything else than can obviate the need to inspect all these images manually to ensure correctness before overwriting. It would be horrible to mistakenly convert and overwrite an image that really did have more than on/off values in it. Most cost thousands of dollars to gather and generate.


a very bad result


update: The biggest offenders are 32bit floats that range up to 100,000px to a side, so ~30GB uncompressed.



Answer



Try converting your rasters to numpy arrays and then check to see if they have the same shape and elements with array_equal. If they are the same, the result should be True:


ArcGIS:


import arcpy, numpy


raster1 = r'C:\path\to\raster.tif'
raster2 = r'C:\path\to\raster.tif'

r1 = arcpy.RasterToNumPyArray(raster1)
r2 = arcpy.RasterToNumPyArray(raster2)

d = numpy.array_equal(r1,r2)

if d == False:

print "They differ"

else:
print "They are the same"

GDAL:


import numpy
from osgeo import gdal

raster1 = r'C:\path\to\raster.tif'

raster2 = r'C:\path\to\raster.tif'

ds1 = gdal.Open(raster1)
ds2 = gdal.Open(raster2)

r1 = numpy.array(ds1.ReadAsArray())
r2 = numpy.array(ds2.ReadAsArray())

d = numpy.array_equal(r1,r2)


if d == False:
print "They differ"

else:
print "They are the same"

boundless suite - OpenGeo SDK modify the getfeatureinfo popup


I would like to modify the opengeo sdk gxp.FeatureEditPopup. I want to construct a clickable hyperlink, I would like to add a new field to the popup that would be "http://someurl.com/" followed by a field from my data: "http://someurl.com/someDataField". Any ideas would be appreciated. Ben




arcgis desktop - Georeferencing CAD files in ArcMap?


I've definitely got at least a couple of grey hairs while trying to convert some dwg. files to shapefiles. I've seen that this subject has been disscused a lot, but hope someone will have the patience to help me.


The procedure I'm following is this:



  1. I add the dwg files to ArcMap.

  2. TOC- right clik-Convert CAD Feature Layer

  3. Data Management Tools-Projections- Define Projection.


The conversion happened but then I realized things are not that simple, looked over the files again and apparently I have some "lost points"....enter image description here


My data is acctually in that middle point and it looks like this: enter image description here



I definitely have to georeference the data before converting it to shapefile. But how do I do that? I tried to create some control points, but I'm probably not doing it the right way cause nothing happens. Plus I have no idea where those lost points should actually be. enter image description here


Forgot to mention that the dwg files have no coordinate system.



Answer



For future reference, it helps to add as much relevant information to your question as possible in order to obtain the best answer. In this case, something that would be useful is to say what projection you think your CAD file should be in, and possibly give a screenshot of the coordinates you are seeing when you have loaded the CAD file into ArcGIS.


You have a few options to ensure that your CAD data will line up with other data you may add. Both of these involve defining a projection for your CAD drawing, without editing the drawing itself.



  1. If you know the features in the CAD drawing are already in a projected coordinate system, you may simply define the projection for the dataset. Here is the Help document describing that process: Defining a coordinate system for a CAD dataset.

  2. If you know the CAD features are not in a projected coordinate system, but want to define a projection for the drawing anyway, thus enabling it to be projected on-the-fly in ArcGIS, then here is the procedure to Create a custom projection file in ArcMap to align CAD data.

  3. If instead of creating a custom projection for your CAD data, but instead want to georeference your CAD file, which creates a .wld file with parameters to shift, scale and rotate your drawing, you will want to use the following help file: Georeferencing a CAD dataset



The reason that I have focused on locating the CAD file with reference to your other data as opposed to converting data from CAD to a different GIS format is the time savings involved. It is a lot faster to define the projection, or georeference your CAD dataset once, knowing that any data you convert out of that drawing will then carry along that projection information, rather than leaving the CAD file as it stands, then having to define the spatial reference separately of each dataset you export from the CAD file.


To complete the process of converting to a shapefile, you would then do the following:


Use the various tools for importing CAD data to convert features from your CAD drawing to a shapefile or any other desired GIS format. This help file contains links to the various options: Importing CAD data
At the time you perform the import, you will have the option of keeping the data in the coordinate system of your CAD file, or using the coordinate system of the Data frame, or the Feature Dataset you are importing the data into.


carto - Using Torque to show rises and falls?


I'm looking at a handful of different mapping projects that I'd like to do something Torque-like with, but I'd like maps to stay on the map for a limited period of time.


Some fake data:


Address                   Start Year         End Year

7706 Heather Limits 2002 2006
5938 Noble Shadow Inlet 2002 2003
9128 Round Pioneer Link 2003 2005
7176 Crystal Via 2004 2005
4475 Harvest Square 2004 2007

Think of a map of buildings participating in a rent subsidy program. One building enters the program in 1988 and exits in 2008. A dozen enter the program in 2002, and will exit in 2032. It seems like I can add points to the map based on a single date, but there's not an easy way to take them off on a specific day, too.


I'd have to turn the "cumulative" option off and then create basically a row for each site for each year that it should be on. Something like


TorqueYear   Address                   Start Year         End Year
2002 7706 Heather Limits 2002 2006

2002 5938 Noble Shadow Inlet 2002 2003
2003 7706 Heather Limits 2002 2006
2003 5938 Noble Shadow Inlet 2002 2003
2003 9128 Round Pioneer Link 2003 2005
2004 7706 Heather Limits 2002 2006
2004 9128 Round Pioneer Link 2003 2005
2004 7176 Crystal Via 2004 2005
2004 4475 Harvest Square 2004 2007
2005 7706 Heather Limits 2002 2006
2005 9128 Round Pioneer Link 2003 2005

2005 7176 Crystal Via 2004 2005
2005 4475 Harvest Square 2004 2007
2006 7706 Heather Limits 2002 2006
2006 4475 Harvest Square 2004 2007
2007 ...

Am I right about that or am I missing something?




python - Adding GeoPandas Dataframe to PostGIS table?


I have a simple GeoPandas Dataframe:


enter image description here


I would like to upload this GeoDataframe to a PostGIS table. I have a Database setup with the PostGIS extension already but can't seem to add this Dataframe as a table.


I have tried the following:



engine = <>
meta = MetaData(engine)
eld_test = Table('eld_test', meta, Column('id', Integer, primary_key=True), Column('key_comb_drvr', Text),
Column('geometry', Geometry('Point', srid=4326)))
eld_test.create(engine)
conn = engine.connect()
conn.execute(eld_test.insert(), df.to_dict('records'))


tiles - Rendering custom OpenStreetMaps style (land=white, water=black): do I need a dedicated computer just to do this?


EDIT: My goal is NOT to use OpenStreetMap or cloudmade as a web map. I described exactly what I want ("a full, rasterized set of tiles for the entire world at zoom=10 (around 68 Gigapixels, water=black, land=white, no labels"). I only meant to use the cloudmade.com link as an example to illustrate what I'm shooting for. (As it is, the cloudmade style is only 90% of the way there.) I need the tiles offline. I need the image (68 Gpix) I described as a mask in a larger raster map task for a custom map I'm building using NASA and USGS data.




As an illustration


I made a custom style that sort of shows what I want:




(No labels, no roads, no features of any kind except all water features (rivers, lakes, oceans, you name it) as black with land being solid white, and the nice Mapnik antialiasing of values in between.)




What I really want


I would like a full, rasterized set of tiles for the entire world at zoom=10 (around 68 Gigapixels).


I considered downloading and installing Planet.osm (11 GB compressed) and running Mapnik (for days/weeks?) to get this data. However, Planet.osm looks like it will take 100GB or maybe even 1TB once the database and index are built. I don't have that kind of disk space on my laptop and since I don't need the full dataset, is there a smaller subset of the data I could download (it says here that the polygon data is only 700 MB)? Actually from what I can tell, the data has doubled in size in the last year so I would need a beefy machine to deal with this.


Is there an easier, direct way to get these tiles (level 10, black and white, just land on water) batched, or do I basically need a dedicated computer to do this?



Answer



Do you have to use OSM? or would a similar but smaller vector dataset work for you? If so consider Natural Earth (http://naturalearthdata.com) which has nice 1:10M scale coastline, land area, ocean, river and lake layers. You could then use GeoServer or MapServer locally (or on a remote server) to create your tiles at any depth you need with any of the usual tile caches (TileCache, GeoWebCache etc.)


From your question it isn't even clear to me that you need tiles - if all you want is a mask then you can probably do this using GRASS - v.toRaster() (I think).


javascript - Leaflet GeoJSON sublayers checked not working


Good morning all,


I have been struggling with division of my geoJSON sublayers, as per the query below:


Leaflet geoJSON sublayers checked on off


I used the proposed code, which seems to bring a nice solution:


    document.querySelector("input[name=vm]").addEventListener('change',  

function() { //VM is a main geoJSON layer
if(this.checked) {
if (!map.hasLayer(job2)) map.addLayer(job2);
if (!map.hasLayer(infill)) map.addLayer(infill);
if (!map.hasLayer(mdu)) map.addLayer(mdu);
if (!map.hasLayer(featurelayer2)) map.addLayer(featuresLayer2);
document.querySelector("input[name=infill]").disabled = false;
document.querySelector("input[name=mdu]").disabled = false;
document.querySelector("input[name=infill]").checked = true;
document.querySelector("input[name=mdu]").checked = true;

}
else {
if (map.hasLayer(job2)) map.removeLayer(job2);
if (map.hasLayer(infill)) map.removeLayer(infill);
if (map.hasLayer(mdu)) map.removeLayer(mdu);
if (map.hasLayer(featuresLayer2)) map.removeLayer(featuresLayer2);
document.querySelector("input[name=infill]").disabled = true;
document.querySelector("input[name=mdu]").disabled = true;
document.querySelector("input[name=infill]").checked = false;
document.querySelector("input[name=mdu]").checked = false;

}
});

document.querySelector("input[name=infill]").addEventListener('change',
function() { // Infill is 1st geoJSON sublayer, belonging to VM layer

if(this.checked) { if (!map.hasLayer(infill)) map.addLayer(infill); if (!map.hasLayer(featurelayer)) map.addLayer(featuresLayer); } else { if (map.hasLayer(infill)) map.removeLayer(infill); if (map.hasLayer(featuresLayer2)) map.removeLayer(featuresLayer2); } });


    document.querySelector("input[name=mdu]").addEventListener('change', 
function() { // Infill is 1st geoJSON sublayer, belonging to VM layer


if(this.checked) { if (!map.hasLayer(mdu)) map.addLayer(mdu); if (!map.hasLayer(featurelayer2)) map.addLayer(featuresLayer2); } else { if (map.hasLayer(mdu)) map.removeLayer(mdu); if (map.hasLayer(featuresLayer2)) map.removeLayer(featuresLayer2); } });


Unfortunately it returns not useful result. I would like to show it in the pics below: 1 - Initial map stage (layer and both sublayers visible) 2 - Main Layer switched off (works fine) 3 - Main layer switched on again (return black dots instead of orange as per to the style. Sublayers still disabled).


enter image description here


And now I considered it in other hand - start disable sublayers instead of the main layers: 1 - Initial map stage 2 - MDU layer switched off (unfortunately Infill is gone instead apart from these ones, gripped by geoJSOn feature) 3 - Both sublayers switched off (but MDU still visible), what should be equalled to whole GeoJSON layer disabled 4 - Infill layer disabled, MDU switched off again (returning black dots as an additional MDU layer - black dots comes from the style not customised yet). enter image description here


What way whould I take to solve this problem?




python - Clustering polygons by sum of values


I would like to group polygons that (1) are contiguous and (2) when the combined values of field X >= Y. Specifically I am grouping zipcodes that are connected. a group is when 1 or more zipcodes 'population' field equals 100,000.


Can someone suggest an existing ArcGIS function or python method for doing this?



Many thanks in advance!


---update----


I wrote a python algorithm to group the polygons using their respective neighbors as grouping members until threshold is met. it works nicely.




postgresql - PostGIS: incorrect interpretation of a polygon that intersects the 180th meridian


I try to find out all geopoints which intersect the polygon set as a parameter.


The problem is when I pass polygon that roughly covers an area of Bering Strait (nearby 180 longitude): enter image description here


So I use the query:


SELECT ST_AsText(l.geo_point)
FROM "lightnings" "l"
WHERE (ST_Intersects(ST_GeomFromText('Polygon((132.0 40.0, -148.0 40.0, -148.0 -8.0, 132.0 -8.0, 132.0 40.0) )', 4326), geo_point));


As you can see, vertexes are set in a correct order, clockwise, from the North-West. But the result covers outside area and including whole other world.


For example, in result:


POINT(75.5637 40.0434)

The problem doesn't touch the 0th meridian.


Elementary test:


SELECT ST_Area(ST_GeomFromText('Polygon((0.0 60.0, 10.0 60.0, 10.0 40.0, 0.0 40.0, 0.0 60.0) )', 4326))
UNION ALL
SELECT ST_Area(ST_GeomFromText('Polygon((-5.0 60.0, 5.0 60.0, 5.0 40.0, -5.0 40.0, -5.0 60.0) )', 4326))

UNION ALL
SELECT ST_Area(ST_GeomFromText('Polygon((175.0 60.0, -175.0 60.0, -175.0 40.0, 175.0 40.0, 175.0 60.0) )', 4326))

gives the result:


200
200
7000

Is there any simple trick to force PostGIS understanding me? I don't like an idea to divide the polygon...




Wednesday 24 February 2016

Import shp to Postgis using Python and ogr



I just export a Postgis Table to shp using this tips but I'm not able to import a shp to Postgis using the same library(ogr). Any idea? Thanks a lot f.



Answer



In pure Python, without using the subprocess module (os.system is deprecated) to call ogr2ogr or shp2pgsql, for example):


1) with ogr



2) with ogr and psycopg2 from the book Python Geospatial Development (Eric Westra), Chapter 7, p.219


import os.path  
import psycopg2
import osgeo.ogr
connection = psycopg2.connect("dbname=... user=...")

cursor = connection.cursor()
cursor.execute("DELETE FROM countries")
srcFile = os.path.join("DISTAL-data", "TM_WORLD_BORDERS-0.3","TM_WORLD_BORDERS-0.3.shp")
shapefile = osgeo.ogr.Open(srcFile)
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
name = feature.GetField("NAME").decode("Latin-1")
wkt = feature.GetGeometryRef().ExportToWkt()
cursor.execute("INSERT INTO countries (name,outline) " +"VALUES (%s, ST_GeometryFromText(%s, " +"4326))", (name.encode("utf8"), wkt))


connection.commit()

3) with psycopg2 only



4) with psycopg2 and other spatial libraries



qgis - How do you symbolise a polygon to look like an AutoCAD revision cloud?


I work for an engineering firm and the standard for annotating areas which require review or areas of interest is to add AutoCAD revision clouds. There is a dedicated tool for this in AutoCAD but I want to create them in QGIS.


Drawing them freehand will take too long so is there a way of symbolising a simple polygon to look like a revision cloud or alternatively is there a way of drawing a polygon which will automatically take on the shape of a revision cloud?



AutoCAD revision cloud examples:


enter image description here


enter image description here


enter image description here




qgis - Sorting points to flight path route


Original question: I would like to sort my points in QGIS. See image


Restated question: The grid tool currently numbers by starting at R1C1 and going to R1Cx, then R2C1 to R2Cx. How do I number the points from R1C1 to R1Cx, then R2Cx to R2C1? Can the solution be applied after a regular grid has been destroyed (points removed) or does it need to be applied to a whole grid, assuming no option for this ordering scheme exists in the grid creation tool itself?


The point is to avoid the long 'row reset' or 'back to start of row' line as this wastes time and fuel, rather than moving directly down to the next row.


enter image description here


enter image description here




leaflet - Client-Side Polygon split


I wanted to know if there is a way to split a polygon using JavaScript. So I have a web map which uses Leaflet API, and the user needs to cut polygons which are shown in the map using a straight line. I was hoping there would be a JavaScript library something JSTS or Turf but if that isn't possible an algorithm would do just fine.




geoserver - Labels disappear while zooming


We are using GeoServer 2.10.0. While displaying the labels in the maps, the labels which are there on the border of the tiles(4*4) are not getting displayed properly. the labels can be seen after further zooming in or out(when they fall in the center of the tile and not in the border). in this image you can see the labels are being cut


On zooming further, the labels come properly, as expected


In the first image, the labels are cut.


1.png



But in the second image, we see on zooming further, the labels are shown as expected.


2.png


Please see the SLD.


SLD:



xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
xmlns="http://www.opengis.net/sld"
xmlns:ogc="http://www.opengis.net/ogc"
xmlns:xlink="http://www.w3.org/1999/xlink"

xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

pxs_gna_mobile_site

Geo_Style_pxs_gna_mobile_site


Geo_Rule_pxs_gna_mobile_site
350000.0





image/png

40







geom




Arial

12
normal
bold




0.5
0.5



5
5




2

#FFFFFF




#000000

-1
true
yes
true
false




Geo_Rule_pxs_gna_mobile_site
350000.0




image/png


40









Invalid GeoJSON when inserting data to a CartoDB (PostGIS) table from Leaflet



I've been following along to this tutorial from CartoDB.


I create the GeoJSON using:


var drawing = "'"+JSON.stringify(layer.toGeoJSON())+"'";

and then inserting using a custom SQL function which includes:


_the_geom := ST_SetSRID(ST_GeomFromGeoJSON(_geojson),4326);

I get the error: Invalid GeoJSON representation. I've tried entering the sql my javascript produces directly into the CartoDB console and get the same error. Example below.


SELECT insert_data(
'{"type":"Feature","properties":{},

"geometry":{"type":"Point",
"coordinates":[-71.11295700073242,42.37896312335581]}}')

Reading the documentation on the function I read "ST_GeomFromGeoJSON works only for JSON Geometry fragments. It throws an error if you try to use it on a whole JSON document." Could the error I get be because of the empty properties object?



Answer



The ST_GeomFromGeoJSON really only expects to be fed the geometry portion of the object, so this works:


select st_astext(st_geomfromgeojson(
'{"type": "Point",
"coordinates": [-71.11295700073242,42.37896312335581]}'
));


So you'd want something more like this in your JS:


var drawing = "'"+JSON.stringify(layer.toGeoJSON().geometry)+"'";

pyqgis - Setting layer visibility in QGIS Python API?


Im having a problem toggling the visibility of map layers in pyqgis (2.4). I have a simple interface, just the map canvas, no legend, and basically a radio button that toggles between two raster layers.


Layers are initialized like this:


lyr = QgsRasterLayer(os.getcwd() + '/data/rasters/dof/dof.vrt', 'Ortophoto')
QgsMapLayerRegistry.instance().addMapLayer(lyr, False)
ml = QgsMapCanvasLayer(lyr)
ml.setVisible(False)

lyr2 = QgsRasterLayer(os.getcwd() + '/data/rasters/chm/chm.vrt', 'Chm')
QgsMapLayerRegistry.instance().addMapLayer(lyr2, False)

ml2 = QgsMapCanvasLayer(lyr2)
ml2.setVisible(True)

When i want to toggle their visibility i do:


ml.setVisible(True)
ml2.setVisible(False)

But unfortunately, the rendering on the map canvas stays the same.


What am i missing?



Answer




At the moment (QGIS 2.x) (and I know it's a bit gross) but you also have to update the layer set on the canvas to update the layer state change


canvas.setLayerSet([ml2, ml])

This will also control the render order. m1 will reneder first then ml2 on top of it.


raster - How to create shadows of buildings with r.sunmask?


I need to create image-raster with shadow of relief and shadow of buildings in one. I know how to create shadow of relief but how i can create shadow of buildings with r.sunmask? I dont know how from 2D polygon layer of buildings create shadow of complete buildings.



Answer




You need to run r.sun, not r.sunmask for this task. There are also convenient GRASS GIS addons to run it in hourly or daily mode.


You first need to extrude your buildings with v.extrude, then run the solar computation. See for extrusion "Extrude 2D polygons to 3D" and for a shadow example here.


enter image description here


(image courtesy: Markus Neteler)


enter image description here


(image courtesy: Vaclav Petras)


qgis - CSV with XY coordinates into a grid


I'm a beginner QGIS user and am struggling with creating a grid. I've used the MMQGIS geometry import function to import data from the CSV file shown below but they pop up as points on QGIS. I want to make it into a grid so that I can pretty much have each of the points be in a box format. I've also tried to use the join function but am not sure what variable both the grid and the points have so that I can join them.


enter image description here



Answer




In QGIS3 you can create a layer from a delimited text file using [Layer] > [Data Source Manager] > [Delimited Text] which allows you to select the delimiter (CVS is one of the options) and to select 2 of the fields (columns) to be used as the X and Y position of the point. The other fields are stored in the Attribute Table.


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


Tuesday 23 February 2016

geodjango - Certain MultiPolygons cause “BOOM! Could not generate outside point!” in PostGIS?


I've just posted this at StackOverflow and someone suggested I try here. Original: https://stackoverflow.com/q/5863859/246265



I'm trying to represent a rectangular area which crosses 180 degrees longitude. For more background see https://stackoverflow.com/q/5737506/246265


Here's my test case:


from django.contrib.gis.geos import Polygon, MultiPolygon
from my_project.my_app.models import Photo

a = Polygon.from_bbox((30, -80, 180, 80)) # the part to the east of 180
b = Polygon.from_bbox((-180, -80, -160, 80)) # a part to the west of 180
c = Polygon.from_bbox((-180, -80, -100, 80)) # a larger part to the west of 180

ok = MultiPolygon(a,b)

ok2 = MultiPolygon(c)
boom = MultiPolygon(a,c)

# This works
Photo.objects.filter(location__coveredby=ok)[:1]
# This also works so c is ok
Photo.objects.filter(location__coveredby=ok2)[:1]
# This gives "BOOM! Could not generate outside point!"
Photo.objects.filter(location__coveredby=boom)[:1]


# splitting c doesn't help
c1 = Polygon.from_bbox((-180, -80, -140, 80))
c2 = Polygon.from_bbox((-140, -80, -100, 80))
test = MultiPolygon(a,c1,c2)
Photo.objects.filter(location__coveredby=test)[:1]
# BOOM! Could not generate outside point!

By changing the numbers I can make this error come and go. (-180, -80, x, 80) works where x <= -140 for example. For every number there is a threshold like this but I can't find a pattern. For boxes with the same area, some will work and others not. For boxes with the same width some will work and others not.


I can look at the SQL being generated but the areas are represented in binary (EWKB) and I'm not sure how to read it.


Can anyone explain this?




Answer



I think I've figured this out, please correct me if my assumptions are wrong.


The MultiPolygon(a,c) covers more than half the earth's surface so postGIS thinks I'm referring to the area outside my MultiPolygon, because it's smaller. This area includes the poles, as Mike Toews pointed out, postGIS doesn't like these extreme values.


My solution is to use an SQL OR to join each Polygon, rather than putting them in a MultiPolygon. My area is split up in the first place to allow for searching areas bigger than half the earth's surface.


Checking if z values of contour lines are correct using ArcGIS Desktop?


What is the easiest way to check if z values of contour lines are correct in ArcGIS?


I need just to look at the values, to see if equidistance between them is correct and unchanged.



Answer



For many contour maps a check of spatial consistency among the contour levels is not possible unless you supply additional information. Here's an example pulled arbitrarily from Google Images:



enter image description here


The criterion for consistency is that the neighboring lines of each contour line should have values differing by at most one contour interval. But if you draw a line segment from the eastern terminus of the 4.0 contour to the eastern terminus of the 4.8 contour, this segment will not intercept any other contour line: thus the 4.8 line is also a "neighbor" of the 4.0 line, even though it is two intervals away. To cope with this problem you need to specify the region being contoured. In this example, it's the polygon bounded by the thick cyan curve.


Assuming you have specified this region as a polygon, then one of the easier ways to do the consistency check is to press Spatial Analyst into service to perform a Euclidean Allocation calculation of the contour lines, using the region polygon as a mask. This expands all lines until they meet their neighbors, but it does not allow the expansion to extend beyond the mask. At any meeting point, only two contour levels will be involved (otherwise you do indeed have a problem!). (You can verify that only two levels ever meet, by computing a focal variety grid and checking that its maximum value is 2.) If you compute a 3 x 3 focal range of the Euclidean allocation grid, it should therefore equal either zero or the contour interval (assuming the interval is constant). That is extremely easy to check by looking at the attribute table or histogram of the focal range grid.


qgis - How to call a method by BOTH a button AND a key shortcut



I have been trying to call a method by BOTH a key shortcut (as described here http://www.qgis.org/en/docs/pyqgis_developer_cookbook/snippets.html) and a button I created on the toolbar.


The button works but when I add the code for the shortcut, nothing happens.


def initGui(self):
# Code for the button
self.my_action = QAction(QIcon(":/path/to/the/icon.png"), u"toolName", self.iface.mainWindow())
self.my_action.triggered.connect(self.myFunction)
self.iface.addToolBarIcon(self.my_action)

# Code for the shortcut
self.iface.registerMainWindowAction(self.my_action, "Ctrl+1")

QObject.connect(self.my_action, SIGNAL("triggered()"),self.myFunction)

def unload(self):
# Unload the button
self.iface.removeToolBarIcon(self.my_action)

# Unload the shortcut
self.iface.unregisterMainWindowAction(self.my_action)

def myFunction(self):

doSomething

Can anybody tell me where (and how) to put the code explained in the link?


Thanks,


m.


EDIT#1


I don't know about my system (if other shortcuts work, I don't see why those set on my tools shouldn't) and as far as my code it is pretty straightforward: in the initGui method I set up the buttons (what I shared is a "place holder" for the real path of the icons. However, I can see them on both the toolbar and the settings form and everything works fine when clicking the buttons) add them to the toolbar; the unload method, unloads them (see the code above WITHOUT the redundant signal/slot connection as dakcarto pointed out in the comment); the run method sets up the environment (it points at the table -- I am editing a PostGIS table -- and starts an editing session + some tuning)


def run(self):
try:
# Reference to the layer

active_layer = self.iface.activeLayer()
# If vector type load a qml file from plugin dir
if active_layer.type() == QgsMapLayer.VectorLayer:
dirPlug = os.path.dirname(os.path.abspath(__file__))
active_layer.loadNamedStyle(dirPlug+"\\file.qml")
# Start editing session
active_layer.startEditing()
# Set CRS and zoom to extent
canvas = self.iface.mapCanvas()
canvas.mapRenderer().setProjectionsEnabled(True)

canvas.mapRenderer().setDestinationCrs(QgsCoordinateReferenceSystem(3857))
canvas.setExtent(active_layer.extent())
# Refresh
canvas.refresh()
self.iface.legendInterface().refreshLayerSymbology(active_layer)
except:
QMessageBox.warning(self.iface.mainWindow(),'WARNING','Please, make sure to select a valid shapefile in the TOC')

Then I have four buttons calling four identical methods that change values in three different fields of selected features (they differ from each other only for the values they are writing into the DB). Everything works fine when using buttons! I just can't connect a shortcut to them! One of these methods:


def oneOfThem(self):

try:
# Reference to the layer
active_layer = self.iface.activeLayer()
# selected features
sel = active_layer.selectedFeatures()
# datetime now
now = datetime.datetime.now().strftime("%Y-%m-%d %H:%M %S")
# username
user = getpass.getuser()
# if something is selected

if sel:
# Change values for selected features
for s in sel:
selID = s.id()
selUser = s.fieldNameIndex('user_lbl')
selRat = s.fieldNameIndex('rating')
selTime = s.fieldNameIndex('tile_done_')

active_layer.changeAttributeValue(selID, selUser, user)
active_layer.changeAttributeValue(selID, selRat, 1)

active_layer.changeAttributeValue(selID, selTime, now)

self.iface.messageBar().pushMessage("Info","Message.", QgsMessageBar.INFO, 3)
else:
self.iface.messageBar().pushMessage("WARNING","No features selected", QgsMessageBar.WARNING, 3)
except:
QMessageBox.warning(self.iface.mainWindow(),'WARNING','Please, make sure to select a valid shapefile in the TOC')

That's pretty much it! I didn't touch the _init_ method and I only added my code to the one coming from the Plugin Builder Plugin.


m.




Answer



You have a redundant signal/slot connection:


def initGui(self):
# Code for the button
self.my_action = QAction(QIcon(":/path/to/the/icon.png"), u"toolName", self.iface.mainWindow())
# Code for the shortcut
self.iface.registerMainWindowAction(self.my_action, "Ctrl+1")
self.my_action.triggered.connect(self.myFunction) # only need this new-style connection
self.iface.addToolBarIcon(self.my_action)


# duplicate old-style connection
# QObject.connect(self.my_action, SIGNAL("triggered()"),self.myFunction)

See if removing the duplicate helps. Also, try a different, unique key sequence in case there is a conflict with the "Ctrl+1" sequence.


Also, check Settings -> Configure shortcuts dialog to see if your registered QAction and key sequence are showing up.


qgis - PostGIS Shapefile Importer Projection SRID


I have a shapefile that I have imported into PostgreSQL using the PostGIS Shapefile Importer and it is currently projected in a State Plane Coordinate System. When I imported it, I selected an import SRID of 4326 (WGS84). Afterwards when I was trying to manipulate, I was getting some odd results because of the fact that my information was still seemingly in a projected coordinate system.


Does the Shapefile Importer not re-project when you assign the SRID while importing?



Answer



you can do a few things





  1. set your srid to your state plane CS then use the st_transform function to convert to WGS84 on your geometry during whatever calculation you are performing -this will only reporject or transform the SRID on the query manipulation you are performing, it will not change the underlining SRID for the original table


    st_transform(ST_SetSRID(geom,'state plane CS'),4326)


  2. change your underlining SRID before you perform the calculations


    SELECT UpdateGeometrySRID('table name','geom',4326);

    ALTER TABLE 'table name'
    ALTER COLUMN geom

    TYPE Geometry(point/line..?, 4326)
    USING ST_Transform(geom, 4326);

    Update 'table name' set geom= st_transform(geom,4326);

    select st_srid(geom) from 'table name';


python - How to create a Shapely LineString from two Points


If have two points, from which I want to create a straight LineString object:


from shapely.geometry import Point, LineString
A = Point(0,0)
B = Point(1,1)


The Shapely manual for LineString states:



A sequence of Point instances is not a valid constructor parameter. A LineString is described by points, but is not composed of Point instances.



So if I have two points A and B, is there a shorter/better/easier way of creating a line AB than my current "best" guess...


AB = LineString(tuple(A.coords) + tuple(B.coords))

... which looks rather complicated. Is there an easier way?


Update


With today's released Shapely 1.3.2, the above statement from the manual is no longer correct. So from now on,



AB = LineString([A, B])

works!



Answer



Since Shapely 1.3, you can create a LineString from Points:


>>> from shapely.geometry import Point, LineString
>>> LineString([Point(0, 0), Point(1, 1)]).wkt
'LINESTRING (0 0, 1 1)'

Apologies for the contradiction in the manual.



arcpy - How to truncate file geodatabase tables with Python?



I need to truncate file geodatabase tables (1 or all) with Python. What is the code for that?



Answer



AFAIK, you can use Delete Rows method in arcpy. from Arcgis Resource Center:


Delete Rows (Data Management)



Summary


Deletes all or the selected subset of rows from the input.


If the input rows are from a feature class or table, all rows will be deleted. If the input rows are from a layer or table view with no selection, all rows will be deleted.



consider this caution:




If run against a layer or table view that does not have a selection, the operation cannot be undone using undo/redo.



Example Code:


import arcpy
from arcpy import env

env.workspace = "C:/data"
arcpy.CopyRows_management("accident.dbf", "C:/output/accident2.dbf")
arcpy.DeleteRows_management("C:/output/accident2.dbf")


i hope it helps you...


pgrouting - Split line at points where other lines touch it in PostGIS


I am trying to create an application that will require a network to be rebuilt several times in an iterative process, adding connections as it progresses. Every iteration needs to use the network created in the last iteration to create the next network.


I am using pgrouting to calculate routes along the network during every iteration, and hoping that building a pgrouting topology does not try to node every intersection, because there are several under/overpasses in the data that I would prefer to simply leave as they are.


However I need to node the new connections that are being made to the existing network, so that I can include them when rebuilding the next network iteration. I have been able to draw a line connecting new points to the network, as well as creating a set of multipoints where the line meets the existing network. I want to split the existing network only at these points, but have all segments retain the attributes of the original line.



Essentially, this is the situation:


enter image description here


I want to manually split the pink line at all the pink dots, so that I can rebuild a pgrouting topology, and I want all the subsequent segments to retain the attributes of the pink line itself.


This is the code that I am using. I just keep receiving the original line's geometry, un-split.


SELECT ST_AsEWKT(ST_Split(N, P)) AS geom
FROM (SELECT
starting_network.geom as N,
cut_points_multi.geom as P
FROM starting_network, cut_points_multi) AS foo;


openlayers 2 - Qgis OSM and layer projection problem


I am using QGIS 2.1.0, and trying to import my vector layer(parcels) with OpenStreetMap layer.
My vector layer (parcels) use EPSG 31277 coordinate system.


This is the definition for EPSG 31277:


(+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs)

OpenStreetMap layer use WGS 84/Pseudo Mercator EPSG 3857.



My project use WGS 84 / Pseudo Mercator EPSG 3857. In project properties I enable 'on the flay' crs transformation.


The problem is that my parcels do not match with OpenStreetMap layer! The difference is about 300 m. Here is a picture that shows that.


enter image description here



Answer



There was some discussion about the right +towgs84 parameters for projections based on MGI Ferro, see Problem with reprojecting raster from MGI 6 to WGS


Projection parameters for Gauss Kruger 7 zone Serbia


As a result, we now have two similar projections in QGIS 2.0.1:


EPSG:3909 +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=682,-203,480,0,0,0,0 +units=m +no_defs

EPSG:31277 +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs


Maybe your data fits better if you assign the EPSG:3909 projection to it with Rightclick -> Set CRS for layer.


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...