R

Spatial Cheatsheet

calsalmon April 1, 2015

This cheatsheet is an attempt to supply you with the key functions and manipulations of spatial vector and raster data. It does not have examples for you to cut and paste, its intention is to provoke the “Oh yes, that’s how you do it” thought when stuck.

Packages

library(sp)  # vector data
library(raster)  # raster data
library(rgdal)  # input/output, projections
library(rgeos)  # geometry ops
library(spdep)  # spatial dependence

Links

Points

# points from scratch
coords = cbind(x, y)
sp = SpatialPoints(coords)
# make spatial data frame
spdf = SpatialPointsDataFrame(coords, data)
spdf = SpatialPointsDataFrame(sp, data)
# promote data frame to spatial
coordinates(data) = cbind(x, y)
coordinates(data) = ~lon + lat
# back to data
as.data.frame(data)

##      lon   lat Z
## 1 11.515 24.52 d
## 2  7.056 27.11 a
## 3 12.945 30.09 c
## 4 12.793 24.72 e
## 5 12.888 28.24 b

data@data

##   Z
## 1 d
## 2 a
## 3 c
## 4 e
## 5 b

bbox(spdf)

##      min   max
## x  7.056 12.94
## y 24.520 30.09

Lines

c1 = cbind(x1, y1)
c2 = cbind(x2, y2)
c3 = cbind(x3, y3)

# simple line strings
L1 = Line(c1)
L2 = Line(c2)
L3 = Line(c3)

# single/multiple line strings
Ls1 = Lines(list(L1), ID = "a")
Ls2 = Lines(list(L2, L3), ID = "b")

# with spatial nature
SL1 = SpatialLines(list(Ls1))
SL12 = SpatialLines(list(Ls1, Ls2))

# made into spatial data frame
SLDF = SpatialLinesDataFrame(SL12, data.frame(Z = c("Road", "River"), row.names = c("a",
    "b")))

as.data.frame(SLDF)

##       Z
## a  Road
## b River

SpatialLinesLengths(SLDF)

## [1] 2.414 4.828

Polygons

# single ring feature
c1 = cbind(x1, y1)
r1 = rbind(c1, c1[1, ])  # join
P1 = Polygon(r1)
Ps1 = Polygons(list(P1), ID = "a")

# double ring feature
c2a = cbind(x2a, y2a)
r2a = rbind(c2a, c2a[1, ])
c2b = cbind(x2b, y2b)
r2b = rbind(c2b, c2b[1, ])

P2a = Polygon(r2a)
P2b = Polygon(r2b)
Ps2 = Polygons(list(P2a, P2b), ID = "b")

# Spatial Polygons Data Frame
SPs = SpatialPolygons(list(Ps1, Ps2))

SPDF = SpatialPolygonsDataFrame(SPs, data.frame(N = c("one", "two"), row.names = c("a", "b")))
SPDF@data

##     N
## a one
## b two

# single ring with hole
hc1 = cbind(xh1, yh1)
hr1 = rbind(hc1, hc1[1, ])
H1 = Polygon(hr1, hole = TRUE)
P1h = Polygons(list(P1, H1), ID = "c")
SP1h = SpatialPolygons(list(P1h))
# plot(SP1h,usePolypath=TRUE)

Raster

# from x,y,z-matrix
r1 = raster(list(x = x, y = y, z = z))
# rows and columns values
r1[, 3]

##  [1] 0.7377 0.3342 0.6924 0.3482 0.2972 0.8148 0.8212 0.5362 0.8750 0.9808 0.2729

r1[2, ]

##  [1] 0.40396 0.79350 0.33422 0.25095 0.64577 0.88173 0.50432 0.73244 0.98500 0.13277 0.59993 0.04035

extent(r1)

## class       : Extent 
## xmin        : 22.95 
## xmax        : 29.45 
## ymin        : 44 
## ymax        : 46.3 

dim(r1)

## [1] 11 12  1

# create empty, then fill
r2 = raster(nrows = nrows, ncols = ncols, xmn = xmn, xmx = xmx, ymn = ymn, ymx = ymx)
r2[] = runif(nrows * ncols)

# create from extent, then set values
r3 = raster(extent(r2), nrows = nrows, ncols = ncols)
values(r3) = runif(nrows * ncols)

# multi-band stack
s1 = stack(list(r1, r2, r3))
dim(s1)

## [1] 11 12  3

# multi-band brick
b1 = brick(list(r1, r2, r3))

Coordinates

# EPSG strings
latlong = "+init=epsg:4326"
ukgrid = "+init=epsg:27700"
google = "+init=epsg:3857"

# Spatial*
proj4string(SPDF)

## [1] NA

proj4string(SPDF) = CRS(latlong)
SL1 = SpatialLines(list(Ls1), proj4string = CRS(latlong))

# Raster CRS
projection(r1)

## [1] "NA"

# - assign or set on creation
projection(r1) = CRS(latlong)
r1 = raster(list(x = x, y = y, z = z), crs = latlong)

# Transform Spatial*
SPtrans = spTransform(SPDF, CRS(google))

# Transform/Warp Raster
rTrans = projectRaster(r1, crs = google)

I/O

# -- vectors

# avoid - doesn't read CRS
library(maptools)
shapes = readShapeSpatial("data.shp")

# read/write shapefiles (and others)
#  - list formats
ogrDrivers()
shapes = readOGR(".", "data")
writeOGR(shapes, ".", "data", "ESRI Shapefile")
writeOGR(shapes, "towns.kml", "towns", "KML")

# -- rasters

# creates SpatialGrid objects
r = readGDAL("data.tif")

# create Rasters/Brick objects from files
r = raster("data.tif")
# - write Raster to GeoTIFF
writeRaster(r, "data2.tif", "GTiff")
# - supported formats
writeFormats()
# or for Google Earth
KML(r, "r.kmz")

Manipulation

# Spatial*DataFrames
#
# subset(Towns,pop>29000) doesn't work
BigTowns = subset(Towns, Towns$pop > 29000)
BigTowns = Towns[Towns$pop > 29000, ]
BigAndSmall = rbind(BigTowns, SmallTowns)

# points in zones
Towns[1:10, ] %over% Zones

##  [1]  1 NA  2 NA NA  2 NA NA NA  1

# rasters
# - sample points
vpt = extract(r1, sPoints)
# - sample polygons
vpoly = extract(r1, sPolys)
# - crop an area
e = extent(raster(xmn = 25, xmx = 27, ymn = 44.5, ymx = 45.5))
cr1 = crop(r1, e)

Plotting

# scale colour
library(RColorBrewer)
palette(brewer.pal(6, "YlOrRd"))
plot(Towns, col = plotrix:::rescale(Towns$pop, c(1, 6)), pch = 19)

# scale size
plot(Towns, cex = plotrix:::rescale(Towns$pop, c(1, 4)), pch = 19)

# polygons
plot(Zones, col = fillColour, border = outlineColour)

# sp colours
sp.theme(set = TRUE, regions = list(col = colours))
spplot(Towns, "pop")

# rasters
plot(r1)

# true colour from bricks
plotRGB(b1, scale = 1)