How can Postgis and R be used as a GIS?
Introduction
This is a quick demo showing how to extract data from postgis and run most GIS operations in R, ending up with a simple web map for a site.
Lidar data
A lot of lidar data in various forms and resolutions is available here.
http://environment.data.gov.uk/ds/survey/index.jsp#/survey?grid=SZ39
I downloaded 2m DSM and DTM data for several grid squares using the interface.
For each grid square I got 4 zip files. When each one is unzipped I got a set of files in *.asc format.
If students did this they would need to be put back together to form a mosaic. Then to work on a specific site they would need to clip a part out of the mosaic. So there would be a lot of steps needed before getting the data ready for work.
This is time consuming for students and messy. It would be far nicer just to be able to get a single joined up raster for a site in one single step. This is one thing postgis can help with.
If all the files are downloaded and unzipped into the same folder using commandline options they can then all be uploaded all together into a postgis database with a few lines of code. Data for the grid squares containing Hengistbury head and Arne are in the data base.
OSM Vector data
Likewise Open street map vector data can be downloaded and added in to the database in a fully automated way from http://download.geofabrik.de/europe/great-britain/england
Natural England layers for the whole of England can be downloaded from here and added using shp2pgsql, again in an automated way.
Packages used
library(rpostgis)
library(sp)
library(raster)
library(rgeos)
library(leaflet)
library(leaflet.extras)Connecting to the database
Once all the data are set up and loaded into the database making a connection to it in R makes it accessible with queries.
conn <- dbConnect("PostgreSQL", host = "postgis", dbname = "bu" ,
user = "docker", password = "docker")
## NOTE THAT A GENUINE INSTITUTIONAL DATABASE WOULD USE SECRET HIDDEN PASSWORDS FOR SECURITY!Extracting site data
We have the Natural England data for the whole of England. However we are working at a site level So we can pull in just the site if we know its name. This can be checked easily by looking at the database using another program such as pgadmin or QGIS that browses through all the tables with a GUI interface.
query<-"select * from sssi_units_england where sssi_name = 'Arne'"
arne <- pgGetGeom(conn, query=query)
plot(arne,axes=TRUE)
Now to get the lidar data for this area. We want the DSN and DTM data we uploaded for a wider area, clipped to the site. This is easy enough. Just pass the bounding box of the polygon as a boundary.
dsm<- pgGetRast(conn, "dsm-arne-hh",boundary=arne)
dtm<- pgGetRast(conn, "dtm-arne-hh",boundary=arne)Deriving new layers from the Lidar data
Notice how easy it is to derive slope, aspect and hillshade from the DTM with a line of R.
slope<-terrain(dtm, opt='slope', unit='degrees')
aspect<-terrain(dtm, opt='aspect')
hillshade<-hillShade(slope,aspect,normalize = TRUE)The syntax in R for basic raster algebra is very simple. For example we can get the (vegetation) height simply by subtracting the dtm from the dsm.
height<-dsm-dtm
height[height<0.1]<-NAForming contour lines
Contour lines are also easily produced.
dsm[dsm<0]<-NA
dtm[dtm<0]<-NA
contours<-rasterToContour(dtm,nlevels=10)Plot to check
plot(dtm)
plot(contours,add=T)
plot(height)
plot(arne,add=T)
plot(slope)
plot(arne,add=T)
plot(hillshade, col=grey(0:100/100))
Getting the OSM line data
We can download all the vector lines for paths, roads etc from OSM within a five km radius of the centre of the site with another query.
query<-"select osm_id, way geom FROM dorset_line
WHERE ST_DWITHIN(way::geography, ST_MAKEPOINT(-2.04818,50.68983)::geography,5000) "
paths<-pgGetGeom(conn, query=query)This is the result
plot(paths,axes =T)
Notice that the OSM data are held in ESPG 4326 while the Lidar and Natural England data were in EPSG 27700 (National grid). Not a problem. Both R and Postgis use the Proj4 library to transform between any CRS. There is also a problem with lines extending beyond the boundary. We’ll sort this out in a moment.
Transforming between CRSs
To demonstrate the transform in PostGIS I’ll download the arne data again, but this time in EPSG 4326 (Lat Long) so that it matches the paths data.
query<-"select gid, st_transform(geom,4326) geom from sssi_units_england where sssi_name = 'Arne'"
arne <- pgGetGeom(conn, query=query)
dbDisconnect(conn) ## Not going to need this again now## [1] TRUEarne@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0Vector intersection
Notice that the path lines extended beyond the boundaries of the SSSIs. For analysis and visualisation we only want the paths within the reserve. This can be sorted out by taking the intersection. Again, it would be possible in postgis or in R itself. We’ve already got the data so I’ll do it in R.
library(rgeos)
library(raster)
paths<-gIntersection(paths, arne, byid = T)
plot(arne,axes=T)
plot(paths,add=T)
Buffering
OK we now have the paths in EPSG 4326. However if we want to buffer around them we need a projected layer (i.e. use the British National Grid CRS). Again, this is not a problem. We’ve clipped the layers now, so we’ll do the transform in R rather than the data base. If we only want to look at areas 10 and 50 meters from the paths we can do this with a buffer.
path_proj <- spTransform(paths, CRS("+init=epsg:27700"))
pathbuf10<-gBuffer(path_proj,width=10)
pathbuf50<-gBuffer(path_proj,width=50)
plot(pathbuf50,axes=T)
Spatial overlay
One of the exercises in last years practical was to extract the characteristics such as slope and aspect at random points. In this case I’ll select points within the buffer around the paths.
randpoints <- spsample(pathbuf50, n=100, type='random')
plot(pathbuf50,axes=TRUE)
plot(randpoints,add=TRUE,pch=21,bg="red")
The extract function runs a spatial overlay of the points on the rasters.
This produces a quick summary of the data values at the random points within 50m from the path.
Vegetation height, elevation and slope at random points 50m from paths
d<-data.frame(vegheight=extract(height,randpoints),elev=extract(dtm,randpoints),slope=extract(slope,randpoints))
summary(d)## vegheight elev slope
## Min. : 0.112 Min. : 0.617 Min. : 0.4677
## 1st Qu.: 0.181 1st Qu.: 4.059 1st Qu.: 1.8851
## Median : 0.321 Median : 8.073 Median : 4.0902
## Mean : 2.002 Mean :11.076 Mean : 5.0394
## 3rd Qu.: 0.910 3rd Qu.:16.121 3rd Qu.: 6.6712
## Max. :17.395 Max. :46.693 Max. :20.3449
## NA's :40Transforming between CRSs
We can back transform the buffers into lat and long if we need to. The units used for the buffer remain as meters, not degrees.
pathbuf10 <- spTransform(pathbuf10, CRS("+init=epsg:4326"))
pathbuf50 <- spTransform(pathbuf50, CRS("+init=epsg:4326"))
contours<- spTransform(contours, CRS("+init=epsg:4326"))And we can also transform the rasters, although this step is not needed for the leaflet map in the case of rasters as they are transformed anyway when they are added to the map (This doesn’t occur with polygons)
dsm <- projectRaster(dsm, crs = pathbuf50@proj4string)
dtm <- projectRaster(dtm, crs = pathbuf50@proj4string)
height <- projectRaster(height, crs = pathbuf50@proj4string)Leaflet map
To design a really nice leaflet map takes a bit of time, as you need to think about colour schemes and possible legends etc. There are a lot of possible options. I’ll just quickly dump out the layers we have now with some of the more or less default (i.e ugly) colours to demonstrate.
Taking a little bit more time would find some nice looking palettes for them. Note that the layers you actually want to see can be chosen using the control on the right. The default makes them all visible at once, so unclick the ones you don’t want. Transparency etc can be set programaticaally, or if the map is in a shiny app sliders can be added. This all takes a bit of time to get right, but would be well worth doing if this were a real life cartographic project. Once the map is right it can be embedded in a web site.
The control on the right hand side switches layers on and off: The control box on the left takes the map full screen.
cls<-terrain.colors(100)[100:1]
library(leaflet)
map<-leaflet(arne) %>%
#setView(-2.04818, 50.68983) %>%
addTiles(group="osm") %>%
addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
addPolygons(data=arne,group="ssi") %>%
addPolylines (data=paths,group="paths",color="red",weight=3) %>%
addPolylines (data=contours,group="contours",color="black",weight = 1,opacity=1) %>%
addPolygons (data=pathbuf50,group="path_buffer50") %>%
addRasterImage(dtm,colors=cls,group="dtm") %>%
#addRasterImage(slope,colors= grey(0:100/100),group="slope") %>%
addRasterImage(height,colors=cls,group="height") %>%
addFullscreenControl() %>%
addLayersControl(baseGroups = c("osm","Satellite"),
overlayGroups = c("ssi","dtm","height","paths","path_buffer50","contours"),
options = layersControlOptions(collapsed = TRUE))
mapWriting out to Arc format
All the files stored in memory in R can be written out in shapefile format or in any raster format and loaded into Arc for further work in the desktop GIS. Or Arc itself could be connected to the postgis data base.
writeOGR(arne, getwd(), "arne", driver="ESRI Shapefile",overwrite_layer = TRUE)
writeOGR(paths, getwd(), "paths", driver="ESRI Shapefile",overwrite_layer = TRUE)
writeRaster(dsm, filename="dsm.tif", format="GTiff")
### etc ...Conclusion
The demo has shown how data can be loaded into R using simple queries in Postgis, reprojected and clipped to the boundaries of the reserve. All the usual operations can be carried out to extrect derived layers from the surface and terrain models. Leaflet allows the results to be displayed as an embeddable web map.