Robinlovelace
commented
7 months ago Here's some pavements code, a bit rough round the edges...
---
format: gfm
# jupyter: python3
---
```{r}
#| label: setup
#| message: false
library(tidyverse)
library(sf)
```
The input dataset is OS's MM Topo product, provided in a zip file as follows:
First we'll check the releases and create a new one:
```{r}
#| eval: false
system("gh release list")
# system("gh release create mm --title 'MM Topo' --notes 'MM Topo release'")
# f = list.files("~/Downloads", pattern = "Download_ed", full.names = TRUE)
# file.copy(f, "mm_topo.zip")
# Download the mm_topo.zip with gh release download
system("gh release download mm --pattern 'mm_topo.zip'")
fs::file_size("mm_topo.zip")
unzip("mm_topo.zip")
```
```{r}
dir_mm = list.dirs(recursive = FALSE) |>
data.frame(directory = _) |>
filter(grepl("mastermap", directory)) |>
pull(directory)
f_gpkg = list.files(dir_mm, pattern = "gpkg", full.names = TRUE)
```
Let's read the data:
```{r}
layers = st_layers(f_gpkg)
layers
mm_ed_area = read_sf(f_gpkg, layer = "Topographicarea")
table(mm_ed_area$descriptivegroup)
names(mm_ed_area)
pavements = mm_ed_area |>
filter(stringr::str_detect(descriptivegroup, "Roadside")) |>
filter(make == "Manmade")
road_or_track = mm_ed_area |>
filter(stringr::str_detect(descriptivegroup, "Road Or Track")) |>
filter(make == "Manmade")
pavements_1 = pavements |>
slice(1)
plot(pavements_1$geom)
paste(unlist(pavements_1[1, ] |> sf::st_drop_geometry()))
names(pavements)
```
# Describe other features available from OS MM Topo
- Traffic calming?
- Road polygons?
# Import OS Open Roads data
```{r}
#| eval=FALSE
if (!file.exists("simplified_network.geojson")) {
OS_road = sf::read_sf("../npt/outputdata/2024-04-09/edinburgh_and_lothians/simplified_network.geojson")
sf::write_sf(OS_road, "simplified_network.geojson")
}
OS_road = sf::read_sf("simplified_network.geojson")
# select a random road
OS_road_1 = OS_road |>
slice(1)
plot(OS_road_1$geometry)
```
We'll create a buffer around the roads with a width of 20m:
```{r}
#| eval=FALSE
OS_road_projected = OS_road |>
sf::st_transform("EPSG:27700")
pavements_bb = stplanr::bb2poly(pavements)
OS_road_projected = OS_road_projected[pavements_bb, ]
OS_road_buffer = OS_road_projected |>
sf::st_buffer(20, endCapStyle = "FLAT")
OS_roads_length = sf::st_length(OS_road_projected)
OS_roads_max = OS_road_projected |>
slice_max(OS_roads_length)
OS_road_buffer_1 = OS_road_buffer |>
slice_max(OS_roads_length)
OS_road_1 = OS_road_projected |>
slice_max(OS_roads_length)
# Find the pavements that intersect
pavements_intersection = pavements[OS_road_buffer_1, ]
plot(pavements_intersection$geom)
plot(OS_road_buffer_1$geometry)
plot(OS_road_1$geometry, col = "grey", add = TRUE)
plot(pavements_intersection$geom, col = "red", add = TRUE)
mapview::mapview(pavements_intersection) +
mapview::mapview(OS_road_buffer_1) +
mapview::mapview(OS_road_1)
pavements_intersection_within = sf::st_intersection(
pavements_intersection,
OS_road_buffer_1
)
mapview::mapview(pavements_intersection_within) +
mapview::mapview(OS_road_buffer_1) +
mapview::mapview(OS_road_1)
# Calculate average width
p1_width = sf::st_area(pavements_intersection_within) |>
sum() |>
as.numeric()
length_section = sf::st_length(OS_road_1) |>
as.numeric()
average_width = p1_width / length_section
average_width
```
# Calculate widths of pavements and carriageways associated with roads
```{r}
#| eval=FALSE
library(sf)
library(dplyr)
library(units)
library(ggplot2)
# Read and transform GeoJSON files
roads = st_read("data/road_example.geojson") |>
st_transform("EPSG:27700")
roadside_polygons = st_read("data/roadside_example.geojson") |>
st_transform("EPSG:27700") # Ensure both datasets are in the same CRS
# Create buffer zones for roads
roads_buffered = st_buffer(roads, dist = 20, endCapStyle = "FLAT")
# Identify and process intersections
intersections = st_intersects(roads_buffered, roadside_polygons, sparse = FALSE)
# Initialize a vector to store calculated average widths
average_widths = numeric(length = nrow(roads))
nrow(roads)
system.time({
for (i in seq(nrow(roads))) {
# for (i in 41:49) {
relevant_indices = which(intersections[i, ])
plot(roads_buffered$geometry[i])
# Add the road geometry:
plot(roads$geometry[i], col = "grey", add = TRUE)
plot(roadside_polygons$geometry[relevant_indices], col = "red", add = TRUE)
if (length(relevant_indices) > 0) {
relevant_polygons = roadside_polygons[relevant_indices, ]
pavements_intersection_within = st_intersection(relevant_polygons, roads_buffered[i, ])
intersecting_area = sum(st_area(pavements_intersection_within))
road_length = st_length(roads[i, ])
# Convert road_length to numeric for comparison
numeric_road_length = as.numeric(road_length, "meters")
if (numeric_road_length > 0) {
average_widths[i] = as.numeric(intersecting_area / road_length, "meters")
} else {
average_widths[i] = NA # Assign NA for cases where road length is zero or invalid
}
} else {
average_widths[i] = NA # Assign NA where no polygons intersect the buffer
}
}
})
# Add average_widths to the roads geojson
roads$average_width = average_widths
# mapview::mapview(roads_buffered) +
mapview::mapview(roads, zcol = "average_width") +
mapview::mapview(roadside_polygons)
# mapview::mapview(roads_buffered[42, ]) # single example
average_widths[42]
ggplot() +
geom_sf(data = roadside_polygons, fill = "green", color = NA, alpha = 0.5) +
geom_sf(data = roads, aes(color = average_width), size = 2) +
scale_color_viridis_c(option = "C", na.value = "grey50", guide = "colourbar", name = "Average Width (m)") +
labs(title = "Road Network by Average Roadside Width") +
theme_minimal()
```
```{r}
#| eval=FALSE
library(rgdal)
library(sf)
library(dplyr)
library(ggplot2)
library(geos)
# Read and transform GeoJSON files using sf
roads_sf = st_read("data/road_example.geojson") |>
st_transform("EPSG:27700")
roadside_polygons_sf = st_read("data/roadside_example.geojson") |>
st_transform("EPSG:27700")
# Convert sf objects to geos geometries for geos operations
roads_geos = as_geos_geometry(roads_sf)
roadside_polygons_geos = as_geos_geometry(st_geometry(roadside_polygons_sf))
# # Create buffered roads using geos
roads_buffered = geos_buffer(roads_geos, dist = 20, params = geos_buffer_params(end_cap_style = c( "flat")))
# Initialize average widths
average_widths = numeric(length = length(roads_geos))
system.time({
for (i in seq_along(roads_geos)) {
# Intersect buffered road with roadside polygons
intersections = geos_intersects(roads_buffered[i], roadside_polygons_geos)
relevant_indices = which(intersections)
plot(roads_buffered[i])
# Add the road geometry:
plot(roads_geos[i], col = "grey", add = TRUE)
plot(roadside_polygons_geos[relevant_indices], col = "red", add = TRUE)
if (length(relevant_indices) > 0) {
relevant_polygons = roadside_polygons_geos[relevant_indices]
pavements_intersection_within = geos_intersection(relevant_polygons, roads_buffered[i])
plot(pavements_intersection_within)
intersecting_area = sum(geos_area(pavements_intersection_within))
road_length = geos_length(roads_geos[i])
# Check road_length > 0 to avoid division by zero
if (road_length > 0) {
average_widths[i] = intersecting_area / road_length
} else {
average_widths[i] = NA # Assign NA for cases where road length is zero or invalid
}
} else {
average_widths[i] = NA # Assign NA where no polygons intersect the buffer
}
}
})
# Convert geos geometries back to sf for attaching average widths and further processing
roads_sf$average_width = average_widths
# Proceed with sf-based plotting or further analysis
# mapview::mapview(roads_buffered) +
mapview::mapview(roads_sf, zcol = "average_width") +
mapview::mapview(roadside_polygons_sf)
# mapview::mapview(roads_buffered[42, ]) # single example
average_widths[42]
ggplot() +
geom_sf(data = roadside_polygons_sf, fill = "green", color = NA, alpha = 0.5) +
geom_sf(data = roads_sf, aes(color = average_width), size = 2) +
scale_color_viridis_c(option = "C", na.value = "grey50", guide = "colourbar", name = "Average Width (m)") +
labs(title = "Road Network by Average Roadside Width") +
theme_minimal()
```
Prepare the sample data in edinburgh:
```{r}
#| eval=FALSE
dir_mm = list.dirs(recursive = FALSE) |>
data.frame(directory = _) |>
filter(grepl("mastermap", directory)) |>
pull(directory)
f_gpkg = list.files(dir_mm, pattern = "gpkg", full.names = TRUE)
layers = st_layers(f_gpkg)
layers
mm_ed_area = read_sf(f_gpkg, layer = "Topographicarea")
pavements = mm_ed_area |>
filter(stringr::str_detect(descriptivegroup, "Roadside")) |>
filter(make == "Manmade")
edinburgh = sf::st_read("data/simplified_network.geojson") |>
sf::st_transform("EPSG:27700")
# create bbbox of pavements and then find edinburgh within
bb_pavements = sf::st_bbox(pavements)
bb_polygon <- st_as_sfc(st_bbox(bb_pavements))
edinburgh = edinburgh |>
sf::st_intersection(bb_polygon) |>
sf::st_transform(st_crs(pavements))
# save edinburgh as geojson
sf::st_write(edinburgh, "data/edinburgh.geojson")
sf::st_write(pavements, "data/pavements.geojson")
```
```{r}
source("R/calculateAverageWidths.R")
road_geojson_path = "data/edinburgh.geojson"
roadside_geojson_path = "data/pavements.geojson"
system.time({
roads_sf = calculateAverageWidths(road_geojson_path, roadside_geojson_path, crs = "EPSG:27700", buffer_dist = 15)
})
# save roads_sf
sf::st_write(roads_sf, "data/roads_sf_width_2.geojson", delete_dsn = TRUE)
ggplot2::ggplot() +
geom_sf(data = roads_sf, aes(color = average_widths_right), size = 2) +
scale_color_viridis_c(option = "C", na.value = "grey50", guide = "colourbar", name = "Average Width (m)") +
labs(title = "Road Network by Average Roadside Width") +
theme_minimal()
summary(roads_sf$average_width)
summary(roads_sf$average_widths_right)
summary(roads_sf$average_widths_left)
```
Heads-up @jananiiU FYI here's some code...