2

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")
h.l.m
  • 13,015
  • 22
  • 82
  • 169

1 Answers1

3

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.

Edzer Pebesma
  • 3,814
  • 16
  • 26
  • I don't see a .grd type file...i only see .asc .gml and .prj – h.l.m Feb 28 '15 at 23:53
  • Also how do i get the data.frame with the three columns of lat, long, altitude? – h.l.m Mar 01 '15 at 00:03
  • 1
    you are right: extension is `.asc`, not `.grd`, I corrected the answer above; you can get the `data.frame` by `as.data.frame(r100, xy = TRUE)` – Edzer Pebesma Mar 01 '15 at 14:12
  • Thanks for that, that is very useful. So I noticed that when carrying out `as.data.frame(r100, xy = TRUE)` the x and y values are not latitude and longitude, how would I convert it to be so... – h.l.m Mar 01 '15 at 15:29
  • 1
    OK, I added this to the answer. – Edzer Pebesma Mar 01 '15 at 20:15