Online Trick

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

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

I'll use the below script that retrieves data for threatened bird species from the World Bank via its API and does some processing, plotting and analysis. There is a package (WDI) that allows you to access the data easily.

# world bank indicators for species - 
# I'll check bird species:
code <- as.character(WDIsearch("bird")[1,1])
bird_data <- WDI(country="all", indicator=code, start=2010, end=2012)

# remove NAs and select values in the range 50 - 1000:
bird_data_sub <- bird_data[!is.na(bird_data$EN.BIR.THRD.NO)&
                           bird_data$EN.BIR.THRD.NO < 1000&
                           bird_data$EN.BIR.THRD.NO > 50, ]

# change in numbers across years 2010 and 2011:
change.no <- aggregate(EN.BIR.THRD.NO ~ country, diff,
                       data = bird_data_sub)
# plot:
par(mar = c(3, 3, 5, 1))
plot(x = change.no[,2], y = 1:nrow(change.no),
     xlim = c(-12, 12), xlab = "", ylab = "",
     yaxt = "n")
abline(v = 0, lty = 2, col = "grey80")
title(main = "Change in Threatened Bird Species in\nCountries with Rich Avifauna (>50)")
text(y = 1:nrow(change.no), 
     x = -2, adj = 1,
     labels = change.no$country)
segments(x0 = 0, y0 = 1:nrow(change.no),
         x1 = change.no[, 2], y1 =  1:nrow(change.no))

# test hypothesis that probability of species decrease is
# equal to probability of increase:
binom.test(sum(change.no < 0), sum(change.no != 0))

For generating the report you can source the script from dropbox.com and stitch it in this fashion:

stitch("http://dl.dropbox.com/s/ga0qbk1o17n17jj/Change_threatened_species.R")

..this is one line of code - can you dig it?..
BTW, for simplicity I use knitr::stitch with its default template...

You should get something like THIS PDF.

EDIT, MARCH 2013
OUTDATED! you can use this approach instead:

library(knitr); library(RCurl); library(WDI)

destfile = "script.txt"
x = getBinaryURL("https://dl.dropbox.com/s/ga0qbk1o17n17jj/Change_threatened_species.R", followlocation = TRUE, ssl.verifypeer = FALSE)
writeBin(x, destfile, useBytes = TRUE)
source(paste(tempdir(), "/script.txt", sep = ""))

stitch(paste(tempdir(), "/script.txt", sep = ""))

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

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://onlinetrickpdf.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) 

You can query global climate data from the KNMI Climate Explorer (the KNMI is the Royal Netherlands Metereological Institute) with R.


Here's a little example how I retreived data for my hometown Innsbruck, Austria and plotted annual total precipitation. You can choose station data by pointing at a map, by setting coordinates, etc.

# get climate (precipitation) data from url:
# http://climexp.knmi.nl/selectstation.cgi?id=someone@somewhere

# station INNSBRUCK, FLUGHAFEN (11120), 47.27N, 11.35E:
ibk_dat <- read.table("http://climexp.knmi.nl/data/pa11120.dat", sep = "", 
                      row.names = 1, col.names = 0:12)

# cut off first and last yr, due to missing data..
ibk_dat <- ibk_dat[c(-1, -50,]

# plot yearly sums:
windows(width = 15, height = 5)
plot(rowSums(ibk_dat), type = "s", ylab = "Annual Total Precipitation (mm)", 
     xlab = NA, col = "blue", xaxt = "n", lwd = 1.5, las = 2, cex.axis = 0.8,
     main = "INNSBRUCK FLUGHAFEN, 47.27N, 11.35E, 593m, WMO station code: 11120") 
axis(1, labels = rownames(ibk_dat), at = 1:nrow(ibk_dat), las = 2, cex.axis = 0.85)

abline(h = mean(rowSums(ibk_dat)), col = 1, lty = 2, lwd = 1.2)
text(1250, adj = 0, "Long-term average", cex = 0.75)
arrows(x0 = 2.5, y0 = 1220,
       x1 = 2.5, y1 = 930, length = 0.05)