21 Oct 2015

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 )

25 Apr 2015

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
    
    
  • 8 Apr 2015

    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
    

    18 Feb 2015

    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
    

    6 Feb 2015

    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
    

    5 Feb 2015

    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")