QGIS Processing Script for Quick Feature Selection and Zoom

Here's a little QGIS processsing script which might come in handy for you as well, if you've got to visually check many features and need to find/select/zoom them manually by attribute values. Using the attribute table's functions (find by expression, etc.) is rather cumbersome, if you have to iterate over loads of features for visual control - that's where the below script kicks in: You call the processing script from the processing toolbox panel (I like to have this one always visible) and just type in the feature's attribute value you're searching for, then, with one go, this feature is selected and zoomed. The zoom-value (2500) is hard coded into the script as well as the default layer (the one you're currently working with), for the sake of not having to put all of this manually each time (so you'll have to change this for your own purpose;).

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

#===============================================
##[User scripts]=group
##Gst_Nr=string
##Name_Layer=string 81127GST_V2
#===============================================
# Gst_Nr is the field in which we are searching
# 81127GST_V2 is the default layer for searching
#===============================================

canvas = iface.mapCanvas()

# First zoom to desired scale
canvas.zoomScale( 2500 )

allLayers = canvas.layers()
n = len(allLayers)
for i in range(0, n):
    if allLayers[i].name() == Name_Layer:
        break
tarL = allLayers[i]

# Get a featureIterator from an expression
expr = QgsExpression( "\"GNR\"='" + Gst_Nr+ "'" )
it = tarL.getFeatures( QgsFeatureRequest( expr ) )

# Build a list of feature Ids from the result obtained above
ids = [j.id() for j in it]

# Select features with the ids
tarL.setSelectedFeatures( ids )

# Zoom to selected features
canvas.zoomToSelected( tarL )

Offline Map Tiles with Mapbox-Studio and MOBAC

I just found I neat workaround on how to get your Mapbox projects (Mapbox is a free, fully customizeable Tool for map design) to any local tile storage format that it also available in MOBAC (MOBAC is a free software for the creation of local map tiles from variable sources, in a bunch of different formats, to be used on almost any mobile device)!

Get Mapbox-Studio and design your custom map

Upload your project to mapbox

Create an account and connect to it

Get your username, project ID and an access token, which are put to the URL for fetching the tiles back from the Mapbox server: http://api.tiles.mapbox.com/v4/your_username.your_project_id/{$z}/{$x}/{$y}.png?access_token=your_private_or_public_api_key

Use the following XML (replacing user, project ID and access token) for a custom mapsource in MOBAC.

 MAPBOX
 0
 18
 png
 None
 http://api.tiles.mapbox.com/v4/gimoya.ganok9hk/{$z}/{$x}/{$y}.png?access_token=pk.eyJ1IjoiZ2ltb3lhIiwiYSI6IkZrTld6NmcifQ.eY6Ymt2kVLvPQ6A2Dt9zAQ 
 #000000

QspatiaLite Use Case: Find Number of Species from Point Data

Here's a short follow up on some previous posting about the use of QspatiaLite for the aggregation of species distribution data. In this case the species data comes as a point layer. For each cell of a 1000 x 1000 m grid (1) the number of individuals per species, (2) the total number of individuals and (3) the number of different species should be calculated.


There is a layer with 10 different species (variabel name is "sp") across the whole extent with names "1", "2", "3", .. , "10" and a layer with the grid cells (variable name = "id") numbered consecutively, from 1 to 150.

In the attribute table of the below screenshot you see that I selected the grid cell with id=1 and the points (=species) within this cell. There are 8 individuals - "4", "5", "6" and "10" occure once, whereas "2" and "8" occure twice.

The query table in the screenshot is the result for (3).

For (1) you have to query for points/species within grid cells and group over grid cells and species and take the count from that aggregation

SELECT
 t.gID AS gID,
 t.sp AS Sp,
 count(*) AS NrIndSp
FROM (SELECT
  g.id AS gID, 
  s.sp AS sp
FROM grid AS g JOIN Sp_distr AS s 
ON within(s.Geometry, g.Geometry)
) as t GROUP BY t.gId, t.sp

For (2) you simple need to query for points/species within grid cells and aggregate over grid cells:

Select 
 t.gID,
 count(*) as NrInd
From (SELECT
  g.id AS gID, 
  s.sp AS sp
FROM grid AS g JOIN Sp_distr AS s 
ON within(s.Geometry, g.Geometry)
) as t 
GROUP BY t.gID 
ORDER BY t.gID

For (3) you'll first need to aggregate over grid cells and points/species, and then again aggregate over this query table by grid cells which will finally give you the distinct species!

SELECT 
 v.gID,
 count(*) AS SpNr
FROM (SELECT 
 t.gID,
 t.sp
FROM (SELECT
  g.id AS gID, 
  s.sp AS sp
FROM grid AS g JOIN Sp_distr AS s 
ON within(s.Geometry, g.Geometry)
) as t GROUP BY t.gId, t.sp
) as v GROUP BY v.gId

Python Processing Script for Splitting Layer and Saving each Feature to Shapefile

A script for a common task: Split a Layer by distinct values of an attribute-field. The script will use the value for the name of the files and remove the field from the new files.

Put this processing script to the appropiate director (in my case it's "C:\Users\Kay\.qgis2\processing\scripts") to make it work. Then it will be easily accessible from your processing toolbox.

##[User scripts]=group
##input=vector
##class_field=field input
##output=output file

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

layer = processing.getObject(input)
provider = layer.dataProvider()
fields = provider.fields()
class_field_index = layer.fieldNameIndex(class_field)

fields.remove( class_field_index )

inFeat = QgsFeature()
outFeat = QgsFeature()
inGeom = QgsGeometry()
nElement = 0
writers = {}

feats = processing.features(layer)
nFeat = len(feats)

for inFeat in feats:
    progress.setPercentage(int(100 * nElement / nFeat))
    nElement += 1
    featAttr = inFeat.attributes()
    classField = featAttr[class_field_index]

    if classField not in writers:
        outputFile = output + '_' + classField + '.shp'
        writers[classField] = VectorWriter(outputFile, None, fields, provider.geometryType(), layer.crs())
        
    inGeom = inFeat.geometry()
    outFeat.setGeometry(inGeom)
    
    del featAttr[class_field_index]
    outFeat.setAttributes(featAttr)
        
    writers[classField].addFeature(outFeat)

for writer in writers.values():
    del writer

QspatiaLite Use Case: Connecting Lines

With QSpatiaLIte you can connect disjoint lines quite easily. With the below SQL you can allow for a grouping variable, in this case the field 'name' within the layer 'segments', by which the group vertices are collected and connected as lines! With this approach the vertices are connected in the order in which they were digitized and existing gaps are closed.

select 
  name as name,
  makeLine(t.ptgeom, 0) as geom 
from (
     select
        name as name,
        DissolvePoints(Collect(ST_Reverse(geometry)))  as ptgeom
     from segments group by name )
as t

Usecase for KML-Parsing: Make New KML-File from File-Collection

In this usecase I had collected several KMLs from the internet but wanted to strip them down for only the relevant parts (the Linestrings inside the Placemark-nodes) and put them all inside one final File. In my script I create a new KML file and populate a folder-node inside it with Linestrings from the collection of KML-files which all reside in the same source directory. For this one needs to parse each file and grab the appropiate nodes and add them to the target kml file. In addition I alter some oroginal values, i.e. I use the file names of the single KML-files as Placemark names inside the new KML-file.

Here is the final file as seen after opening in Google Earth:

library(XML)

# new kml file... needs to be well-formed
z <-
  '<kml xmlns="http://www.opengis.net/kml/2.2">
      <Document>
         <Folder>
            <name>ROUTES</name>
         </Folder>
      </Document>
    </kml>'
new_xmlDoc <- xmlInternalTreeParse(z, useInternalNodes = TRUE)

# important add all namespace definitions...
ns <- c(gx="http://www.google.com/kml/ext/2.2",
        kml="http://www.opengis.net/kml/2.2",
        atom="http://www.w3.org/2005/Atom")
ensureNamespace(new_xmlDoc, ns)

# get the root off the new file for latter processing
new_root <- xmlRoot(new_xmlDoc)

# loop over files from folder
# and insert Placemark content of each file as children nodes into 
# the new file

setwd("C:/Users/Kay/Google Drive/SKI-BIKE/Gastein")
files <- dir(pattern="bergfex*")

for (f in files) { 
   
   # get placemark node of each file
   doc <- xmlInternalTreeParse(f, useInternalNodes = TRUE)
   root <- xmlRoot(doc)
   plcm_node <- root[["Document"]][["Folder"]][["Folder"]][["Placemark"]]

   # insert file name as Placemark name
   xmlValue(plcm_node[["name"]]) <- sub('bergfextour_(.*)[.]kml', '\\1', f)

   # add placemark node to new doc
   addChildren(new_root[["Document"]][["Folder"]], plcm_node)

}

# save it...
saveXML(new_xmlDoc, "collapsed_ROUTES.kml")

How to Integrate OpenStreetMap Nominatim Search into a Webpage or OpenLayers Map

Here' a concise example on how to integrate Nominatim search into a webpage:
http://jsfiddle.net/gimoya/yp8hmog7

And here's a live Openlayers example:
http://gimoya.bplaced.net/Openlayers/Bergfex/bergfex.html