Creating DEMs from contour lines using R

UPDATE: I made a slight change to the script, because I was getting an error when I ran the script on Windows, but not on Linux.  Previously, the script read:

# 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))

But I found, that the resolution = c(12,12) argument was giving an error on my Windows computer.  Fortunately, there is an easy fix. Just remove the argument and add the line res(contourRaster) <- 12:

# 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"))
res(contourRaster) <- 12

Everything should work smoothly now!

Quick Note: A few years ago I did a post on using R for spatial analyses. I recently thought I’d return to the topic, mainly to show how versatile R is and how many of the things you can do in dedicated GIS software (which can cost thousands of dollars to license!) can be re-created easily and quickly using R.

DEMs, or digitial elevation models, are one of the most common types of data used in spatial analysis. DEMs are raster files that contains elevation information for each cell in the raster. I recently was give a shapefile of contour lines for a particular area and needed to create a DEM from them. As such, I thought I’d demonstrate how I created a DEM from the contour lines that were in a vector shapefile (.shp), one of the most common GIS file types.

I won’t go into all the differences between rasters and vectors (ArcGIS has some good online resources in case you need to brush up) but the basic steps are:

1. Read in the files to R.
2. Create a raster “grid” to interpolate the elevation data onto.
3. Convert the contour lines to points so you can interpolate between elevation points.
4. Interpolate the elevation data onto the raster grid.
5. Crop the DEM raster to the boundary of the area of interest.

Before I begin, I obtained my data from the Mass GIS webportal, which is a really nice web-based GIS for querying all types of spatial data within the Commonwealth of Massachusetts. Since I’m interested in elevation data only, I looked for a sub-basin (ie. watershed) that had a nice topographic range and shape. I then downloaded the watershed boundaries, and the 3m contour lines that covered the watershed, as separate shapefiles.

Step 1. This is quite simple, since there are some really nice R packages for dealing with spatial data. To read in our data, we will use the maptools library. You’ll notice that the readShapeLines() and readShapePoly() has a proj4string argument. This is where you tell R what the projection and datum are for your data. MassGIS actually asks which projection you’d like your data in, so I knew ahead of time what EPSG I needed to input:

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

#Read in shape files
setwd("C:/LocationOfYourContourData")
contour3m <- readShapeLines("GISDATA_CONTOURS5K_ARC",
                            proj4string = CRS("+init=epsg:26918"))

setwd("C:/LocationOfYourBoundaryData")
ws_bound <- readShapePoly("GISDATA_SUBBASINS_POLY",
                           proj4string = CRS("+init=epsg:26918"))

Lets plot up the files to see what we are working with. I’m interested in the watershed in the center, but you can see that the downloaded data included adjacent watersheds. Also, the contours extend beyond the watershed of interest. We’ll deal with this later when we crop and clean up our data:

# 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(HudsonWS_contourSub, col = "black")
       })
watershed boundaries downloaded from MassGIS.  The watershed of interest is highlighted in magenta.

Watershed boundaries downloaded from MassGIS. The watershed of interest is highlighted in magenta.

Added 3m contour lines to plot to show extent of contour lines outside watershed of interest.

Added 3m contour lines to plot to show extent of contour lines outside watershed of interest.

Step 2. I need to create an empty raster grid that spans the range of my sub-basin. Since I’m using the contour lines to create the DEM, I used them as the basis for my raster grid. Notice that I manually set the resolution (size) of each cell in the raster to be square (12m x 12m). This is required if you ever want to import your data into ArcGIS (as far as I’m aware). I also set the resolution so that I didn’t have too many cells. The resolution should depend on your needs and computing power:

# 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))

Step 3. Since interpolation methods were conceived to work with point data, we need to convert the elevation contour lines to elevation points. Fortunately, this is quite easy to do in R using as(). Essentially we are creating points along each contour line that has as its value the elevation of the associated contour line:

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

Step 4. Now we have points from which we can interpolate a elevation raster. There are probably many different intrpolation packages out there, but I’m going to use gstat, since it seems pretty robust for spatial statistics, and there are nice examples using the gstat and raster packages together for interpolating point data to rasters. We will use gstat() to create a model for interpolating our point elevation data using inverse distance weighting. Then, we use interpolate() to predict the elevation value at each cell in contourRaster. Note that this step can take some time depending on the computing power of your computer:

#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)

Step 5. Finally, we can clean up the data by clipping our DEM to the boundaries of the watershed. We can also clip our contour lines to the boundaries of our watershed:

# 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)

That wasn’t too bad. We can plot it up nicely using spplot() like I did in a previous post:

# 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")
       })
Final DEM with 30m contour lines.

Final DEM with 30m contour lines.

Not too shabby! Anyway, keep an eye on Chit Chat R for more spatial analysis in near future…

Advertisements

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