• gotchas abound
• there's more than one way to do it

• plotting data on a map
• controlling plot region and orientation
• get better map data
• scaling plots by data value
• GDAL raster manipulation
• build the same map in QGIS

f <- "SAZSENSE_TM_20120209.txt"
szt <- read.table(f, comment.char = "/", sep = "\t", header = TRUE)
summary(szt)

##     Cruise       Station   Type                  date         Lon       
##  AA0307:148   10     :12   B:148   0007-01-20T16:53:12   Min.   :140.5  
##               12     :12           0007-01-22T18:25:12   1st Qu.:145.9  
##               2b     :12           0007-01-29T23:15:12   Median :149.4  
##               4a     :12           0007-02-01T08:53:12   Mean   :148.5  
##               5      :12           0007-02-08T06:52:12   3rd Qu.:152.5  
##               6a     :12           0007-02-09T00:45:12   Max.   :153.7  
##               (Other):76           (Other)         :76                  
##       Lat            BotDepth          Depth              Al        
##  Min.   :-53.99   Min.   : 250.0   Min.   :  15.0   Min.   :0.0970  
##  1st Qu.:-49.00   1st Qu.: 500.0   1st Qu.:  50.0   1st Qu.:0.4233  
##  Median :-45.58   Median :1000.0   Median : 125.0   Median :0.8215  
##  Mean   :-47.22   Mean   : 764.9   Mean   : 225.6   Mean   :1.6875  
##  3rd Qu.:-44.95   3rd Qu.:1000.0   3rd Qu.: 262.5   3rd Qu.:2.3213  
##  Max.   :-43.69   Max.   :1000.0   Max.   :1000.0   Max.   :7.3430  
##                                                     NA's   :2       
##        Co              Ni              Cu               Zn         
##  Min.   : 5.20   Min.   :2.510   Min.   :0.2170   Min.   : 0.1000  
##  1st Qu.:15.60   1st Qu.:3.527   1st Qu.:0.4622   1st Qu.: 0.9275  
##  Median :25.35   Median :4.045   Median :0.5430   Median : 1.9600  
##  Mean   :23.40   Mean   :4.258   Mean   :0.6068   Mean   : 3.2236  
##  3rd Qu.:29.80   3rd Qu.:4.718   3rd Qu.:0.6432   3rd Qu.: 4.1150  
##  Max.   :42.90   Max.   :7.310   Max.   :1.7930   Max.   :17.7400  
##                                                                    
##        Cd               Pb               Fe        
##  Min.   : -6.20   Min.   :  7.40   Min.   :0.0900  
##  1st Qu.: 47.85   1st Qu.: 17.82   1st Qu.:0.1900  
##  Median :164.10   Median : 23.15   Median :0.2700  
##  Mean   :198.23   Mean   : 29.00   Mean   :0.2849  
##  3rd Qu.:247.88   3rd Qu.: 32.90   3rd Qu.:0.3600  
##  Max.   :900.40   Max.   :128.40   Max.   :0.6900  
## 

plot(Lat ~ Lon, data = szt)

library(maptools)
data(wrld_simpl)
plot(wrld_simpl)
points(Lat ~ Lon, data = szt)

plot(wrld_simpl, xlim = range(szt$Lon) + c(-5, 5), 
                  ylim = range(szt$Lat) + c(-3, 3))
points(Lat ~ Lon, data = szt)

plot(Lat ~ Lon, data = szt, xlim = range(szt$Lon) + c(-5, 5), 
     ylim = range(szt$Lat) + c(-12, 5), axes = FALSE)
plot(wrld_simpl, add = TRUE)
degAxis(1);degAxis(2);box()

plot(subset(wrld_simpl, NAME == "Australia"), asp = "")

plot(subset(wrld_simpl, NAME == "Australia"), asp = 1)

plot(subset(wrld_simpl, NAME == "Australia"), asp = 1/cos(-30 * pi/180))

xx <- c(140, 155); yy <- c(-45, -20)
plot(wrld_simpl, xlim = xx, ylim = yy, asp = 1.5); axis(1);axis(2)
rect(xx[1], yy[1], xx[2], yy[2], border = "dodgerblue", lwd = 2)

source("aspectplot.R") ## thanks Spoon
op <- aspectplot(xlim = xx, ylim = yy, asp = 1.5)
plot(wrld_simpl, add = TRUE); axis(1);axis(2)
rect(xx[1], yy[1], xx[2], yy[2], border = "dodgerblue", lwd = 2)

par(op)

plot(wrld_simpl, xlim = c(100, 160), ylim = c(-85, -40), asp = "")
axis(1);axis(2)

op <- aspectplot(xlim = c(100, 160), ylim = c(-85, -40), asp = 1/cos(-60 * pi/180))
plot(wrld_simpl, add = TRUE); axis(1);axis(2)

par(op)

library(maptools)
data(wrld_simpl)
cmatch <- c("Australia", "New Zealand", "Antarctica")
aanz <- subset(wrld_simpl, NAME %in% cmatch)
plot(aanz)

library(raster)
expoly <- as(extent(100, 175, -78, -30), "SpatialPolygons")
library(rgeos)
aanzclip2 <- gIntersection(aanz, expoly, byid = TRUE)
plot(aanzclip2, col = c("red", "green", "blue"))
title("worst colour scheme ever")

gIntersection(aanz, expoly, byid = TRUE)

## class       : SpatialPolygons 
## features    : 3 
## extent      : 100, 175, -78, -30  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

gIntersection(aanz, expoly, byid = FALSE)

## class       : SpatialPolygons 
## features    : 1 
## extent      : 100, 175, -78, -30  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

Raster has a booty function:

library(raster)
x0 <- getData("GADM", country = "AUS", level = 0)  ## 6.8  Mb
plot(wrld_simpl, xlim = c(147.96, 148.02), 
     ylim = c(-43.26, -43.20), col = "grey")
plot(x0, add = TRUE); title("wrld_simpl vs GADM level 0")

Better than wrld_simpl, but juggling.

oz <- getData("GADM", country = "AUS", level = 0)  
nz <- getData("GADM", country = "NZ", level = 0)  
ant <- getData("GADM", country = "ATA", level = 0)

library(raster)
help(getData)
library(rgdal)
help(readOGR)

GDAL can read vector data from mostly anywhere, in R and QGIS.

Vector drivers: ^{_{ARCGEN, AVCBin, AVCE00, AeronavFAA, BNA, CSV, CartoDB, CouchDB, DGN, DXF, EDIGEO, ESRI Shapefile, ElasticSearch, GFT, GME, GML, GMT, GPKG, GPSBabel, GPSTrackMaker, GPX, GeoJSON, GeoRSS, Geoconcept, Geomedia, HTF, Idrisi, Interlis 1, Interlis 2, KML, MSSQLSpatial, MapInfo File, Memory, NAS, ODBC, ODS, OSM, OpenAir, OpenFileGDB, PCIDSK, PDF, PDS, PGDump, PGeo, PostgreSQL, REC, S57, SDTS, SEGUKOOA, SEGY, SQLite, SUA, SVG, SXF, TIGER, UK .NTF, VFK, VRT, WAsP, WFS, Walk, XLSX, XPlane}}

scl <- function(x) (x - min(x, na.rm = TRUE))/diff(range(x, na.rm = TRUE))
plot(Lat ~ Lon, data = szt, cex = scl(sqrt(Al)) * 2 + 1)

plot(Al ~ Depth, data = szt, cex = scl(sqrt(Al)) * 2 + 1)

library(rgdal)
prj <- "+proj=laea +lat_0=-52 +lon_0=147 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
oz2 <- spTransform(oz, CRS(prj))
nz2 <- spTransform(nz, CRS(prj))
ant2 <- spTransform(ant, CRS(prj))
xy <- project(cbind(szt$Lon, szt$Lat), prj)
ex <- 0.5 * union(extent(oz2), union(extent(nz2), extent(ant2)))
plot(ex, type = "n", asp = 1, axes = FALSE, xlab = "", ylab = "")
plot(oz2, add = TRUE, col = "grey", border = "transparent")
plot(nz2, add = TRUE, col = "firebrick", border = "transparent")
plot(ant2, add = TRUE, col = "dodgerblue", border = "transparent")
points(xy, cex = scl(sqrt(szt$Al)) * 1.5 + 1)

• Use GDAL to regrid Etopo2 for background in efficient format
• Import data and maps to QGIS
• Format data: data-scaled point colour and size
• Forget aspect ratio: use a map projection!
  
• Print composer

gdal_translate ETOPO2v2c_f4.nc etopo_sh_015.tif \
 -a_srs "+proj=longlat +datum=WGS84" \ 
 -projwin -180 0 180 -90 \
 -outsize 20% 20%  \ 
 -co TILED=YES -co COMPRESS=LZW

• open gdalinfo screendump

alt text

• start QGIS

+proj=laea +lat_0=-52 +lon_0=147 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

alt text

World outlines source: http://www.gadm.org/

QGIS: www.qgis.org

GDAL: www.gdal.org

OSGEO installer for Windows: http://trac.osgeo.org/osgeo4w/

Better R data and graphics: dplyr and ggvis

• GIS data structures
• Cropping / editing polygons
• Problems with Antarctica world shape
• Rasterizing polygons and lines
• Coercion between classes
• Interpolating from irregular data (Y~X, Depth~Lat, etc.)
• Tools for tracking data
• Calculating distance, distance-to-stuff
• Projections - sinusoidal, circumpolar, conformal conic, …
• Drawing graticules (lon-lat lines)

mdsumner@gmail.com Australian Antarctic Division