Workflow for Replacing DOM With DEM Values in a Zone Specified by Vector Layer

1. Create polygon file with same CRS as DOM and DEM and draw polygon mask (obv. polygons with holes don't work!) 

2. Use the "Rasterize (overwrite with fixed value)" algorithm, using the polygon mask as vector layer and the DOM as raster layer (set the DOM as reference layer!). Set the burn value to -9999. BEAWARE: The original file will be overwritten - if you don't want to alter the original file, make a copy of the DOM first!

3. Use the raster calculator to rewrite the values in the DEM with values from the DOM by using this formula: 

(DOM@1 = -9999) * DGM@1 + (DOM@1 != -9999) * DOM@1

Windows Batch Script for Zipping Shapefile Components

Windows Batch File to zip ESRI shapefile-components and delete original files. Put the code to a text file and save with .bat extension. Save the file to the folder you want to run the commands.

@ECHO OFF
set /p $dum="Hit enter to zip Shapefiles in %~dp0 ..."
FOR %%F IN (*.shp) DO "C:\Program Files\7-Zip\7zG.exe" a "%%~nF.zip" "%%~nF.shp" "%%~nF.dbf" "%%~nF.prj" "%%~nF.shx"
for /f "delims=" %%F in ('dir /b /a-d ^| findstr /vile ".bat .zip"') do Echo "%%F"
set /p $dum="Hit enter to delete original files, listed above.."
for /f "delims=" %%F in ('dir /b /a-d ^| findstr /vile ".bat .zip"') do del "%%F"
Echo Done !!
set /p $dum="Hit [Enter] to exit..."

QGIS 3: Symbology with Geometry Generator - Draw One Convex Hull For All Features With Same Attribute

I have a layer named "Trails" with an attribute "Trail ID", which contains unique, consecutive Feature/Trail IDs and an attribute "Schwierigkeit" (trail difficulty). For all features with the value "black" for the attribute "Schwierigkeit" I want to render one convex hull. The below code will select the last element of all the features in the layer and apply the code for drawing the polygon only once. The first

array_foreach()

in the code will create an array of all features (series generated from 1 to feature count number). Over this array, the second array_foreach() will apply the geometry function on each element that meets the condition of the

array_filter()

function. The

collect_geometries()

function finally puts all those single geometries within the resulting array into one multiline geometry, for which I then draw the hull. The purpose of this procedure, is to check if the trails in my dataset show a spatial aggregation according to their trail difficulty..

if($id = maximum($id),
convex_hull(
collect_geometries(
with_variable('my_arr', 
array_foreach(
generate_series(1,  layer_property( 'trails', 'feature_count'), 1),
get_feature('trails', 'Trail ID', @element)
),
array_foreach(
array_filter(@my_arr, attribute(@element, 'Schwierigkeit')='black') , geometry(@element)))))
, NULL)

QGIS 3: Layout expression to get all features of a layer within the current map view

In a map layout template you can insert an expression which gets all feature's names of a certain layer, which are contained in the current layout's view/extent. For wxample you could concatenate all feature's names contained in the current map view like so:

aggregate(layer:='Corridors_86ef2ea9_d1a1_4d9f_9735_7b7b2fb54cb2', 
aggregate:='concatenate', 
expression:="Name", filter:=within($geometry, map_get( item_variables('Main Map'), 
'map_extent')), concatenator:=', ')

Aggregation of Different Layers in QGIS 3.x with Field Calculator Expression Alone!!

Here's an example of how to "spatially aggregate" Polygon attribute values over features of a point layer by intersecting the two layers. Go in "Edit mode" with the target layer and check "Virtual Field" with type "decimal". Then use this in expression builder, where "layer" is the point layer, and aggregate the attribute values that intersect the target layer's geometry. With "Expression" you can define the attribute, that should aggregated!

aggregate(layer:='BIKFFH_Karwendel BIKFFH_PL', aggregate:='sum', expression:="LNUMMER", 
filter:=intersects($geometry,geometry(@parent))) 

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

QspatiaLite Quicktip: Convert MULTILINESTRING to LINESTRING

One often encounters the problem, that after digitizing or running processing algorithms, the output geometry is MULTILINESTRING, but we rather wished to have the geometrytype LINESTRING. Until know I used a quite cumbersome, multistep workflow for conversion between these geometry-types - however, as we will see, all of this becomes ridicously easy with spatial SQL:

select 
   replace(replace(replace(replace(replace(replace(astext(Collect(t.geometry)), 'MULTILINESTRING((','§'), '))', '%'), '(', ''), ')', ''), '§', 'LINESTRING('), '%', ')'
) as geom
from (
    select MultiLinestringFromText('MULTILINESTRING((-1 -1, 0 0), (1 1, 4 4))') as geometry
) as t

resulting in:

LINESTRING(-1 -1, 0 0, 1 1, 4 4)

However, if your orginal line was something like MULTILINESTRING((-1 -1, 0 0), (0 0, 4 4))
you'd end up with:

LINESTRING(-1 -1, 0 0, 0 0, 4 4)

which contains double vertices, which we certainly don't want!

So be aware, that the ordering / direction of the linestring will be as in the segments of the original layer! And as we saw, gaps between subsequent end-/startnodes will be closed in the new geometry!! It is adviseable to doublecheck before / after conversion!

If you deal with a multilinestring (or a combination of any type of linesstrings) which share end/startnodes nodes things are even easieruse this SQL:

SELECT AsText(Linemerge(MultiLinestringFromText('MULTILINESTRING((-1 -1, 0 0), (0 0, 4 4))')))

resulting in:

LINESTRING(-1 -1, 0 0, 4 4)

QspatiaLite Use Case: Query for Species Richness within Search-Radius

Following up my previous blogpost on using SpatiaLite for the calculation of diversity metrics from spatial data, I'll add this SQL-query which counts unique species-names from the intersection of species polygons and a circle-buffer around centroids of an input grid. The species number within the bufferarea are joined to a newly created grid. I use a subquery which grabs only those cells from the rectangular input grid, for which the condition that the buffer-area around the grid-cell's centroid covers the species unioned polygons at least to 80%.

Example data is HERE. You can use the shipped qml-stylefile for the newly generated grid. It labels three grid-cells with the species counts for illustration.

Import grid- and Sp_distr-layers with QspatiaLite Plugin

Run query and choose option "Create spatial table and load in QGIS", mind to set "geom" as geometry column

select 
    g1.PKUID as gID,
    count (distinct s.species) as sp_num_inbu, 
    g1.Geometry AS geom
from (
 select g.*
 from(select Gunion(geometry) as geom
           from Sp_distr) as u, grid as g
 where area(intersection(buffer(centroid(g.geometry), 500), u.geom)) > pow(500, 2)*pi()*0.8
) as g1 join Sp_distr as s on intersects(buffer(centroid( g1.Geometry), 500), s.Geometry)
group by gID

QspatiaLite Use Case: Get Subselection of Grid which Covers Polygon

Here's another short SQL-query which I used to get a subselect from a rectengular grid. Aim is to keep only the grid-cells that fully cover the area of a second polygon-layer - cells which do not overlap the polygon's area completely will be skipped from the new grid-layer.

select 
  g.*
from(select Gunion(geometry) as geom
           from MYPLGN) as u, grid as g
where area(intersection(g.geometry, u.geom)) = area(g.geometry)

QspatiaLite Use Case: Connect Points with Same ID with Line Using the QspatiaLite Plugin

Another short example illustrating the effectiveness of geoprocessing with SpatiaLite, using the great QGIS-plugin QspatialLite.

We have a point-layer with an ID column ("Birds"), with each ID occuring twice, each ID representing an individual. The Ids should be used as start- & end-nodes for the connecting lines. Note that this also would apply if there were more than two points - then the same query could be used to connect all bird individual's points to a line by the order in each group!

We want each set of points, grouped by ID, to be connected. This is easily achieved by importing the points to a SpatiaLite-DB with the QspatiaLite plugin and running a very simple query:

SELECT 
    ID,
    makeline(Geometry) AS geom
FROM Birds
GROUP BY ID

Load the result to QGIS and that's it!

QspatiaLite Use Case: Find Dominant Species and Species Count within Sampling Areas Using the QspatiaLite Plugin

This blogpost shows how to find the dominant species and species counts within sampling polygons. The Species-layer that I'll use here is comprised of overlapping polygons which represent the distribution of several species. The Regions-layer represents areas of interest over which we would like to calculate some measures like species count, dominant species and area occupied by the dominant species.

Since QGIS now makes import/export and querying of spatial data easy, we can use the spatiaLite engine to join the intersection of both layers to the region table and then aggregate this intersections by applying max- and count-function on each region. We'll also keep the identity and the area-value of the species with the largest intersecting area.

For the presented example I'll use

Regions, which is a polygon layer with a areas of interest

Species, which is a polygon layer with overlapping features, representing species

Do the calculation in 2 easy steps:

Import the layers to a spatiaLite DB with the Import function of the plugin (example data: HERE)

Run the query. For later use you can load this table to QGIS or export with the plugin's export button.

SELECT   
  t.region AS region,
  t.species AS sp_dom,
  count(*) AS sp_number,
  max(t.sp_area) / 10000 AS sp_dom_area
  FROM ( 
      SELECT
          g.region AS region, s.species AS species,
          area(intersection(g.Geometry, s.Geometry)) AS sp_area
          FROM Regions AS g JOIN Sp_Distribution AS s 
          ON INTERSECTS(g.Geometry, s.Geometry)  
  ) AS t
GROUP BY t.region
ORDER BY t.region

Addendum:
If you wish to calculate any other diversity measures, like Diversity- or Heterogenity-Indices, you might just run the below query (which actually is the subquery from above) and feed the resulting table to any statistic-software!

The output table will contain region's IDs, each intersecting species and the intersection area.
The intersection area, which is the species' area per polygon, is the metric that would be used for the calculation of diversity / heterogenity measures, etc. of regions!

SELECT 
  g.region AS regID, 
  s.species AS sp,
  AREA(INTERSECTION(g.geometry, s.geometry)) AS sp_area
FROM Regions AS g JOIN Sp_Distribution AS s 
ON INTERSECTS(g.Geometry,s.Geometry)
ORDER BY regID, sp_area ASC

I tested this on

QGIS 2.6 Brighton

with the QspatiaLite Plugin installed

QspatiaLite Use Case: SpatiaLite Aggregation over Points within Polygons using the QspatiaLite Plugin

Here's a nice example for aggregation of points per polygon areas, which I grabbed from an Answer on SO, by user @Micha. The polygons could be regions of interest, a sampling grid, etc.
Say you want to do maximum, minimum, averages, etc. per polygon using the spatial database SpatiaLite.

You'll first need to import both of your layers into a spatialite DB, called "sensors" (the point layer) here, with a "pollution" column and "SHAPE1" (the polygons) with a "plgnID" column. You can do this easily with the QspatiaLite-plugin "Import" button (example data is HERE).

Now this query will give you various statistics from the sensors for each polygon:

SELECT g.plgnID AS "plgn_ID",
   AVG(s.pollution) AS "Average Pollution", 
   MAX(s.pollution) AS "Maximum Pollution",
   COUNT(*) AS "Number of Sensors"
FROM sensors AS s JOIN SHAPE1 AS g 
ON contains(g.geometry, s.geometry)
GROUP BY g.plgnID

QspatiaLite Use Case: SpatiaLite Aggregation over Intersections of Polygons with QspatiaLite Plugin

This applies to several usecases: Imagine you have a grid or polygon-layer of sampling areas and want to know the dominant feature of another polygon layer under each grid cell / sampling polygon - this could be soiltypes, landuse classes, etc. Other than the dominant feature you might be interested in the diversity of features (i.e. number of soils, etc.) per grid cell / sampling area.

QGIS alone does not provide handy tools for aggregation of features of one layer combined with other layers, but the spatiaLite engine is tailored for this! Since QGIS now makes import/export and querying of spatial data easy, it seems very worthy to dive into spatiaLite and utilize its powerful tools!


For the presented example I'll use:

SHAPE1, which is a polygon layer with a sampling grid/areas

Soils, which is a polygon layer with soiltypes

I tested this on

QGIS 2.6 Brighton

with the QspatiaLite Plugin installed

Import the above layers to a spatiaLite DB with the Import function of the plugin (example data: HERE)

Run the query and choose "create spatial table and load in QGIS" and put geom as geometry column! (I chose SHAPE2 as name for the newly created layer..)

SELECT t.geom AS geom, 
    t.plgnID AS plgnID, 
    t.soiltype AS soiltype, 
    max(t.soil_area) AS MaxArea, count () AS n_soiltypes
       FROM (SELECT 
          g.Geometry AS geom, g.plgnID AS plgnID, s.Soiltype AS soiltype,
          AREA(INTERSECTION(g.geometryO, s.geometry)) AS soil_area
          FROM SHAPE1 AS g JOIN Soils AS s 
          ON INTERSECTS(g.Geometry,s.Geometry)
       ) AS t
GROUP BY t.plgnID
ORDER BY t.plgnID

That's it!

Make a KML-File from an OpenStreetMap Trail

Ever wished to use a trail on OSM on your GPS or smartphone? With this neat little R-Script this can easily be done. You'll just need to search OpenStreetMap for the ID of the trail (way), put this as argument to osmar::get_osm, convert to KML and you're good to go!

# get OSM data
library(osmar)
library(maptools)
  
rotewandsteig <- get_osm(way(166274005), full = T)
sp_rotewandsteig <- as_sp(rotewandsteig, what = "lines")
  
# convert to KML 
kmlLine(sp_rotewandsteig@lines[[1]], kmlfile = "rotewandsteig.kml",
        lwd = 3, col = "blue", name = "Rotewandsteig") 

# view it
shell.exec("rotewandsteig.kml")

Custom Feature Edit Forms in QGIS with Auto-Filling Using PyQt

For anyone interested in the capabilities of customized feature edit forms in QGIS I'd like to reference the following GIS-Stackexchange posting: http://gis.stackexchange.com/questions/107090/auto-populate-field-as-concatenation-of-other-fields

Quick GIS-Tip: Permanentely Sort Attribute Table by Attribute

Here's a short shell-script (I have C:\OSGeo4W64\bin\ on my PATH) for sorting GIS data, a sqlite-db in this case, and saving it in the newly created order to a file. The DB is sorted in ascending order by the attribute 'Name' and written to a file with the KML Driver.

cd C:/Users/Kay/Documents/Web/Openlayers/Docs/Digitizing
ogr2ogr -sql "SELECT * FROM 'trails-db' ORDER BY Name ASC" -f "KML" trails.kml trails.sqlite

R GIS: Function to Reverse KML Paths

This is a function I wrote up for reversing KML-paths. The paths within a KML can be partially matched by their name-tags

## name:          ReverseKmlPath   
## use:           Reverse KML-pathsby matching their Name tags
## arguments:     PATH_TO_DOC, the path to the KML-file
##                NAME, the value of the name tag, function uses partial matching!
##                'Trail_xyz' will be matched by 'rail'
## requirements:  KML-structure with Placemarks containing a  and a  tag
## author:        Kay Cichini
## date:          01-05-2014
## license:       CC-BY-NC-SA

ReverseKmlPath <- function(PATH_TO_DOC, NAMES) {
    
    require(XML)

    doc <- xmlInternalTreeParse(PATH_TO_DOC)
    
    if (xmlNamespaceDefinitions(doc)[[1]]$uri == "http://www.opengis.net/kml/2.2") {
        namespaces <- c(kml = "http://www.opengis.net/kml/2.2")
        flag <- 1
    } else {
        if (xmlNamespaceDefinitions(doc)[[1]]$uri == "http://earth.google.com/kml/2.0") { 
                namespaces <- c(kml0 = "http://earth.google.com/kml/2.0")
                flag <- 0
            } else {
                stop ("Stopped!: Check namespace issue..")
            }
    }
        
    
    for (NAME in NAMES) {
        
        if (flag) { 
              query <- paste0("//kml:Placemark[contains(kml:name,'", sprintf("%s", NAME), "'", ")]//kml:coordinates")
          } else {
              query <- paste0("//kml0:Placemark[contains(kml0:name,'", sprintf("%s", NAME), "'", ")]//kml0:coordinates")
          }

        coords <- tryCatch(getNodeSet(doc, query, namespaces), 
                           error = function(e) message(paste("\nError: *", NAME, "* was NOT successfully matched\n")))
        
        for (i in length(coords)) {

            #grab coordinates from node and reverse order
            rev_coord_vector <- rev(unlist(strsplit(gsub("\\t|\\n", "", xmlValue(coords[[i]])), "\\s")))
            rev_coord_string <- paste(rev_coord_vector, collapse = " ")

            # re-insert reversed line-string:
            xmlValue(coords[[i]]) <- rev_coord_string

            # message
            if (flag) { 
                  query <- paste0("//kml:Placemark[contains(kml:name,'", sprintf("%s", NAME), "'", ")]//kml:name")
              } else {
                  query <- paste0("//kml0:Placemark[contains(kml0:name,'", sprintf("%s", NAME), "'", ")]//kml0:name")
            }
            match <- xmlValue(getNodeSet(doc, query, namespaces)[[i]])
            message(paste0("matched name: ", match, "\n..."))

        }
    }

    # save:
    message("Reversed paths saved to:")
    saveXML(doc, paste0(dirname(PATH_TO_DOC), "/reversed_", basename(PATH_TO_DOC)),
            prefix = newXMLCommentNode("This file was created with the R-package XML::saveXML, see: "))
}

## not run: 
tf <- tempfile(fileext = ".kml")
download.file("http://dev.openlayers.org/releases/OpenLayers-2.13.1/examples/kml/lines.kml", tf, mode = "wb")
ReverseKmlPath( PATH_TO_DOC = tf, NAMES = c("Absolute", "Relative") )

shell.exec(tf)
shell.exec(paste0(dirname(tf), "/reversed_", basename(tf)))

R GIS: Generalizer for KML Paths

I'm posting a recent project's spin-off, which is a custom line-generalizer which I used for huge KML-paths. Anyone with a less clumpsy approach?

## line generalizing function: takes two vectors of with x/ycoords 
## and return ids of x/y elements which distance to its next element
## is shorter than the average distance between consecutive vertices
## multiplied by 'fac'
check_dist <- function(x, y, fac) {
    dm <- as.matrix(dist(cbind(x, y)))
    
    ## supradiagonal holds distance from 1st to 2nd, 2nd to 3rd, etc. element
    d <- diag(dm[-1, -ncol(dm)])
    mean_dist <- mean(d)
    keep <- logical()
    
    ## allways keep first..
    keep[1] <- T
    for (i in 1:(length(x) - 2)) {
        keep[i + 1] <- (d[i] > mean_dist * fac)
        message(paste0("Distance from item ", i, " to item ", i + 1, " is: ", d[i]))
    }
    message(paste0("Treshold is: ", mean_dist * fac))
    cat("--\n")
    ## .. and always keep last
    keep[length(x)] <- T
    return(keep)
}

## Testing function check_dist:
x <- rnorm(5)
y <- rnorm(5)
(keep <- check_dist(x, y, 1.2))

plot(x, y)
lines(x[keep], y[keep], lwd = 4, col = "green")
lines(x, y, lwd = 1, col = "red")
text(x, y + 0.1, labels = c(1:length(x)))


## exclude vertices by generalization rule. coordinate-nodes with low number of vertices, 
## segments with less than 'min_for_gen' vertices will not be simplified, in any case coordinates will be
## rounded to 5-th decimal place

generalize_kml_contour_node <- function(node, min_for_gen, fac) {
    
    require(XML)
    
    LineString <- xmlValue(node, trim = T)
    
    LineStrSplit <- strsplit(unlist(strsplit(LineString, "\\s")), ",")
    
    # filter out empty LineStrings which result from strsplit on '\\s'
    LineStrSplit <- LineStrSplit[sapply(LineStrSplit, length) > 0]
    
    # all 3 values are required, in case of error see for missing z-values:
    x <- round(as.numeric(sapply(LineStrSplit, "[[", 1, simplify = T)), 5)
    y <- round(as.numeric(sapply(LineStrSplit, "[[", 2, simplify = T)), 5)
    z <- round(as.numeric(sapply(LineStrSplit, "[[", 3, simplify = T)), 5)
    
    # for lines longer than 'min_for_gen' vertices, generalize LineStrings
    if (length(x) >= min_for_gen) {
        keep <- check_dist(x, y, fac)
        x <- x[keep]
        y <- y[keep]
        z <- z[keep]
        xmlValue(node) <- paste(paste(x, y, z, sep = ","), collapse = " ")
        
        # for all other cases, insert rounded values
    } else {
        xmlValue(node) <- paste(paste(x, y, z, sep = ","), collapse = " ")
    }
}

## mind to use the appropiate namespace definition: alternatively use: 
## c(kml ='http://opengis.net/kml/2.2')
kml_generalize <- function(kml_file, min_for_gen, fac) {
    doc <- xmlInternalTreeParse(kml_file)
    nodes <- getNodeSet(doc, "//kml:LineString//kml:coordinates", c(kml = "http://earth.google.com/kml/2.0"))
    mapply(generalize_kml_contour_node, nodes, min_for_gen, fac)
    saveXML(doc, paste0(dirname(kml_file), "/simpl_", basename(kml_file)))
}

## get KML-files and generalize them
kml_file <- tempfile(fileext = ".kml")
download.file("http://dev.openlayers.org/releases/OpenLayers-2.13.1/examples/kml/lines.kml", 
              kml_file, mode = "wb")
kml_generalize(kml_file, 5, 0.9)
shell.exec(kml_file)
shell.exec(paste0(dirname(kml_file), "/simpl_", basename(kml_file)))

Use Case: Make Contour Lines for Google Earth with Spatial R

Here's comes a script I wrote for creating contour lines in KML-format to be used with Google Earth http://github.com/gimoya/theBioBucket-Archives/blob/master/R/contours_for_google_earth.R

If you want to check or just use the datasets I created for the Alps region, you can download it here: http://terrain-overlays.blogspot.co.at/index.html

Use GDAL from R Console to Split Raster into Tiles

When working with raster datasets I often encounter performance issues caused by the large filesizes. I thus wrote up a little R function that invokes gdal_translate which would split the raster into parts which makes subsequent processing more CPU friendly. I didn't use built-in R functions simply because performance is much better when using gdal from the command line..

The screenshot to the left shows a raster in QGIS that was split into four parts with the below script.

## get filesnames (assuming the datasets were downloaded already. 
## please see http://thebiobucket.blogspot.co.at/2013/06/use-r-to-bulk-download-digital.html 
## on how to download high-resolution DEMs)
setwd("D:/GIS_DataBase/DEM")
files <- dir(pattern = ".hgt")

## function for single file processing mind to replace the PATH to gdalinfo.exe!
## s = division applied to each side of raster, i.e. s = 2 gives 4 tiles, 3 gives 9, etc.
split_raster <- function(file, s = 2) {
    
    filename <- gsub(".hgt", "", file)
    gdalinfo_str <- paste0("\"C:/OSGeo4W64/bin/gdalinfo.exe\" ", file)
      
    # pick size of each side
    x <- as.numeric(gsub("[^0-9]", "", unlist(strsplit(system(gdalinfo_str, intern = T)[3], ", "))))[1]
    y <- as.numeric(gsub("[^0-9]", "", unlist(strsplit(system(gdalinfo_str, intern = T)[3], ", "))))[2]
    
    # t is nr. of iterations per side
    t <- s - 1
    for (i in 0:t) {
        for (j in 0:t) {
            # [-srcwin xoff yoff xsize ysize] src_dataset dst_dataset
            srcwin_str <- paste("-srcwin ", i * x/s, j * y/s, x/s, y/s)
            gdal_str <- paste0("\"C:/OSGeo4W64/bin/gdal_translate.exe\" ", srcwin_str, " ", "\"", file, "\" ", "\"", filename, "_", i, "_", j, ".tif\"")
            system(gdal_str)
        }
    }
}

## process all files and save to same directory
mapply(split_raster, files, 2) 

Use Case: Spatial R & Google Earth for Terrain Analyses

I'd like to share code that uses spatial R and Google Earth for terrain analyses. In this example I took SRTM data at 1" resolution from http://www.viewfinderpanoramas.org/dem3.html#alps read it into R did a little processing and finally wrapped it up in a KML-file that I use as ground-overlay in Google Earth. In fact I eventually converted it into a sqlitedb with MAPC2MAPC for usage on a mobile device.

See the code here on github..

R GIS: Polygon Intersection with gIntersection{rgeos}

A short tutorial on doing intersections in R GIS. gIntersection{rgeos} will pick the polygons of the first submitted polygon contained within the second poylgon - this is done without cutting the polygon's edges which cross the clip source polygon. For the function that I use to download the example data, url_shp_to_spdf() please see HERE.

library(rgeos)
library(dismo)

URLs <- c("http://gis.tirol.gv.at/ogd/umwelt/wasser/wis_gew_pl.zip",               # all water bodies in Tyrol
          "http://gis.tirol.gv.at/ogd/umwelt/wasser/wis_tseepeicher_pl.zip")       # only artificial..

y <- lapply(URLs, url_shp_to_spdf)
z <- unlist(unlist(y))
a <- getData('GADM', country = "AT", level = 2)

b <- a[a$NAME_2=="Innsbruck Land", ]                                               # political district's boundaries
c <- spTransform(b, z[[1]]@proj4string)                                            # (a ring polygon)    
z1_c <- gIntersection(z[[1]], c, byid = TRUE)                                      
z2_c <- gIntersection(z[[2]], c, byid = TRUE)

plot(c)
plot(z1_c, lwd = 5, border = "red", add = T)
plot(z2_c, lwd = 5, border = "green", add = T)
plot(z[[1]], border = "blue", add = T)              # I plot this on top, so it will be easier to identify
plot(z[[2]], border = "brown", add = T)

Batch Downloading Zipped Shapefiles with R

Here's a function I use to download multiple zipped shapefiles from url and load them to the workspace:

URLs <- c("http://gis.tirol.gv.at/ogd/umwelt/wasser/wis_gew_pl.zip",
          "http://gis.tirol.gv.at/ogd/umwelt/wasser/wis_tseepeicher_pl.zip")

url_shp_to_spdf <- function(URL) {

  require(rgdal)

  wd <- getwd()
  td <- tempdir()
  setwd(td)

  temp <- tempfile(fileext = ".zip")
  download.file(URL, temp)
  unzip(temp)

  shp <- dir(tempdir(), "*.shp$")
  lyr <- sub(".shp$", "", shp)
  y <- lapply(X = lyr, FUN = function(x) readOGR(dsn=shp, layer=lyr))
  names(y) <- lyr

  unlink(dir(td))
  setwd(wd)
  return(y)
  }

y <- lapply(URLs, url_shp_to_spdf)
z <- unlist(unlist(y))

# finally use it:
plot(z[[1]])

Follow Up on Spatial Overlays with R - Getting Altitude for a Set of Points

A short follow up on a previous post on spatial overlays with R.

library(sp)
library(dismo)

# some addresses in Austria
pts <- geocode(c("Aldrans, Grubenweg", "Wien, Stephansdom", "Salzburg, Mozartplatz"))
 
# make pts spatial
coords <- SpatialPoints(pts[, c("longitude", "latitude")])
spdf_pts <- SpatialPointsDataFrame(coords, pts)

# assign CRS/projection (which is WGS 1984)
crs <- CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") 
proj4string(spdf_pts) <- crs
 
# spatial data to extract (altitude)
alt <- getData('alt', country = "AT")

# convert alt from raster to grid (needed for over::sp)
# and assign CRS (which is the same as spdf_pts, see > alt@crs)
# don't mind warning - the CRS is the same..
spdf_alt <- as(alt, 'SpatialGridDataFrame')
proj4string(spdf_alt) <- crs

# view
plot(alt)
# plot pts on top
plot(spdf_pts, cex = 2, col = 2, add = T)
 
# check data
str(spdf_pts)
str(spdf_alt)
 
# get the raster/pixel/grid data (> ?over):
cbind(spdf_pts$interpretedPlace, over(spdf_pts, spdf_alt))

# result:
#                                       spdf_pts$interpretedPlace AUT_msk_alt
# 1                              Grubenweg, 6071 Aldrans, Austria         736
# 3 Saint Stephen's Vienna, Stephansplatz 1, 1010 Vienna, Austria         183
# 2                           Mozartplatz, 5020 Salzburg, Austria         450

Loading Multiple Shapefiles to the R-Console Simultaneously

A quick tip on how to load multiple shapefiles (point shapefiles, i.e.) to the R console in one go:

library(maptools)

# get all files with the .shp extension from working directory 
setwd("D:/GIS_DataBase/GIS_Tirol/Tirol_Verbreitungskarten/Verbreitungs_Daten")

shps <- dir(getwd(), "*.shp")

# the assign function will take the string representing shp and turn it into a variable
# which holds the spatial points data
for (shp in shps) assign(shp, readShapePoints(shp))
plot(get(shp[1])) # i.e.
# ...done

Quick Tip for Austrian QGIS Users: How to Use Bergfex (OEK 50) Topographic Map as Base Map in QGIS

In the below stackexchange postings you will learn how to get a custom tile source running on QGIS using the OpenLayers Plugin.

http://gis.stackexchange.com/questions/69907/easiest-way-to-create-historical-street-maps-using-qgis
http://gis.stackexchange.com/questions/56615/how-to-add-another-layer-origin-in-my-installed-openlayers-plugin
http://gis.stackexchange.com/questions/70648/using-custom-source-for-qgis-openlayers-plugin

For the Bergfex tile source you can find the working html-file here:
https://gist.github.com/gimoya/6450983

Quick Tip for Austrian MOBAC Users: A Link to a Mapsource for Bergfex Topographic Map!

Here's a link to a Beanshell script enabling Bergfex topographic maps (OEK 50) as mapsource in MOBAC. With this mapsoursce and MOBAC you can create really nice offline topographic maps for Austria and use the created tiles for your special purposes.

Spatial Overlays with R - Retrieving Polygon Attributes for a Set of Points

A short tutorial for spatial overlays using R-GIS..

library(sp)
library(dismo)
 
# spatial data (political districts of Austria)
gadm <- getData('GADM', country = "AT", level = 2)

# view
plot(gadm)
 
# some addresses
pts <- geocode(c("Aldrans, Grubenweg", "Wien, Stephansdom", "Salzburg, Mozartplatz"))
 
# make pts spatial
coords <- SpatialPoints(pts[, c("longitude", "latitude")])
spdf_pts <- SpatialPointsDataFrame(coords, pts)

# assign CRS/projection (which is WGS 1984)
crs <- CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") 
proj4string(spdf_pts) <- crs
 
# check data
str(spdf_pts)
str(gadm)
 
# plot pts on top
plot(spdf_pts, cex = 2, col = 2, add = T)
 
# do an intersection (points in polygon)
# yielding the polygon's attribute data
over(spdf_pts, gadm)

R GIS: Terrain Analysis for Polygons as Simple as it Gets!

library(rgdal)
library(raster)

alt <- getData('alt', country = "AT")
gadm <- getData('GADM', country = "AT", level = 2)
gadm_sub <- gadm[sample(1:length(gadm), 5), ]

plot(alt)
plot(gadm_sub, add=T)

asp <- terrain(alt, opt = "aspect", unit = "degrees", df = F)
slo <- terrain(alt, opt = "slope", unit = "degrees", df = F)

> extract(slo, gadm_sub, fun = mean, na.rm = T, small = T, df = T)
  ID     slope
1  1  9.959053
2  2  1.047443
3  3  7.456165
4  4  1.673786
5  5 11.946553

> extract(asp, gadm_sub, fun = mean, na.rm = T, small = T, df = T)
  ID   aspect
1  1 170.8065
2  2 184.0130
3  3 190.7155
4  4 136.8953
5  5 205.2115

Use R to Bulk-Download Digital Elevation Data with 1" Resolution

Here's a little r-script to convenientely download high quality digital elevation data, i.e. for the Alps, from HERE:

require(XML)

dir.create("D:/GIS_DataBase/DEM/")
setwd("D:/GIS_DataBase/DEM/")

doc <- htmlParse("http://www.viewfinderpanoramas.org/dem3.html#alps")
urls <- paste0("http://www.viewfinderpanoramas.org", xpathSApply(doc,'//*/a[contains(@href,"/dem1/N4")]/@href'))
names <- gsub(".*dem1/(\\w+\\.zip)", "\\1", urls)

for (i in 1:length(urls)) download.file(urls[i], names[i]) 

# unzip all files in dir and delete them afterwards
sapply(list.files(pattern = "*.zip"), unzip)
unlink(list.files(pattern = "*.zip"))

p.s.: Also check raster::getData which pulls SRTM data at 90m resolution for a location / region!

QGIS: Curing Small Aesthetical Flaw

Procrastination ahead! ..When starting QGIS, does the popping up of the cmd prompt window also annoy you like me? If you want to solve this, put the below vbs script in your PATH/bin folder (or anywhere else, if you wish).

Check the path to qgis.bat in the script and change it if yours is different. Then, go to the QGIS-Desktop shortcut and in the options dialogue point to the vbs script as target. Your done - no more popping cmd windows when starting QGIS!

Set WshShell = CreateObject("WScript.Shell")
WshShell.Run chr(34) & "C:\OSGeo4W\bin\qgis.bat" & Chr(34), 0
Set WshShell = Nothing

Creating a QGIS-Style (qml-file) with an R-Script

How to get from a txt-file with short names and labels to a QGIS-Style (qml-file)?
I used the below R-script to create a style for this legend table where I copy-pasted the parts I needed to a txt-file, like for the WRB-FULL (WRB-FULL: Full soil code of the STU from the World Reference Base for Soil Resources). The vector data to which I applied the style is freely available at ESDAC - you just need to submit a form to get access to the data. BTW, thanks to a helping hand on SO.

You can find the QGIS-styler script in theBioBucket-Repository on GitHub.

Rasterize CORINE Landcover 2006 Seamless Vector Data with OGR-rasterize

Following up my latest postings (HERE & HERE) on Corine Landcover I'll share how to rasterize this data at a desired resolution with OGR rasterize:

cd D:/GIS_DataBase/CorineLC/shps_app_and_extr/

# grab extracted and appended / merged file, created previously:
myfile=extr_and_app.dbf

name=${myfile%.dbf}

# rasterize y- and x-resolution = 50 map units: 
gdal_rasterize -a code_06 -tr 50 50 -l $name $myfile D:/GIS_DataBase/CorineLC/shps_app_and_extr/clc_tirol_raster_50.tif

Processing CORINE Land Cover 2006 Seamless Vector Data with OGR/GDAL's ogr2ogr and BASH

Lately I posted about using R to programmatically download all (43 zip-files) seamless 2006 vector data of CORINE Land Cover (see here). In this follow-up I share two bash-scripts with OGR/GDAL's ogr2ogr which I used to extract features by extent (Austria/Tyrol) and to combine all seperate shp-files into one by appending the dbfs. The style I used for the map is also available with german labels HERE.

##################################################
# name: extr_by_extent.sh
# Author: Kay Cichini
# Date: 30-04-2013
# Purpose: Extract shp-file features in {infolder} by input extent 
#          and save new shp-files to {outfolder} that is created on the fly
# Extent Tirol in EPSG:3035 - ETRS89 / ETRS-LAEA
# from D:\GIS_DataBase\GIS_Tirol\Tirol_Extent_LEAE-ETRS.shp
# xMin,yMin 4328054.73,2617730.23 : xMax,yMax 4547691.32,2739524.35

infolder=D:/GIS_DataBase/CorineLC/shps
outfolder=D:/GIS_DataBase/CorineLC/shps_extracted

ext='4328054.73 2617730.23 4547691.32 2739524.35'

ogr2ogr -overwrite $outfolder -spat $ext $infolder
##################################################

##################################################
# name: app_shps.sh
# Author: Kay Cichini
# Date: 30-04-2013
# Purpose: Append shp-files to newly created file

# working directory with shp-files to be appended into one file
mydir=D:/GIS_DataBase/CorineLC/shps_extracted
cd $mydir

# directory where final shp-file will be saved
mydir_final=D:/GIS_DataBase/CorineLC/shps_app_and_extr
mkdir $mydir_final

# get dbfs, which are the actual files to which append the data to
declare -a dbfs=(*.dbf)

# loop through dbfs in dir and append each to the dbf of 
# 'extr_and_app.shp' that will be created by ogr2ogr on the fly 
# and saved to {mydir_final}
for dbf in ${dbfs[@]}; do
  echo appending $dbf to $mydir_final/extr_and_app.dbf
  ogr2ogr -append $mydir_final/extr_and_app.dbf $dbf
done
##################################################

Programmatically Download CORINE Land Cover Seamless Vector Data with R

Thanks to a helpful SO-Answer I was able to download all CLC vector data (43 zip-files) programmatically:

require(XML)

path_to_files <- "D:/GIS_DataBase/CorineLC/Seamless"
dir.create(path_to_files)
setwd(path_to_files)

doc <- htmlParse("http://www.eea.europa.eu/data-and-maps/data/clc-2006-vector-data-version-2")
urls <- xpathSApply(doc,'//*/a[contains(@href,".zip/at_download/file")]/@href')

# function to get zip file names
get_zip_name <- function(x) unlist(strsplit(x, "/"))[grep(".zip", unlist(strsplit(x, "/")))]

# function to plug into sapply
dl_urls <- function(x) try(download.file(x, get_zip_name(x), mode = "wb"))

# download all zip-files
sapply(urls, dl_urls)

# function for unzipping
try_unzip <- function(x) try(unzip(x))

# unzip all files in dir and delete them afterwards
sapply(list.files(pattern = "*.zip"), try_unzip)

# unlink(list.files(pattern = "*.zip"))

Displaying any Google Maps Content in Google Earth via Network-Links

This is a handy option for displaying Google Maps in Google Earth. This could be any of your own maps created in Google Maps, a route or whatever. One nice feature of GE's Network-Links is that if you or any other contributor of a Google map modifies it (in Google Maps), it will be refreshed also in Google Earth.To do so, just grab the link from Google Maps. Like this one for a route:

https://maps.google.at/maps?saddr=Tivoli+Stadion+Tirol,+Innsbruck&daddr=Haggen,+Sankt+Sigmund+im+Sellrain&hl=de&ll=47.227496,11.250687&spn=0.12007,0.338173&sll=47.230526,11.250825&sspn=0.120063,0.338173&geocode=FagR0QIdYR-uACGQHB8W0fIkjCkfsEUhQ2mdRzGQHB8W0fIkjA%3BFZhY0AId6CypACnvxXqyRz2dRzG6da9r3hgbNA&oq=hagge&t=h&dirflg=h&mra=ls&z=12


Then go to GE and choose "Add" from main menu and "Network-Link" from the submenu. Then paste the above link and append "&output=kml" to it. That's it!

Myricaria Occurrence Map for Tyrol, Austria

This is a map of the current occurrence data of Myricaria germanica (courtesey of the Tiroler Landesmuseum Ferdinandeum) created with the freeware QGIS.

Calculate Single Contour-Line from DEM with QGIS / GDAL

In QGIS:

- from menu: Raster / Extraction / Contour

- define output name path/to/output.shp

- alter GDAL call for single contour line at 900 m asl:

gdal_contour -fl 900 "path/to/dem_raster.asc" "path/to/output.shp"

- for removing small poplygons or lines add area or length field (attr table / field calc or vector / geometry / add geometry)

- query by length or area to deselect unwanted iso-lines


Finally, you can export the contours as KML and check it in Google Earth:

Use IUCN API with R & XPath

Thanks to a posting on R-sig-eco mailing list I learned of the IUCN-API. Here's a simple example for what can be done with it (output as pdf is HERE):

Line Slope Calculation in ArcGis 9.3 (using XTools)

Ojective:
You have a polyline, say a path, river, etc., and want to know average slope of each single line.

Approach:
Use z-values of nodes of polyline and calculate percentual slope by
(line-segment shape_length / z-Difference) * 100

..A Quick Geo-Trick for GoogleMaps in R (using dismo)

... I thought this geocoding-bit might be worth sharing (found HERE when searching the web for dismo-documentation).

Convert Date Field into Year in ArcGis

Objective: a date field [date] with format dd.mm.yyyy, i.e.,  should be converted to format yyyy.

Solution:

(1) Add an integer field [YEAR] to the attribute table.

(2) Field Calculator: YEAR =

               year([date])

Using IUCN-Data, ArcMap 9.3 and R to Map Species Diversity

I'm overwhelmed by the ever-growing loads of good data that's made available via the web. For instance IUCN collects and hosts spatial species data which is free for download. I'm itching to play with all this data and in the end there may arise some valuable outcome:)

In the below example I made a map for amphibian species richness without much effort. I used a 5x5 km grid (provided by Statistik Austria) and amphibian species ranges and intersected this data to yield species numbers per grid-cell. For aggregation of the tables from the polygons, that were produced by intersection, I needed to switch from ArcMap to R because ESRI's "Summary Statistics" only cover MIN, MAX, MEAN, FIRST/LAST... and there is no easy way to apply custom functions. So I just accessed the dbf that I wanted to aggregate with R and did the calculations there.. Then reloaded the layer in ArcMap and was done..

A Simple Example for the Use of Shapefiles in R

A simple example for drawing an occurrence-map (polygons with species' points) with the R-packages maptools and sp using shapefiles.
HERE is the example data.

Nice Species Distribution Maps with GBIF-Data in R

Here's an example of how to easily produce real nice distribution maps from GBIF-data in R with package maps...

ArcGis Style with Marker Symbols for Tetrao tetrix & T. urogallus

...lately I designed two nice symbol markers for Tetrao spp. in vector fromat - they are for free download here.

Picture Marker Symbols in ArcGis with Vector Format

Ever wanted to add pictures as marker symbols in ArcGis? ..then you may have noticed that raster images look poor in an ArcGis map no matter how high the resolution of the input image is...

Neighborhood-Statistics for Samplepoints Based on Raster-Data

There's no ArcGis-tool for calculating neighborhood-statistics for samplepoints based on raster-data: say you had a raster with presence/absence data of a given feature and for a set of samplepoints you wanted to know the feature's occurrence rate (percentage of cells with presence of that feature) in the samplepoints neighborhood, then you will find this model (saved in toolbox NESTPTS.tbx) useful.

Visual Basic - IF THEN Usage in ArcGis Field Calculator

Recently I had to do some statistics on a numerical field in an attribute table and I had to classify this field according to a set of intervals for this purpose. I had to learn that there is no tool for this. So I had to use the fieldcalculator with the appropiate IF THEN VBA-Syntax:

Here is how it can be done (one of many ways):