Skip to content

Instantly share code, notes, and snippets.

@SwampThingPaul
Last active April 21, 2022 19:04
Show Gist options
  • Save SwampThingPaul/ef0229d76d1253935f9ef1801dc51127 to your computer and use it in GitHub Desktop.
Save SwampThingPaul/ef0229d76d1253935f9ef1801dc51127 to your computer and use it in GitHub Desktop.
generalized code for plotting rasters/interpolated maps
# 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
@Ignimbrit
Copy link

Thanks a lot. Very cool.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment