Simulate sea level change using R

Question

Dropbox linked files I am trying to simulate sea level change using R but have no idea how to do that. Similar desired result. Similar desired result I have seen an article which says the person used contour() but dont know how to use this function. I have also seen a lot of people using ggplot2 but am not sure how to implement a dataframe using my data set for that to work.

I have a dtm which is a combined version of 5 asc tiles. using the code below.

setwd()
f <-list.files(pattern = ".asc")  
r <- lapply(f, raster) 
x <- do.call("merge",r) 
writeRaster(x,"DTM_combine.asc", overwrite=TRUE)  
library(rgdal)
r = raster("DTM_combine.asc")
plot(r)

I then reclassified the raster using this.

image(r,zlim=c(0,70), main="DEM Findhornbay", col=col) m=c(0,5,1,5,10,2,10,15,3,15,20,4,20,25,5,25,30,6,30,35,7,35,40,8,40,45,9,45,50,10,50,55,11,55,60,12,60,65,13,65,70,14)
mat=matrix(m,ncol=3,byrow=TRUE)
r=reclassify(r,mat)
rcat=reclassify(r,mat)
plot(rcat)
col <-terrain.colors(14)
plot(rcat)
brk <-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
plot(rcat, col=col, breaks=brk, main="Recalssed DEM Findhorn Bay")
plot(r, col=col, breaks=brk, main="Recalssed DEM Findhorn Bay")

which produced this. I am not sure if this is the correct way for this to work, but the next objective is to simulate sea level change.

plot3D(r)

Answer 1

Raster package and base plot

To draw contours using rasters and the capabilities provided by the raster package you proceed as follows:

First load packages, declare your crs, load your raster and assign your crs:

library(raster)
library(magrittr)

#We define the crs of the project
crs <- "+proj=utm +datum=WGS84 +no_defs"

#read your raster DTM_combine.asc and assign it the projection
r <- raster::raster("PathToFolder/DTM_combine.asc")
raster::projection(r) <- crs

Insofar as we are working with a digital elevation model, I'm assuming that elevation units are in meters. We can define a vector of elevations niveaux as a sequence between 0 and max elevation value every 5 meters:

niveaux <- (raster::extract(r, raster::extent(r)) %>% 
              range(finite = TRUE))[2] %>% 
                seq(0,.,by=5)

Then you can plot the raster:

raster::plot(r)

Then add your contours. Here I'm adding the coastal line (elevation =0) in blue and a retreat of such coastal line (in red) in the assumption that the sea level rises 5 meters.

raster::contour(r, add=TRUE, col="blue", lwd=0.2, levels=niveaux[1], drawlabels=FALSE)
raster::contour(r, add=TRUE, col="red", lwd=0.2, levels=niveaux[2], drawlabels=FALSE)

Alternatively, you could work only with the contours if you plot them as follows:

raster::contour(r, col="blue", lwd=0.2, levels=niveaux[1], drawlabels=FALSE)
raster::contour(r, add=TRUE, col="red", lwd=0.2, levels=niveaux[2], drawlabels=FALSE)

The ggplot2 way

First, convert your raster in a dataframe:

r.df <- raster::as.data.frame(r, xy = TRUE)

Then plot the raster with geom_raster and make the contours with geom_contour in the same way we did before:

  ggplot(data = r.df, aes(x = x, y = y)) +
  geom_raster(aes(fill=DTM_combine)) +
  geom_contour(aes(z=DTM_combine, color = "Current"), breaks = niveaux[1], size=0.5) + 
  geom_contour(aes(z=DTM_combine, color = "New"), breaks = niveaux[2], size=0.5) +
  scale_fill_gradientn(colours = terrain.colors(10, rev=TRUE), na.value = "transparent") +
  theme_void() +
  labs(fill="Elevation", color="Coastal lines")

The resulting map:

If you want to avoid the white square in the upper left corner you can fix your na.value in scale_fill_gradientn() to #f2f2f2 for example.

*Bonus: If you want to fill the space between coastal lines you can have a look at the function geom_contour_fill() in the metR package.

Hope it helps