This section is an ongoing repository for spatial functions and code snippets that I’ve found useful, for my benefit as much as anyone else’s. Feel free to use without attribution.

Operation {library} Code
open a shapefile {maptools} pts = readShapePoints(‘shapefile.shp’)
lns = readShapeLines(‘shapefile.shp’)
plys = readShapePoly(‘shapefile.shp’)
extract polygon centroids as a spatialPoints object {rgeos} gCentroid(spatPolygonsObj, byid=T)
extract point coordinates/polygon centroids
{sp / maptools}
coordinates(spatial_obj)
data.frame from spatial d.f. spatial_df@data
query Coordinate Reference System (CRS) spatial_obj@proj4string
assign a CRS proj4 string{maptools} (where absent), for:
1. latlongs (WGS84)
2. standard Mercator

3. UK (British National Grid / OSGB36)

4. USA (I think)

5. UTM (e.g. Zone 16 Belize)

spatial_df@proj4string = . . .

1. CRS(“+proj=longlat +datum=WGS84″)
2. CRS(“+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs”)
3. CRS(“+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs”)
4. CRS(“+proj=utm +zone=19 +datum=NAD27 +ellps=clrk66 +towgs84=-3,142,183,0,0,0,0″)
5. CRS(“+proj=utm +zone=16 ellps=WGS84″)

useful table of proj4 strings for spatial projections / transformations {rgdal} epsg = make_EPSG()
View(epsg)
re-project a spatial object {rgdal} spTransform(spatial_df, CRS(“new proj4 string”))
function to clip a spatial object to the bounding box of another (e.g. ‘point in polygon’) {rgeos}

[N.B. you may need to change the quotation marks]

clipshape = function(shapetoclip, baselayer){
e = extent(baselayer)
bbx = readWKT(paste(“POLYGON((“,e@xmin,e@ymin,”,”,e@xmax,e@ymin, +
“,”,e@xmax,e@ymax,”,”,e@xmin,e@ymax,”,”,e@xmin,e@ymin,”))”))
gIntersection(shapetoclip, bbx)
}
aggregate a raster object {raster} (e.g. by a factor of 10) aggregate(raster_obj, 10)
merge adjacent raster tiles {raster} region = merge(N17W089, N17W090, N18W089, N18W090)
clip a raster tile {raster} crop(region, extent([xmin], [xmax], [ymin], [ymax]))
re-project a raster object projectRaster(region, crs=“new proj4 string”)
calculate great circle between 2 points
{geosphere + sp}
– set CRS e.g. to “+proj=longlat +datum=WGS84” for latlons
– (opt) n vertices / string ID
greatcircle = function(p1, p2, n=50, ID=NA, crs=NA){
gcLine = Line(gcIntermediate(p1, p2, n=n))
gcLines = Lines(gcLine, ID)
SpatialLines(list(gcLines), proj4string = CRS(crs))
}
Lon = c(0,52); NY = c(-75,43)
plot(greatcircle(Lon, NY, n=100), axes=T)
great circles between point lists
{maps + geosphere}
– also accepts spatial objects
data(world.cities)
(fm = world.cities[sample(1:nrow(world.cities), 15),c(1,5,4)])
(to = world.cities[sample(1:nrow(world.cities), 15),c(1,5,4)])
gcircles = gcIntermediate(fm[,2:3], to[,2:3], breakAtDateLine=T)
plot.new(); plot.window(xlim=c(-180,180), ylim=c(-90,90))
for(path in gcircles) if(length(path) != 2) lines(path) # not crossed meridian
text(rbind(fm[,2:3], to[,2:3]), c(fm$name, to$name), cex=0.6)
vector of raster values {raster}
– i.e. elevation profile of a path
extract(raster_obj, spatialLines_obj)
locating/grouping raster pixels by value {raster} r = raster(volcano, crs=CRS(‘+proj=longlat +datum=WGS84′))
r@extent = extent(-76,-74,42,44)
p = data.frame(rasterToPoints(r))
peak = p[p[,3] == r@data@max,1:2] # max value
low = p[p[,3] == r@data@min,1:2] # min value
plot(r)
points(peak, pch=15, col=’red’)
points(low, pch=15, col=’blue’)
Get profile of raster values along a (poly)line
{raster, sp}
[using ‘r’ from above]
line_obj = Lines(Line(matrix(c(-75,-74.02,43.55,42.13),2)), ID=NA)
path = SpatialLines(list(line_obj), proj4string = r@crs)
profl = extract(x=r, y=path)[[1]];
plot(profl, type=’l’, col=’red’, xlab=”, ylab=’meters’, xaxt=’n’)