Creating DEMs from contour lines using R

Download 3m contour line data files here.
Download sub-basin data files here.
Download R script here.

library(sp)
library(maptools)
library(rgdal)
library(raster)
library(rgeos)

#Read in shape files
setwd("~/Documents/chitchatr/data/Contours_3m_Lines")
contour3m <- readShapeLines("GISDATA_CONTOURS5K_ARC",
                            proj4string = CRS("+init=epsg:26918"))

setwd("~/Documents/chitchatr/data/Subbasins_Outlines")
ws_bound <- readShapePoly("GISDATA_SUBBASINS_POLY",
                          proj4string = CRS("+init=epsg:26918"))

# plot of ws_bound with watershed of interest highlighted
spplot(ws_bound, zcol = "SUB_ID", col = "black", colorkey = F,
       fill = "grey80", 
       panel = function(x, y, ...){
         panel.polygonsplot(x,y, ...)
         sp.polygons(ws_bound[as.character(ws_bound$SUB_ID) == 1007, ],
                     fill = "magenta")
       })

# same plot as above, but include contour3m lines
spplot(ws_bound, zcol = "SUB_ID", col = "black", colorkey = F,
       fill = "grey80", 
       panel = function(x, y, ...){
         panel.polygonsplot(x,y, ...)
         sp.polygons(ws_bound[as.character(ws_bound$SUB_ID) == 1007, ],
                     fill = "magenta")
         sp.lines(contour3m, col = "black")
       })

# Create raster grid where we are going to interpolate our
# elevation data.
contourBBox <- bbox(contour3m)
contourRaster <- raster(xmn=contourBBox[1,1], 
                        xmx=ceiling(contourBBox[1,2]),
                        ymn=contourBBox[2,1], 
                        ymx=ceiling(contourBBox[2,2]),
                        crs = CRS("+init=epsg:26918"),
                        resolution = c(12,12))

# Create points from contour lines so we can
# interpolate between points.
p <- as(contour3m, 'SpatialPointsDataFrame')

#interpolate between points using IDW into contourRaster
library(gstat)
g <- gstat(id="ELEV_M", formula = ELEV_M~1, data=p,
           nmax=7, set=list(idp = .5))
interpDEM <- interpolate(contourRaster, g)

# extract the watershed boundary of interest
HudsonWS_bound <- ws_bound[as.character(ws_bound$SUB_ID) == 1007, ]
# clip raster to watershed boundaries
interpDEM_mask <- mask(interpDEM, HudsonWS_bound)
# clip contours to watershed boundary
HudsonWS_contour <- crop(contour3m,HudsonWS_bound)

# subset 3m contours for nice plotting
HudsonWS_contourSub <- HudsonWS_contour[(HudsonWS_contour$ELEV_M) %in% 
                                          seq(354, 993, 30), ]
# Plot it up!
spplot(interpDEM_mask, col.regions = heat.colors(100),
       xlim=c(bbox(interpDEM_mask)[1,1] + 900,
              bbox(interpDEM_mask)[1,2] - 900),
       ylim=c(bbox(interpDEM_mask)[2,1] + 1150,
              bbox(interpDEM_mask)[2,2] - 900),
       colorkey = list(space = "top", width = 1), main = "Elev. (m)",
       panel = function(x, y, ...){
         panel.gridplot(x,y, ...)
         sp.lines(HudsonWS_contourSub, col = "black")
       })
Advertisements

One thought on “Creating DEMs from contour lines using R

  1. Pingback: Plotting DEMs in 3D using persp(): Part 2 | Chit Chat R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s