16 Sept 2013

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

10 Sept 2013

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

26 Aug 2013

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.

19 Aug 2013

Text Mining with R - Comparing Word Counts in Two Text Documents

Here's what I came up with to compare word counts in two pieces of text. If you got any idea, I'd love to learn about alternatives!

## a function that compares word counts in two texts
wordcount <- function(x, y, stem = F, minlen = 1, marg = F) {

                        require(tm)

                        x_clean <- unlist(strsplit(removePunctuation(x), "\\s+"))
                        y_clean <- unlist(strsplit(removePunctuation(y), "\\s+"))

                        x_clean <- tolower(x_clean[nchar(x_clean) >= minlen])
                        y_clean <- tolower(y_clean[nchar(y_clean) >= minlen])

                        if ( stem == T ) {

                          x_stem <- stemDocument(x_clean)
                          y_stem <- stemDocument(y_clean)
                          x_tab <- table(x_stem)
                          y_tab <- table(y_stem)    

                          cnam <- sort(unique(c(names(x_tab), names(y_tab))))

                          z <- matrix(rep(0, 3*(length(cnam)+1)), 3, length(cnam)+1, dimnames=list(c("x", "y", "rowsum"), c(cnam, "colsum")))
                          z["x", names(x_tab)] <- x_tab
                          z["y", names(y_tab)] <- y_tab
                          z["rowsum",] <- colSums(z)
                          z[,"colsum"] <- rowSums(z)
                          ifelse(marg == T, return(t(z)), return(t(z[1:dim(z)[1]-1, 1:dim(z)[2]-1])))

                          } else { 

                          x_tab <- table(x_clean)
                          y_tab <- table(y_clean)    

                          cnam <- sort(unique(c(names(x_tab), names(y_tab))))

                          z <- matrix(rep(0, 3*(length(cnam)+1)), 3, length(cnam)+1, dimnames=list(c("x", "y", "rowsum"), c(cnam, "colsum")))
                          z["x", names(x_tab)] <- x_tab
                          z["y", names(y_tab)] <- y_tab
                          z["rowsum",] <- colSums(z)
                          z[,"colsum"] <- rowSums(z)
                          ifelse(marg == T, return(t(z)), return(t(z[1:dim(z)[1]-1, 1:dim(z)[2]-1])))
                          }
                        }

## example
x = "Hello new, new world, this is one of my nice text documents - I wrote it today"
y = "Good bye old, old world, this is a nicely and well written text document"

wordcount(x, y, stem = T, minlen = 3, marg = T)

Follow-Up:

Thanks a lot for the comments! As I'm not that much into text mining I was trying to reinvent the wheel (in a rather dilettante manner) - missing the capabilities of existing packages. Here's the shortest code that I was able to find doing the same thing (with the potential to get out much more of it, if desired).
x = "Hello new, new world, this is one of my nice text documents"
y = "Good bye old, old world, this is a text document"
z = "Good bye old, old world, this is a text document with WORDS for STEMMING  - BTW, what is the stem of irregular verbs like write, wrote, written?"

# make a corpus with two or more documents (the cool thing here is that it could be endless (almost) numbers 
# of documents to be cross tabulated with the used terms. And the control function enables you
# to do lots of tricks with it before it will be tabulated, see ?termFreq, i.e.)

xyz <- as.list(c(x,y,z))
xyz_corp <- Corpus(VectorSource(xyz))

cntr <- list(removePunctuation = T, stemming = T, wordLengths = c(3, Inf))

as.matrix(TermDocumentMatrix(xyz_corp, control = cntr))

20 Jun 2013

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)

18 Jun 2013

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!

7 Jun 2013

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

21 May 2013

R Quick Tip: Shutdown Windows after Script Has Finished

Quite often I have long procedures running and want to do this over night. However, my computer would still be running all night after the script has finished. This is easily circumvented by the following lines that I put at the end of such a script:

# set working dir
# setwd("C:/Users/Kay/Desktop")
 
# long procedure:
for(i in 1:1e+5) {cat(i); cat("\n..................\n")}
 
d <- "something"
 
# save history
savehistory()
 
# and worspace
save.image()
 
# then shutdown after 240 s
system("C:/Windows/system32/shutdown.exe -f -s -t 240")
 
# this would abort the shutdown:
# system("C:/Windows/system32/shutdown.exe -a")

6 May 2013

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.

2 May 2013

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

21 Apr 2013

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

18 Apr 2013

Use ZOTERO Bibliographies with Custom CSL Citation Styles in MS Word Processor

Just discovered this nice online CSL editor which is very handy for plugging together your desired citation styles. With Zotero's Word Plug-In for Firefox (Tutorial) you can nicely insert intext citation and bibliographies in an existing or tailored citation style.

12 Apr 2013

Download File from Google Drive/Docs Programmatically with R

Following up my lattest posting on how to download files from the cloud with R..

dl_from_GoogleD <- function(output, key, format) {

## Arguments:
## output = output file name
## key = Google document key
## format = output format (pdf, rtf, doc, txt..)
## Note: File must be shareable!

                        require(RCurl)
                        bin <- getBinaryURL(paste0("https://docs.google.com/document/d/", key, "/export?format=", format),
                                            ssl.verifypeer = FALSE)
                        con <- file(output, open = "wb")
                        writeBin(bin, con)
                        close(con)
                        message(noquote(paste(output, "read into", getwd())))                        
                        }


# Example:
dl_from_GoogleD(output = "dl_test.pdf", 
                key = "1DdauvkcVm5XtRBkQIv1na8PeLAwpCBdW8pALCFpRWeM",
                format = "pdf")
shell.exec("dl_test.pdf")
EDIT: Here's how it can be done for spreadsheet-like data, like HERE, which is a comma seperated file with .txt extension saved to Google Drive. See also this post
library(RCurl)
setwd(tempdir())
destfile = "test_google_docs.csv"
x = getBinaryURL("https://docs.google.com/uc?export=download&id=0B2wAunwURQNsR0I0a0NlQUlJdzA", followlocation = TRUE, ssl.verifypeer = FALSE)
writeBin(x, destfile, useBytes = TRUE)
shell.exec(paste(tempdir(), "/test_google_docs.csv", sep = ""))

10 Apr 2013

Tweaking Movie Subtitles with R

I use R to fix subtitles that are not in sync with my movies. For the example below the subs were showing too early - so I added some time to each sequence in the srt file. For simplicity I used exactly 1 second in the below example.
You'll see that I use my function dl_from_dropbox(), on which I wrote a post previously, to get the example file!



Download Files from Dropbox Programmatically with R

Here is a usefull snippet that I stole from qdap::url_dl to download files from my Dropbox to the working directory.
Argument x is the document name and d the document key.

dl_from_dropbox <- function(x, key) {
                        require(RCurl)
                        bin <- getBinaryURL(paste0("https://dl.dropboxusercontent.com/s/", key, "/", x),
                                            ssl.verifypeer = FALSE)
                        con <- file(x, open = "wb")
                        writeBin(bin, con)
                        close(con)
                        message(noquote(paste(x, "read into", getwd())))                        
                        }
# Example:
dl_from_dropbox("GViewer_Embeds.txt", "06fqlz6gswj80nj")
shell.exec("GViewer_Embeds.txt")
PS: Also see this R-package for interfacing with Dropbox

5 Apr 2013

Embedding Google Docs and Google Viewer

..Here's a short example of the usage of the Google Viewer and Google Docs embeds (see here). If you want to embed a view of any file on Google Drive/Docs you can do this by putting a pdf preview (the file could be in any format, see first example, which is a Google text document) in an iframe or, in the case of spreadsheets, the speradsheet itself into an iframe.

The second embed below shows the latest outcome of an online questionnaire and in the speradsheet that holds the result the new answers to the questoinnaire are updated every 5 minutes after refreshing the view. So, you can, i.e., make an online questoinnaire with a Google form and display the results like below. With a little tweaking of the parameters, which I obviously haven't done yet, the result would of course be more visually appealing.

If you try for one of your files you just need to replace the id/key in the URL..





The Result:



2 Apr 2013

VBS Shutdown Timer

When running time comsuming processes I often wish to shutdown automatically after a given time. There are several applications on the internet but I was not really happy with the ones I found. So, I put together a really small VB-Script which in fact does the job perfectly.

Dim Shell : Set Shell = CreateObject("WScript.Shell")

argument_t = inputbox("Time in Seconds:", "Shutdown", "3600")

Shell.Run "C:/Windows/system32/shutdown.exe -f -s -t " & _
          argument_t, 0

MsgBox "Press OK if you want to cancel the Shutdown..", _
       vbInformation + vbOkayonly, "Cancel Shutdown"

Shell.Run "C:/Windows/system32/shutdown.exe -a", 0

I put the vbs file on my disk and a shortcut to it on the Desktop. You can choose the time until shutdown with an input box.



14 Mar 2013

Example Code for a Customized RSS-Feed

A short code snippet for the customized RSS-Feed on this Google-Blogger site. I simply added the below code to a gadget to yield a tailored RSS-Feed in the sidebar. Most of the code was grabbed from the Google Code Playground.



Land Tirol





    
 

Loading...

11 Mar 2013

MS Word 2010, Find & Replace with Wildcards

Quite often I need to re-format citations in my texts, i.e., missing brackets around the year of a publication. In MS Word I simlpy solve such issues with the find and replace option using wildcards. I.e. for the issue with years without brackets:

short cut crtl-h
then put ([1-2][0-9]{3}) in the find-box and
(^&) in the replace box
hit 'replace' and you're done

The find string will search for any 4-digit number starting with 1 or 2 followed by 3 digits between 0 and 9 and will be replaced by the same number put inbetween round brackets. Note that the brackets in the find-string denote a token which then can be referenced by the string ^& in the replace step.

15 Feb 2013

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!

3 Feb 2013

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.