Search code examples
rgisrgdalr-maptoolsr-sp

R Converting contour lines to elevation plot

I would like to be able to create an elevation plot from contour lines in R. I am very new to using shape files

At the moment I have downloaded data from here which provides .shp files for all of the UK.

It also provides the contour lines, summarising the topology of the UK.

For the elevation plot I would like a data.frame or data.table of evenly spaced points (100m apart from each other) to produce a data output giving an x, y and z value. Where x and y represent the latitude and longitude (or Eastings and Northings), and z represent the height (in meters above sea-level).

I think there are probably some tools that will automatically carry out the interpolation for you, but am unsure how it would work with geo-spatial data.

This is my basic start...

require(maptools)
xx <- readShapeSpatial("HP40_line.shp")
Solution

  • Choose "ASCII Grid and GML (Grid)" as download format for the "OS Terrain 50" product, and download the file. This will give you a zip file containing many directories of zip files, each of which contains portions of a 50 m elevation grid of the UK (the portion I looked at had 200 x 200 cells, meaning 10 km x 10 km). I went into the directory data/su, unzipped the zip file there, and did

    library(raster)
r = raster("SU99.asc")
plot(r)

    to aggregate this to a 100 m grid, I did

    r100 = aggregate(r) # default is factor 2: 50 -> 100 m

    As mentioned above, the advice is to work on the grids as contour lines are derived from grids, working the other way around is a painful and a great loss of information.

    Getting grid values in longitude latitude as a data.frame can be done in two ways:

    df = as.data.frame(projectRaster(r, crs = CRS("+proj=longlat")), xy = TRUE)

    unprojects the grid to a new grid in longitude / latitude. As these grids cannot coincide, it minimally moves points (see ?projectRaster).

    The second option is to convert the grid to points, and unproject these to longitude latitude, by

    df2 = as.data.frame(spTransform(as(r, "SpatialPointsDataFrame"), CRS("+proj=longlat")))

    This does not move points, and as a consequence does not result in a grid.