Introduction

In the last class we looked at the tables held in the shared online postgis database. After a discussion about the nature of spatial data and coordinate reference systems we selected the SSSI polygons for Arne RSPB reserve and loaded them into memory. Let’s pick up from there and get some more data.

## R packages needed
library(rpostgis)
library(RPostgreSQL)
library(sp)
library(raster)
library(rgeos)
load("spatial_class1.rob")
plot(arne,axes=TRUE)

OK. So we have the polygon in memory again. We first need to reconnect to postgis in order to get some more data.

conn <- dbConnect("PostgreSQL", host = "postgis", dbname = "bu" ,
    user = "bu_user", password = "bu_user")

Lidar data

I already loaded 2m resolution lidar data for the grid squares containing Hengistbury head and Arne into the data base from this source http://environment.data.gov.uk/ds/survey/index.jsp#/survey

Lidar data is high resolution so it gets big quickly. The query below will pull out the data for our area by using the bounding box of the polygons we downloaded last time. It takes a few seconds to run.

dsm<- pgGetRast(conn, "dsm_arne_hh",boundary=arne)
dtm<- pgGetRast(conn, "dtm_arne_hh",boundary=arne)

dsm[dsm<0.1]<-NA
dtm[dtm<0.1]<-NA

Plot to check

Let’s just plot one of these out to see if everything is OK.

plot(dsm)
plot(arne,add=T)

So we can see that the data are nicely lined up, as we would expect because I knew that the Lidar layers were in EPSG 27700 CRS and we used the bounding box of the polygons to clip out the rasters from the larger data set stored on the server. This was nice and easy.

To get an idea of the size of the raster we could have a look at the number of rows and columns.

dsm@ncols
## [1] 1931
dsm@nrows
## [1] 2054

So we have 3966274 values. About 4 million. This is quite large, but it can be handled by R without any issues. If rasters are very large there can be problems holding all the data in memory and a different approach is needed. The problem of potentially slow raster processing applies to all GIS programs. Don’t try to process the entire Lidar data set for Dorset in one go! It can be possible to visualise this size dataset in a GIS using techniques such as pyramiding or using mosaics of tiles on a mapserver, but not all raster operations scale up easily. However in this case there are no problems.

The basic statistical properties of a raster can be summarised quickly in R

summary(dsm)
##               layer
## Min.          0.100
## 1st Qu.       0.957
## Median        4.115
## 3rd Qu.      12.438
## Max.         68.680
## NA's    1086336.000

So we have a minimum level for the surface below sea level (as Arne is on the coast) and a maximum height of just 68.7 meters. The NAs are missing values that occur in this case as the Lidar point cloud didn’t cover the water surface.

Terrain modelling

These layers are actually standard derived products that were built from the original lidar point cloud. This is much more convenient than building your own, which can be done. We have a digital terrain model that should represent the ground surface and a digital surface model that represents the top of vegetation and buildings. There is a lot we can do with these data.

For example we can form contours from the digital terrain model.

contours<-rasterToContour(dtm,nlevels=10)
plot(dtm)
plot(contours,add=T)

Terrain modelling derives additional layers from these building blocks by running algorithms over them. The algorithms look at the neighbourhood of each pixel in order to work out elements such as the slope and aspect. The raster package in R has a range of functions for this. Let’s get the slope and aspect from the digital terrain model. This makes more sense than fiding the slope from the surface model, as this would pick up edges of forest and building as slopes.

slope<-terrain(dtm, opt='slope', unit='degrees')
aspect<-terrain(dtm, opt='aspect',unit='degrees')
plot(slope)
plot(arne,add=T)

Raster algebra

Carrying out raster operations in R is very simple and inuitive. To obtain the height of vegetation (and maybe some buildings ) we just need to subract the digital terrain model from the digital surface model.

height<-dsm-dtm

We may want to set all values below a low threshold to NA

height[height<0.1]<-NA
plot(height)
plot(arne,add=T)

Raster aggregation

We often need to coarsen the resolution of a raster map and use some statistic calculated from all the pixels within a window. For example, if we are interested in the vegetation height around a study point we would not want to know the precices height at the point as it might just happen to fall into a small gap in a forest canopy. The orginal raster is at 2m resolution, so if we want mean vegetation height at 10m resolution we use factor of 5.

height10x10<-aggregate(height, fact=5, fun=mean, expand=TRUE, na.rm=TRUE)
height10x10[height10x10<0.1]<-NA
plot(height10x10)
plot(arne,add=T)

Saving the rasters

We might want to do some further work on the rasters in another desktop GIS program such as Arc or QGIS. They can be saved in conventional GeoTiff format locally.

writeRaster(dsm,"arne_dsm.tif",format="GTiff", overwrite=TRUE)
writeRaster(dtm,"arne_dtm.tif",format="GTiff", overwrite=TRUE)
writeRaster(height,"arne_height.tif",format="GTiff", overwrite=TRUE)
writeRaster(height10x10,"arne_height10.tif",format="GTiff", overwrite=TRUE)
writeRaster(slope,"arne_slope.tif",format="GTiff", overwrite=TRUE)
writeRaster(aspect,"arne_aspect.tif",format="GTiff", overwrite=TRUE)

Adding to the Leaflet map

Looking at raster layers as a static image without any reference point and without the ability to zoom in an out is beginning to become frustrating. A GIS such as Arc or QGIS has very good capabilities for visualising several layers and switching between them.

contours_4326<- spTransform(contours, CRS("+init=epsg:4326"))

library(leaflet)
cls<-terrain.colors(100)[100:1]

arne_map<-arne_map %>%
  addRasterImage(height,group="Height",colors=cls) %>%
  addRasterImage(height10x10,group="Height 10m",colors=cls) %>%
  addRasterImage(dtm,group="DTM",colors=cls) %>%
  addPolylines (data=contours_4326,group="Contours",color="black",weight = 1,opacity=1) %>% 
  removeLayersControl() %>%
  addLayersControl(baseGroups = c("OSM","Satellite"),overlayGroups = c("Arne SSSIs","Height","Height 10m","DTM","Contours"))
arne_map

So now we have the basis of a more complete site map. We’ll clean up and save the work for the next class.

Cleaning up

dbDisconnect(conn)
## [1] TRUE
save(list=ls(),file="spatial_class2.rob")