Load packages

The giscourse package contains functions that use the pre loaded open street map derived pgrouting network held on the server.

library(giscourse)
library(sf)
library(mapview)
library(dplyr)
library(sf)
library(leaflet.extras)

Load data

The catchments shapefile consists of lines drawn between the starting point and the end point of the journeys.

catch<-read_sf("Catchments.shp")
#access<-read_sf("GB_AccessPoint.shp")
#access<-st_crop(access,catch)
parks<-read_sf("UGA_singlepart.shp")

Turning the lines into points

First transform the catchment layer to EPSG 4326 (lattitude and longitude).

catch<-st_transform(catch,4326)
catch$Id<-1:length(catch$Id) # Add an unique ID

Now, casting the line geometry to point produces two points for each id. So when the result is grouped by id and summarised it will be turned into a three column data frame with the id, start and end point.

catch_points<-st_cast(catch,"POINT")
catch_points %>% group_by(Id) %>% summarise(start=geometry[1],end=geometry[2]) ->routes

Finding the nodes on the pgrouting network

The giscourse::pnts_to_nodes function finds the nearest node on the network for each point. So we can find the end nodes and the start nodes.

conn<-sconnect11()
end_nodes<-giscourse::pnts_to_nodes(routes$end)
end_nodes$route<-1:299
start_nodes<-giscourse::pnts_to_nodes(routes$start)
start_nodes$route<-1:299
nodes<-rbind(start_nodes,end_nodes)
st_write(nodes,"nodes.shp",delete_dsn=TRUE)
## Deleting source `nodes.shp' using driver `ESRI Shapefile'
## Writing layer `nodes' to data source `nodes.shp' using driver `ESRI Shapefile'
## features:       598
## fields:         2
## geometry type:  Point

Finding the routes

The shortest route function can now be used to find 299 lines representing the shortest routes along the road network between the ends of the lines strings loaded from the catchment.shp file. The calculations take a few minutes so after running the code once it is commented out and the results saved to file.

# d<-data.frame(start=start_nodes$id,end=end_nodes$id)


# short<-shortest_route(c(d[1,1],d[1,2]))
# short$id<-1
# 
# for(i in 2:299){
#   sh<-shortest_route(c(d[i,1],d[i,2]))
#   sh$id<-i
#   short<-rbind(short,sh)
# }

# st_write(short,"shortest_routes_single.shp")

Aggregating and casting to multilinestring routes.

# library(dplyr)
# short %>% group_by(id) %>% summarise(dist=sum(length_m),time=sum(cost_s)) ->routes
# routes<-st_cast(routes,"MULTILINESTRING")
# st_write(routes,"shortest_routes_multi.shp",delete_dsn =TRUE)
short<-st_read("shortest_routes_single.shp")
## Reading layer `shortest_routes_single' from data source `/home/rstudio/webpages/examples/park_catchments/shortest_routes_single.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15569 features and 4 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -1.980732 ymin: 50.71315 xmax: -1.759583 ymax: 50.8467
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
routes<-st_read("shortest_routes_multi.shp")
## Reading layer `shortest_routes_multi' from data source `/home/rstudio/webpages/examples/park_catchments/shortest_routes_multi.shp' using driver `ESRI Shapefile'
## Simple feature collection with 299 features and 3 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: -1.980732 ymin: 50.71315 xmax: -1.759583 ymax: 50.8467
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

This provides networks of either individual lines strings or rather more usefully a layer consisiting of 299 routes.

A problem arises however. The end nodes and start nodes are arbitrary. They do not necessarily coincide with a park entrance and a residential address respectively. One way to work around this is to calculate how many times each node is duplicated in the data set. Nodes that are present more than five times are assumed to be the end nodes.

nodes %>% group_by(id) %>% summarise(n=n()) -> nodes
nodes %>% filter(n>5) ->endpnts

Now find the name of the parks that coincide with the end points. There are many polygons in the park layer without names, so remove these.

parks %>% filter(!is.na(Name)) -> parks

Now use a spatial join to link the attributes of the park layer to the attributes of the endpoints layer. A small buffer is needed to ensure an overlap.

endpnts<-st_join(endpnts,st_buffer(parks,0.002))

Now join the attributes of the endpoints of the routes to the routes themselves.

st_join(routes,endpnts) %>% group_by(Name) ->routes
mapview(routes,zcol="Name",burst=TRUE, hide = TRUE, legend=FALSE) + mapview(endpnts, cex=0.2) + mapview(parks, col.regions="darkgreen") ->mp

mp@map %>% addFullscreenControl()

Density of usage

library(RColorBrewer)
pal<-brewer.pal(12,"YlOrRd")
short<-st_cast(short,"LINESTRING")
short$id<-as.character(short$geometry)
short %>% group_by(id) %>% summarise(n=n()) ->sections
sections$id<-""
mapview(sections,zcol="n", color=pal) + mapview(parks,col.regions="darkgreen") + mapview(nodes,cex=0.1)
sections<-st_transform(sections,27700)
routes<-st_transform(routes,27700)
#st_write(sections,dsn="shapefiles/line_sections.shp")
#st_write(routes,dsn="shapefiles/routes.shp")
#system("zip -r shapefiles.zip shapefiles")
# path<-"~/webpages/examples/park_catchments/SHAPE_1/temp/UK_residential_population_2011_1_km"
# flnm<-dir(path,pattern="shp")
# #add_shp(flnm = flnm, pth = path, srid = 27700, tabnm = "uk_residential_pop_1km",db = "gis_course")
# 
# cn<-sconnect()
# 
# fl<-st_zm(st_read(cn,"uk_residential_pop_1km"))
# str(fl)
# st_crs(fl)<-st_crs(routes)
# fl<-st_crop(fl,routes)
# fl$pop<-as.numeric(as.character(fl$pop))
# mapview(fl,zcol="pop")