Using R I'd like to create a raster from a set of contour shape files, where each cell in the raster contains a value for the derived gradient.

I am using OS Terrain 50 contour lines which is open data from, the Ordnance Survey: https://www.ordnancesurvey.co.uk/opendatadownload/products.html


As requested, here is how I've approached the problem so far, which hasn't really worked and I'm not sure why.

Firstly I load my contour shapefile:

c <- readOGR("osrastersandvectors/terrain_shape_files/HP40_line.shp")

Then the following to create a raster at 10m resolution:

c.raster <- raster()
extent(c.raster) <- extent(c)
res(c.raster) <- 10

Then I rasterize:

c.raster <- rasterize(c,c.raster)

This works in that I get a raster of the original contour lines, where cell values are elevation. Note, there are no elevation values beyond where the original contour lines are. Then I'm thinking I need to calculate a value for each cell, given the contour line values. So some form of interpolation to give slope (gradient?). So something along the lines of:

slope = terrain(c.raster, opt=c('slope'), unit = 'degrees', neighbours=8)

This returns a raster but it's not what I want: a single layer raster where each cells value represents the gradient at that point.

  • Welcome to GIS SE! For questions that involve code we ask that you show us where you are stuck with your own code by including a code snippet in your question. Please check out our short tour for more about how this site works.
    – PolyGeo
    Aug 1, 2019 at 7:11
  • This can be done as a two-step process - first convert the vector contour data to an elevation raster, and then compute the gradient raster from the elevation raster. Pretty sure those processes are doable and are documented and are possibly also answers to questions on here already...
    – Spacedman
    Aug 1, 2019 at 9:11

1 Answer 1


So I found a great blog post here that answers my question.

# Load the contour data and subset by layer required
readOGR("myfolder/contours.shp") %>%
  subset(FEAT_TYPE =="ContourLine") -> c

# Create a target raster for the DEM    
dem_bbox <- extent(c)
dem_raster <-raster(dem_bbox)
projection(dem_raster) <- CRS(projection(c))
# Set resolution at 10m
res(dem_raster) <- 10

# Convert contour to df of spatial points
dem_points <- as(c, "SpatialPointsDataFrame")

# Compute the interpolation function
dem_interp <- gstat(formula = elevation ~ 1, locations = dem_points,
                    set = list(idp = 0), nmax = 5)

# Obtain interpolation values for raster grid
DEM <- interpolate(dem_raster, dem_interp)

# Derive derive additional terrain vars from DEM using terrain fn.
derived_vars <- terrain(DEM, opt = c('slope', 'roughness', 'aspect'))
slope <- derived_vars[["slope"]]
roughness <- derived_vars[["roughness"]]
aspect <- derived_vars[["aspect"]]

This all works fine. I do have an additional question. Is there an equivalent way to implement this using the sf library?

  • Did you find out if there is a method to apply this to sf objects? Would one need to convert it to a raster first following here: stackoverflow.com/questions/64806684/…?
    – QPaps
    Feb 8, 2021 at 15:39
  • A further question - in your gstat() formula - what is elevation?
    – QPaps
    Feb 8, 2021 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.