The master map output is certainly not the easiest to use. The main problem is that it seems to be designed with cartography in mind, rather than analysis. The categorisation of te objects is rather confusing. One aspect that makes filtering the data especially difficult is that several fields hold arrays of options, rather than one single category. This can be handled using a PostGIS query, but it makes life challenging.
The following code loads in a lot of potentially useful libraries. Not all will actually be used.
library(rpostgis)
library(RPostgreSQL)
library(sp)
library(raster)
library(rgeos)
library(rgdal)
library(RColorBrewer)
library(leaflet.extras)
library(mapview)
library(sf)
library(DT)
dt<-function(x) DT::datatable(x,
extensions = c('Buttons'), options = list(
dom = 'Blfrtip',
buttons = c('copy', 'csv', 'excel'), colReorder = FALSE
))
mapviewOptions(basemaps =c("OpenStreetMap","Esri.WorldImagery", "OpenTopoMap", "OpenStreetMap.BlackAndWhite", "OpenStreetMap.HOT"))
conn <- dbConnect("PostgreSQL", host = "postgis", dbname = "bu" ,
user = "docker", password = "docker")
The code runs a postGIS subquery to find the distances based on the criteria used to filter the mastermap layers.
For example the line
Will take all lines from the os2008 totographicline layer where the phrase “Road Or Track” appears in the descritpivegroup array. I am unsure whether this is the best way of selecting the objects that you want at the moment. The code can be modified to use different criteria (don’t try this without guidance, as any slight error in the syntax will break the query. I always test the queries in PGAdmin before adding them to R). The query takes about half a minute to run and adds four fields to the nest sites table.
query<-"SELECT
tid,year,\"NEST\",\"OUTCOME\",geometry geom,
(SELECT
st_distance(wkb_geometry,geometry)
FROM
(select * from os2008.topographicline where 'Road Or Track' = any(descriptivegroup) ) t
ORDER BY
geometry <-> wkb_geometry
LIMIT
1
) AS roaddistance2008,
(SELECT
st_distance(wkb_geometry,geometry)
FROM
(select * from os2016.topographicline where 'Road Or Track' = any(descriptivegroup) ) t
ORDER BY
geometry <-> wkb_geometry
LIMIT
1
) AS roaddistance2016,
(SELECT
st_distance(wkb_geometry,geometry)
FROM
(select * from os2008.topographicarea where 'Building' = any(descriptivegroup) ) t
ORDER BY
geometry <-> wkb_geometry
LIMIT
1
) AS builddistance2008,
(SELECT
st_distance(wkb_geometry,geometry)
FROM
(select * from os2016.topographicarea where 'Building' = any(descriptivegroup) ) t
ORDER BY
geometry <-> wkb_geometry
LIMIT
1
) AS builddistance2016
FROM
curlews.nestsites"
nest_sites <- st_read(conn, query=query)
An alternative is to use the open street map layer. These are the features that show up on the open street map base layer.
query<-"drop table if exists curlews.nestsites_lat_lon;
create table curlews.nestsites_lat_lon as select tid,year,\"NEST\",\"OUTCOME\",st_transform(geometry,4326) geom from curlews.nestsites;
SELECT
tid,year,\"NEST\",\"OUTCOME\",geom,
(SELECT
st_distance(way::geography,geom::geography)
FROM
(select * from england_roads ) t
ORDER BY
geom<-> way
LIMIT
1
) AS osm_distance
FROM
(select * from curlews.nestsites_lat_lon) s"
osm <- st_read(conn, query=query)
osm$osm_distance<-round(osm$osm_distance,0)
nest_sites$roaddistance2008<-round(nest_sites$roaddistance2008,0)
nest_sites$roaddistance2016<-round(nest_sites$roaddistance2016,0)
nest_sites$builddistance2008 <-round(nest_sites$builddistance2008,0)
nest_sites$builddistance2016 <-round(nest_sites$builddistance2016,0)
The map can be opened full screen by clicking on the box at the top. The controls can show several different base maps and switch from the “nest_sites” points, which have distances calculated using os master map to the open street map (osm) layer. There is a distance and area calculating tool included.
mp<- mapview(nest_sites) + mapview(osm)
tile<-"OpenStreetMap" ## Change if needed
zoom<- - 3 ## Larger negative values to zoom out
pos<-"bottomleft"
mp@map %>% addFullscreenControl() %>%
addMiniMap(position=pos,zoomLevelOffset = zoom,tiles=tile,toggleDisplay=TRUE,minimized=TRUE) %>%
addMeasurePathToolbar() %>%
addMeasure( primaryLengthUnit = "meters", primaryAreaUnit = "hectares")
Download the results by clicking the buttons if you want them in a plain Excel spreadheet.
nest_sites <- nest_sites %>% st_set_geometry(NULL)
osm<-osm %>% st_set_geometry(NULL)
osm<-osm[,c(1,5)]
d<-merge(nest_sites,osm)
dt(d)
dbDisconnect(conn)
## [1] TRUE