1 Get OSM data

1.1 Lima administrative boundary

OSM boundaries were downloaded via https://osm-boundaries.com/

adminLima <- st_read(dsn="data_day5/lima.gpkg", layer="limaAdmin")
## Reading layer `limaAdmin' from data source `/home/slautenb/Documents/lehre/HD/ws_2020_21/heikaLab/R4UrbanDataLab_2020/data_day5/lima.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -77.19921 ymin: -12.51993 xmax: -76.62082 ymax: -11.57244
## geographic CRS: WGS 84
limaBbox <- st_bbox(adminLima)
limaBboxStr <- paste(limaBbox, collapse = ", ")
tm_shape(adminLima) + tm_polygons(col="admin_level")

Since the airport is missing it might be useful, to combine the administartive boundary from Lima and Callao at least for some analysis:

adminLimaCallao <- st_read(dsn="data_day5/lima.gpkg", layer="limaCallao")
## Reading layer `limaCallao' from data source `/home/slautenb/Documents/lehre/HD/ws_2020_21/heikaLab/R4UrbanDataLab_2020/data_day5/lima.gpkg' using driver `GPKG'
## Simple feature collection with 2 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -77.25156 ymin: -12.51993 xmax: -76.62082 ymax: -11.57244
## geographic CRS: WGS 84
limaBbox <- st_bbox(adminLimaCallao)
limaBboxStr <- paste(limaBbox, collapse = ", ")
tm_shape(adminLimaCallao) + tm_polygons(col="name")

1.2 Highways for Lima

We will be using the ohsome api to query a number of features. We specify the following parameters: - bounding box - time in the YYYY-MM-DD format (I chose a relative new time stamp - the planet file is frequently updated but always a bit behind given the heavy lifitng involved) - a filter that defines which key-value combinations should be downloaded - with properties=“tags” I specify that I would like to get the tag (to be able to distinguish the differt road classes) - using clipGeometry=TRUE geometries are clipped at the bounding box. Otherwise longer roads might continue outside the bounding box region

For information on the definition of the different OSM keys and values refer to the OSM Wiki. For highway for example see: https://wiki.openstreetmap.org/wiki/Key:highway

limaHighway_result <- RCurl::postForm("https://api.ohsome.org/v1/elements/geometry", bboxes = limaBboxStr, filter = "highway in (motorway, trunk, primary, secondary, tertiary) and type:way", time = "2020-05-01",  binary=FALSE, properties="tags", clipGeometry =TRUE)

content <- geojsonio::as.json(limaHighway_result)
highwayLima <- geojsonio::geojson_sf(content)

Create an additional argument that will be used for plotting:

highwayLima <- highwayLima %>% 
    mutate( cls4map = case_when(
                highway == "trunk" | highway == "motorway" ~  10,
                highway == "primary"        ~   5,
                highway == "secondary"        ~   4,
                highway == "tertiary"        ~   3
))

1.3 Health sites

Healthsites in Lima. Healthsites might be represented by nodes, ways or relationship. Therefore, I use a different end point (centroid instead of geometry) that returns only the centroid of the feature. For relations that involve multiple buildings belonging to the same hospital that could lead to a point outside of the geometries.

limaHealth_result <- RCurl::postForm("https://api.ohsome.org/v1/elements/centroid", bboxes = limaBboxStr, filter = "amenity in (hospital, clinics, doctors, pharmacy) or healthcare in (centre, clinic, doctor, community_health_worker, hospital)", time = "2020-05-01",  binary=FALSE, properties="tags")

content <- geojsonio::as.json(limaHealth_result)
healthsitesLima <- geojsonio::geojson_sf(content)

This returned 2624 healthsites.

The different types are distributed as follows (using the two keys amenity and healthcare which can be both used in OSM).

xtabs(~amenity + healthcare, data=healthsitesLima, addNA = TRUE)
##           healthcare
## amenity    centre clinic doctor hospital pharmacy podiatrist  yes <NA>
##   clinic        2    103      0        1        0          0    0    0
##   doctors       0      0    114        0        0          1    1  249
##   hospital      0      0      0       36        0          0    0   94
##   pharmacy      0      0      0        0      121          0    0 1894
##   <NA>          6      0      1        1        0          0    0    0
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(highwayLima) + tm_lines(col= "highway", palette = "Set1") + tm_basemap("OpenStreetMap.DE") + tm_shape(healthsitesLima) + tm_dots(col="amenity", legend.show = TRUE, palette= "Accent")
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(highwayLima) + tm_lines(col= "highway", palette = c("blue", "red", "orange", "yellow", "blue" ), lwd= "cls4map", legend.lwd.show = FALSE ) +  tm_shape(healthsitesLima %>% filter(amenity == "hospital" | amenity=="clinic")) + tm_dots(col="amenity", legend.show = TRUE, palette= c("black", "green"))

2 Distance calculation

2.1 Distance to nearest hospital

First we creat an empty raster that is used as a template for the distance rasters (definining the extent, the cell size and the CRS of the raster).

limaRas <- raster(xmn= limaBbox[1], xmx=limaBbox[3], ymn= limaBbox[2], ymx= limaBbox[4], res=0.005)

distanceFromPoints calculates the distance to the closest point of the set of points provided. The coordinate system is WGS84 - distanceFromPoints calculates geodesic distance for geographic data. In adition I calculate isolines for a nicer cartographic representation.

dist2Hospital <- distanceFromPoints(limaRas, xy= healthsitesLima %>% filter(amenity == "hospital")) 
# mask by admin boundary
limaMask <- rasterize(adminLimaCallao, limaRas)
dist2Hospital <- mask(dist2Hospital, limaMask)

# create countour lines
dist2HospitalContour5km <- rasterToContour(dist2Hospital, levels= seq(0, 30000, by=5000))  %>%   st_as_sf()
dist2HospitalContour2_5km <- rasterToContour(dist2Hospital, levels= seq(2500, 30000, by=5000))  %>%   st_as_sf()
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(dist2Hospital) + tm_raster(col="layer", alpha=.5, palette = "-plasma", title = "Distance hospital [m]") +
  tm_shape(healthsitesLima %>% filter(amenity == "hospital")) + 
  tm_dots() +
  tm_shape(dist2HospitalContour5km) + tm_lines(lty=1, col="brown") +
  tm_shape(dist2HospitalContour2_5km) + tm_lines(lty=4, col="brown") +
  tm_shape(adminLima) + tm_borders()+
  tm_basemap("OpenStreetMap.DE")

Of course distances at the islande are not to sensitive. We see that most parts of the administrative district of Lima and Callao are within 5km euclidian distance to any hospital. However, some settlements are farer away than 10km from the closest hospital.

2.2 Distance to primary roads

To calculate distance to primary roads, we have to first rasterize the roads and when use distance to calculate the distance to all non NA cells. The selection of primary roads use filter from tidyvers as above.

# convert to raster
primaryRoadsRas <- highwayLima %>% 
  filter(highway== "primary") %>% 
  rasterize(y=limaRas)
# calc euclidian distance
dist2Primary <- distance(primaryRoadsRas) 
# mask areas outside Lima An Callao
dist2Primary <- mask(dist2Primary, limaMask)
# calculate isolines
dist2PrimaryContour1km <- rasterToContour(dist2Primary, levels= seq(1000, 30000, by=2000))  %>%   st_as_sf()
dist2PrimaryContour2km <- rasterToContour(dist2Primary, levels= seq(2000, 30000, by=2000))  %>%   st_as_sf()

Lets see what we got:

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(dist2Primary) + tm_raster(col="layer", alpha=.5, palette = "-plasma", title = "Distance primary road [m]", style="log10_pretty") +
  tm_shape(healthsitesLima %>% filter(amenity == "hospital")) + 
  tm_dots() +
  tm_shape(dist2PrimaryContour1km) + tm_lines(lty=1, col="brown") +
  tm_shape(dist2PrimaryContour2km) + tm_lines(lty=4, col="brown") +
  tm_basemap("OpenStreetMap.DE")
## Warning: non-rounded breaks occur, because style = "log10_pretty" is designed
## for large values

If we want to get statistics on how far the hospitals are from primary roads we can simply extract the information from the distance raster.

# create sp object for hospitals
hospitalSp <- healthsitesLima %>% filter(amenity == "hospital") %>% as("Spatial")
hospitalsDist2Primary <- raster::extract(dist2Primary, hospitalSp)
hist(x=hospitalsDist2Primary, main="Distance hospitals to primary roads in Lima", xlab= "Distance [m]", col="wheat", breaks=seq(0, 12000, by=500), las=1)

2.3 Population in distance of

popRas <- raster("data_day5/per_ppp_2020_UNadj_constrained.tif")
print(popRas)
## class      : RasterLayer 
## dimensions : 21977, 15214, 334358078  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : -81.33125, -68.65292, -18.35208, -0.03791647  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : /home/slautenb/Documents/lehre/HD/ws_2020_21/heikaLab/R4UrbanDataLab_2020/data_day5/per_ppp_2020_UNadj_constrained.tif 
## names      : per_ppp_2020_UNadj_constrained 
## values     : 0, 26471.09  (min, max)

Crop to case study region

popRas <- crop(popRas, adminLimaCallao)
names(popRas) <- "per_ppp_2020_UNadj_constrained"
tmap_mode("view")
tm_shape(popRas) + tm_raster(col="per_ppp_2020_UNadj_constrained", style = "log10_pretty", alpha=.7) + tm_scale_bar() + tm_layout(title="Population density - WorldPop constrained" ) + tm_basemap("OpenStreetMap.DE")
crs(popRas)
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

New template raster with same resolution as popRas and new mask and distance calculation to hospitals.

dist2Hospital <- distanceFromPoints(popRas, xy= healthsitesLima %>% filter(amenity == "hospital")) 
# mask by admin boundary
limaMask <- rasterize(adminLimaCallao, popRas)
dist2Hospital <- mask(dist2Hospital, limaMask)

# create countour lines
dist2HospitalContour5km <- rasterToContour(dist2Hospital, levels= seq(0, 30000, by=5000))  %>%   st_as_sf()
dist2HospitalContour2_5km <- rasterToContour(dist2Hospital, levels= seq(2500, 30000, by=5000))  %>%   st_as_sf()

Zonal statistics

Discretize the distance raster by cut

dist2HospitalZones <- cut(dist2Hospital, breaks=c(0, 500, 1000, 2500, 5000, 10000, 26000))
tm_shape(dist2HospitalZones) + tm_raster(col="layer")
popByDist2Hospital <- zonal(popRas, dist2HospitalZones, fun= 'sum')
popByDist2Hospital <- as.data.frame(popByDist2Hospital)
popByDist2Hospital$zoneName <- c("0-0.5", "0.5-1", "1-2.5", "2.5-5", "5-10", ">10")
# ensue factor levels are in right order
popByDist2Hospital$zoneName <- factor(popByDist2Hospital$zoneName, levels= c("0-0.5", "0.5-1", "1-2.5", "2.5-5", "5-10", ">10"))
ggplot(popByDist2Hospital, mapping=aes(x=zoneName, y=sum)) + geom_col()+ xlab("Distance to nearest hospital [km]") + ylab("Population") + labs("Population in distance to hospital", subtitle = "Lima, OSM and WorldPop constrained")

We could also display the distribution of values in each of the zones by using functionality from the rasterVis package that comes loaded with functions for visualisation of raster data.

require(rasterVis)
## Loading required package: rasterVis
## Loading required package: lattice
## Loading required package: latticeExtra
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
rasStack <- raster::stack(dist2HospitalZones, popRas)
names(rasStack) <- c("dist2HospitalZone", "pop")
bwplot(pop ~ dist2HospitalZone, data=rasStack, scales=list(y=list(log=10)))

2.4 Exercise

Calculate the distance to secondary roads and when calculate the distance at which all amenity=doctors are from secondary roads.

What is the distance of hospitals to either primary or secondary roads?

How are people distributed in distance to amenity=doctors?