generalized code for plotting rasters/interpolated maps

RasterPlotting.R

# GIS libraries 
library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(mapmisc)

## defines the coordinate reference system
utm17=CRS("+init=epsg:26917")

## Example Data
## No data loading, just explainations
lakeO # this is a SpatialPolygonDataFrame outline of a Lake Okeechobee
mud.pts # sampling location

# These are interpolated layers clipped to lakeO
mud1 
mud2
mud3
mud4

# I save the maps/figures using png(..) but you can also use tiff(),pdf(),etc
# png(filename="example.png",width=6.5,height=5,units="in",res=200,type="windows",bg="white")
# set up plotting window spacing and different frames using par() and layout()
par(family="serif",mar=c(0.1,0.1,0.1,0.1),oma=c(0.5,1,1,0.5));
layout(matrix(1:5,1,5),widths=c(1,1,1,1,0.5))

# estblish your breaks points in the data
# this represents the data from 0 to 140 in 10 unit increments
b=seq(0,140,10)
# select your palatte for each break (it needs to be 1 less than breaks)
pal=viridis::plasma(length(b)-1,alpha = 0.8)

# set up bounding box of the map
bbox.lims=bbox(lakeO)
## 1st plot
# alittle hacky but I plot the polygon then add to it.
plot(lakeO,ylim=bbox.lims[c(2,4)],xlim=bbox.lims[c(1,3)],lwd=0.05)
# This is the interpolation image() seems to be the most flexible, makes sure to had add=TRUE
image(mud1,add=T,breaks=b,col=pal)
# this draw contours of the interpolation based on you breaks
# you can also use contour(...) but its less customizable
plot(rasterToContour(mud1,levels=b,nlevels=length(b)),col=adjustcolor("white",0.5),lwd=0.5,add=T)
# adds sampling points
plot(mud.pts,pch=19,add=T,col="grey",cex=0.25)
# addes scale bar and north arrow...there are other packages and tools for this
mapmisc::scaleBar(utm17,"bottom",bty="n",cex=1,seg.len=4,outer=T)
# add text to top or anywhere you want
mtext(side=3,line=-1,1988)

## 2st plot
# repeat for mud2 
## 3st plot
# repeat for mud3
## 4st plot
# repeat for mud4

# the custom legend
# make an empty plot
plot(0:1,0:1,ann=F,axes=F,type="n")
# make the color ramp graphic
legend_image=as.raster(matrix(pal,ncol=1))
# where you want things 
top.val=0.65
bot.val=0.35
x.max=0.50
x.min=0.25
# text on the side
text(x=x.max, y = c(top.val,bot.val), labels = c(0,max(b)),cex=0.75,adj=0,pos=4,offset=0.5)
# draw the color ramp
rasterImage(legend_image,x.min,bot.val,x.max,top.val)
# legend title
text(x=x.min+((x.max-x.min)/2),y=top.val,"Mud Depth (cm)",adj=0,pos=3,cex=0.8,xpd=NA)

# can also do legend(...)
dev.off() # closing the png/pdf/jpeg/etc connection and finalizes the map

Thanks a lot. Very cool.