Dr Philipp Boersch‑Supan

Marine & Quantitative Ecologist

Mapping my trip with GMT and Inkscape

This morning I tweeted some maps of my upcoming trip to sub-antarctic New Zealand. Here is how I made them.

For most of my dealings with spatial data I use either R (with it’s great ecosystem of ), or a full-fledged GIS package such as GRASS or QGIS (if I need a more interactive work environment and/or particular geoinformatics tools). This combination covers pretty much all the bases for my spatial analysis and mapping needs. However, occasionally I want to make a pretty map that isn’t rectangular (usually world maps, or maps of large sections of the world) and/or one that is projected in a certain way. Maybe something like this figure of the circulation in the Southern Indian Ocean I made for my PhD thesis:

SWIO circulation

Making “publication quality” maps in all kinds of projections is a real strength of the Generic Mapping Tools aka GMT, a software suite that is popular with marine geologists and is firmly grounded in the unix paradigm of one tool - one command line utility. As such the syntax and the mapping workflow may seem like a blast from the past for many (cool kids use D3 for fancy maps these days, but apparently I’m not one of them) - but with a moderate amount of pain one can get beautiful maps that put your old Mercator or unprojected WGS84 rectangles to shame.

GMT natively writes all output to postscript files, which makes the resulting maps editable in vector graphics software such as Inkscape. This is great for tweaking or polishing graphics that go beyond maps of spatial data. In the above figure of ocean currents, for example, I used GMT to plot the globe and the bathymetry (subsea topography), but the stylized fronts and currents were added in Inkscape (based on projected coordinates of frontal locations plotted with GMT in an intermediate step).

I did something similar, but much simpler, for the maps I tweeted today, and I’ll show the individual steps below.

Step 1: Draw the map

I adapted one of the GMT examples to show what’s on the opposite side of the world during my trip. The code for this is given below, and also in the script antipodes.sh on github.

This time I was happy with the GMT output as is, so no post-GMT modifications.

#!/bin/bash
#	  This map is a quick and dirty adaptation of GMT EXAMPLE 25
#   (using GMT 5.4.0_r17345)
#
# Purpose:	Display distribution of antipode types around NZ
# GMT modules:	gmtset, grdlandmask, grdmath, grd2xyz, gmtmath, grdimage, pscoast, pslegend
# Unix progs:	cat
#
# Create D minutes global grid with -1 over oceans and +1 over land
D=10
gmt grdlandmask -Rg -I${D}m -Dc -A500 -N-1/1/1/1/1 -r -Gwetdry.nc
# Manipulate so -1 means ocean/ocean antipode, +1 = land/land, and 0 elsewhere
gmt grdmath -fg wetdry.nc DUP 180 ROTX FLIPUD ADD 2 DIV = key.nc
# Calculate percentage area of each type of antipode match.
gmt grdmath -Rg -I${D}m -r Y COSD 60 $D DIV 360 MUL DUP MUL PI DIV DIV 100 MUL = scale.nc
gmt grdmath -fg key.nc -1 EQ 0 NAN scale.nc MUL = tmp.nc
gmt grd2xyz tmp.nc -s -ZTLf > key.b
ocean=`gmt math -bi1f -Ca -S key.b SUM UPPER RINT =`
gmt grdmath -fg key.nc 1 EQ 0 NAN scale.nc MUL = tmp.nc
gmt grd2xyz tmp.nc -s -ZTLf > key.b
land=`gmt math -bi1f -Ca -S key.b SUM UPPER RINT =`
gmt grdmath -fg key.nc 0 EQ 0 NAN scale.nc MUL = tmp.nc
gmt grd2xyz tmp.nc -s -ZTLf > key.b
mixed=`gmt math -bi1f -Ca -S key.b SUM UPPER RINT =`
# Generate corresponding color table
gmt makecpt -Cblue,gray,red -T-1.5/1.5/1 -N > key.cpt
# Create the final plot and overlay coastlines
gmt set FONT_ANNOT_PRIMARY +10p FORMAT_GEO_MAP dddF
range="g"
proj="G-185/-45/4.5i"
psfile='enderby_antipodes.ps'
#plot antipodes raster
gmt grdimage key.nc -R${range} -J${proj} -Ckey.cpt -B10/10 -P -K > $psfile
#overplot landcover
gmt pscoast -R${range} -J${proj} -Dl -A0/0/1 -G35 -W0.4p -O -K >> $psfile
#add legend
gmt pslegend -R${range} -J${proj} -O -DJBC+w4.5i -Y-0.75i >> $psfile << END
N 2
S 0.15i s 0.2i gray  0.25p 0.3i Terrestrial Antipodes
S 0.15i s 0.2i blue 0.25p 0.3i Oceanic Antipodes 
END
#convert to PNG
gmt psconvert -Tg -A $psfile
#remove intermediate results
rm -f wetdry.nc scale.nc key.nc tmp.nc key.cpt key.b gmt.conf

antipodes map


Enderby Trust Scholarship Voyage

Map of New Zealand and travel route

I am beyond thrilled to have been awarded a travel scholarship by the Enderby Trust, to go on a trip to the sub-Antarctic islands of New Zealand. In mid-December I will be visiting Tini Heke (Snares Island), Motu Maha (Auckland Islands), and Motu Ihupuku (Campbell Island), and – weather allowing – some of the seabird colonies on them. I am very excited to get this chance to see the species I am modelling in their natural habitats.

I will not be able to update the blog during the expedition, as there is no regular internet connection on the vessel. The vessel position and cccasional updates from the expedition leader can be found at http://www.heritage-expeditions.com/captains-blog/.


Time series similiarity search: New paper and R package rucrdtw

CRAN_Status_Badge Peer Reviewed Software status

I am very happy to announce a new software paper and R package published today in The Journal of Open Source Software and on CRAN. rucrdtw is a wrapper for the UCR Suite, a C++ library for very fast time series similarity searches under dynamic time warping (DTW).

Dynamic Time Warping (DTW) methods provide algorithms to optimally map a given time series onto all or part of another time series (Berndt and Clifford 1994). The remaining cumulative distance between the series after the alignement is a useful distance metric in time series data mining applications for tasks such as classification, clustering, and anomaly detection.

Calculating a DTW alignment is computationally relatively expensive, and as a consequence DTW is often a bottleneck in time series data mining applications. The UCR Suite (Rakthanmanon et al. 2012) provides a highly optimized algorithm for best-match subsequence searches that avoids unnecessary distance computations and thereby enables fast DTW and Euclidean Distance queries even in data sets containing trillions of observations.

Figure 1: UCR DTW is approximately 3 orders of magnitude faster than a naive sliding-window search using DTW distance.

A broad suite of DTW algorithms is implemented in R in the dtw package (Giorgino 2009). The rucrdtw R package provides complementary functionality for fast similarity searches by providing R bindings for the UCR Suite via Rcpp (Eddelbuettel and Francois 2011). In addition to queries and data stored in text files, rucrdtw also implements methods for queries and/or data that are held in memory as R objects, as well as a method to do fast similarity searches against reference libraries of time series.

References

Berndt, Donald J, and James Clifford. 1994. “Using Dynamic Time Warping to Find Patterns in Time Series.” In KDD Workshop, 10:359–70. 16. AAAI. http://www.aaai.org/Library/Workshops/1994/ws94-03-031.php.

Eddelbuettel, Dirk, and Romain Francois. 2011. “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software 40 (1): 1–18. doi:10.18637/jss.v040.i08.

Giorgino, Toni. 2009. “Computing and Visualizing Dynamic Time Warping Alignments in R: The Dtw Package.” Journal of Statistical Software 31 (7): 1–24. doi:10.18637/jss.v031.i07.

Rakthanmanon, Thanawin, Bilson Campana, Abdullah Mueen, Gustavo Batista, Brandon Westover, Qiang Zhu, Jesin Zakaria, and Eamonn Keogh. 2012. “Searching and Mining Trillions of Time Series Subsequences Under Dynamic Time Warping.” In Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 262–70. ACM. doi:10.1145/2339530.2339576.


deBInfer v0.4.1 now on CRAN!

CRAN Badge

We are very excited to announce the initial CRAN release of our R package for Bayesian inference in differential equation models. At this point we consider the package API stable, and we will strive to preserve backwards compatibility in future releases. Also included in the CRAN release are three vignettes covering example applications for ODEs and DDEs.

The CRAN version of the package can be found at https://cran.r-project.org/package=deBInfer and the package can be installed in R with thecommand install.packages("deBInfer").

We continue to work on the package, and a development version is available at https://github.com/pboesu/debinfer it can be installed directly from github using the devtools package: install.packages("devtools") devtools::install_github("pboesu/debinfer")

HTML help files for the development version, and the most recent versions of the vignettes are available at http://pboesu.github.io/debinfer

Background on the methods implemented in deBInfer can be found in our paper deBInfer: Bayesian inference for dynamical models of biological systems in R


Conference summer part 1: ISEC and ISBA presentations

We have been busy spreading the word about our new R package for implementing Bayesian inference for differential equation models. Our poster from the 2016 meeting of the International Society of Bayesian Analysis and slides from Leah’s talk at the International Statistical Ecology Conference are now available on figshare.  For more detail have a look at the arXiv preprint, or better yet, download the package from github, and try it.

At ISEC we also presented a poster featuring some early results on our work on albatross movement models, and their identifiability when using data from different generations of biologgers.