Let's draw a map

require(sp)
require(rgdal)
scot = readOGR("./Data", "scot")

## OGR data source with driver: ESRI Shapefile 
## Source: "./Data", layer: "scot"
## with 56 features and 5 fields
## Feature type: wkbPolygon with 2 dimensions

plot(scot)

What have we got?

class(scot)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

summary(scot)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##      min     max
## x   7151  468393
## y 529557 1218479
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs]
## Data attributes:
##            NAME          ID           Cases          Expected    
##  Aberdeen    : 1   Min.   : 1.0   Min.   : 0.00   Min.   : 1.10  
##  Angus       : 1   1st Qu.:14.8   1st Qu.: 4.75   1st Qu.: 4.05  
##  Annandale   : 1   Median :28.5   Median : 8.00   Median : 6.30  
##  Argyll-Bute : 1   Mean   :28.5   Mean   : 9.57   Mean   : 9.57  
##  Banff-Buchan: 1   3rd Qu.:42.2   3rd Qu.:11.00   3rd Qu.:10.12  
##  Bearsden    : 1   Max.   :56.0   Max.   :39.00   Max.   :88.70  
##  (Other)     :50                                                 
##       AFF       
##  Min.   : 0.00  
##  1st Qu.: 1.00  
##  Median : 7.00  
##  Mean   : 8.66  
##  3rd Qu.:11.50  
##  Max.   :24.00  
##

Features

par(mfrow = c(2, 2))
hist(scot$Cases, main = "Cases")
hist(scot$Cases/scot$Expected, main = "Rate")
plot(scot[2, ])
title(paste(scot$NAME[2], "\n", scot$Cases[2], " cases"))
plot(scot)
plot(scot[scot$Cases > 20, ], add = TRUE, col = "red")

par(mfrow = c(1, 1))

Colouring

cases = scot$Cases
scale10 = 1 + 9 * (cases - min(cases))/(diff(range(cases)))
require(RColorBrewer)
display.brewer.all()

palette(brewer.pal(10, "Spectral"))
plot(scot, col = scale10)
axis(1)
axis(2)

spplot(scot, "Cases")

Convert coordinates

# what coordinate reference are we?
proj4string(scot)

## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"

# make lat-long
scotLL = spTransform(scot, CRS("+init=epsg:4326"))
extent(scotLL)

## class       : Extent 
## xmin        : -8.621 
## xmax        : -0.751 
## ymin        : 54.63 
## ymax        : 60.85

# make Google Mercator
scotM = spTransform(scot, CRS("+init=epsg:3857"))

Lets have a background

require(dismo)
bg = gmap(extent(scotM))
plot(bg, inter = TRUE)
plot(scotM, add = TRUE)

OpenStreetMap background

require(OpenStreetMap)
e = extent(scotLL)
scotmap = openmap(c(ymax(e), xmin(e)), c(ymin(e), xmax(e)))
plot(scotmap)
plot(scotM, add = TRUE)

Spatial Points and Centroids

scotC = SpatialPoints(scotM)
plot(scotmap)
plot(scotC, add = TRUE, cex = scale10/3, col = "#FF000080", pch = 19)

Alternative

scotT = openmap(c(ymax(e), xmin(e)), c(ymin(e), xmax(e)), type = "stamen-toner")
plot(scotT)
palette(brewer.pal(10, "Spectral"))
plot(scotC, add = TRUE, col = scale10, pch = 19, cex = 3)

Classifying and legends

require(classInt)
# classify rate into groups
classify = classIntervals(scot$Cases/scot$Expected, style = "pretty")
mypal <- rev(brewer.pal(10, "Spectral"))
# assign colours to groups
rateColours = findColours(classify, mypal)
# plot with colours
plot(scot, col = rateColours)
legend("right", fill = attr(rateColours, "palette"), legend = names(attr(rateColours, 
    "table")), title = "Rate")

Doesn't matter if code looks ugly - you only need to type it once

Transparent Colours

Colour #FF0000 is red. #FF000040 is a transparent red...

mypal <- rev(brewer.pal(10, "Spectral"))
rateColours = findColours(classify, mypal)
plot(scotmap)
plot(scotM, add = TRUE, col = paste0(rateColours, "40"))

Lines

birds = readOGR("./Data", "birds")

## OGR data source with driver: ESRI Shapefile 
## Source: "./Data", layer: "birds"
## with 2 features and 3 fields
## Feature type: wkbLineString with 2 dimensions

summary(birds)  # what have we?

## Object of class SpatialLinesDataFrame
## Coordinates:
##      min   max
## x -17.48  1.40
## y  14.68 56.64
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
##       id         age         gender 
##  Arthur:1   Min.   :3.0   Female:1  
##  Bravo :1   1st Qu.:3.5   Male  :1  
##             Median :4.0             
##             Mean   :4.0             
##             3rd Qu.:4.5             
##             Max.   :5.0

par(mfrow = c(1, 2))
plot(birds, col = c("red", "blue"))
axis(1)
axis(2)
birds = spTransform(birds, CRS("+init=epsg:3857"))
plot(scotmap)
plot(birds, add = TRUE, col = c("red", "blue"))

par(mfrow = c(1, 1))

Conversion to points

Suppose I want to map the pre-migration locations of bird 1...

bird1 = as(birds[1, ], "SpatialPoints")
bird1 = bird1[1:93, ]
e = extent(spTransform(bird1, CRS("+init=epsg:4326")))
map = openmap(c(ymax(e), xmin(e)), c(ymin(e), xmax(e)))
plot(map)
plot(bird1, add = TRUE, col = "red")

Density Estimation

require(scales)  # for nice colours
suppressPackageStartupMessages(require(sparr))
e = extent(bird1)
birdsmooth = bivariate.density(coordinates(bird1), xrange = c(xmin(e), xmax(e)), 
    yrange = c(ymin(e), ymax(e)), pilotH = 100, res = 100, comment = FALSE)

## Warning: WIN not specified - using rectangular region defined by xrange
## and yrange arguments

plot(map)
plot(birdsmooth, add = TRUE, col = alpha("red", c(0, seq(0.1, 1, len = 99))))
points(bird1, pch = 19, cex = 0.2)
lines(birds[1, ], col = "#20202080", lty = 2)

Perspective

persp(birdsmooth$X, birdsmooth$Y, birdsmooth$Zm, theta = 23)

require(hexbin)
plot(hexbin(coordinates(bird1)))

Export to Google Earth

writeOGR(birds, "./Data/birds.kml", "osprey", "KML")

Quantum GIS

Rasters

Creation

The raster package and friends

require(raster)
# various construction methods
# 
# from a matrix
m = matrix(1:1200, 30, 40)
r = raster(m)
r = raster(m, xmn = -10, xmx = 10, ymn = 54, ymx = 60)

# from numbers
r = raster(ncols = 30, nrows = 40, xmn = -10, xmx = 10, ymn = 54, ymx = 60)
r[] = runif(1200)

# from GDAL data files
landuse = raster("./Data/cumbriaLandUse.tif")
dem = raster("./Data/cumbriaFullDem.tif")
plot(landuse)

plot(dem)

Raster Properties

dim(dem)

## [1] 956 827   1

extent(dem)

## class       : Extent 
## xmin        : 3442100 
## xmax        : 3524800 
## ymin        : 3514200 
## ymax        : 3609800

projection(dem)

## [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"

summary(dem)

## Warning: summary is an estimate based on a sample of 1e+05 cells (12.65%
## of all cells)

##         cumbriaFullDem
## Min.              0.00
## 1st Qu.          17.67
## Median          133.36
## 3rd Qu.         270.70
## Max.            925.40
## NA's              0.00

res(dem)  # cell size

## [1] 100 100

Raster Ops

thinDem = aggregate(dem, 5)

plot(thinDem * (thinDem > 300))

Values

dem[500, 400]

##       
## 91.13

hist(values(dem))

Point Sampling

cellFromXY(dem, c(3460000, 3540000))

## [1] 577426

dem[cellFromXY(dem, c(3460000, 3540000))]

##       
## 43.53

x = seq(3464670, 3476071, len = 100)
y = seq(3555793, 3567975, len = 100)
plot(dem)
lines(x, y, lwd = 3)

trans = dem[cellFromXY(dem, cbind(x, y))]
plot(trans, type = "l")

Projection

demLL = projectRaster(thinDem, crs = "+init=epsg:4326")
plot(demLL)

Terrain

demT = terrain(thinDem, opt = c("slope", "aspect"))
demShade = hillShade(demT[[1]], demT[[2]])
plot(demShade, col = grey(0:100/100), legend = FALSE)