Pulling in OSMand data for the spatial sampling section

[pbkeating]

Hi!

Awesome work on the spatial sampling! Would it be possible to integrate pulling underlying OSMand data on households (where available) into the script based on the osmand package https://github.com/ropensci/osmdata ?

Suggestion from Jorieke, GIS specialist on the team. Thanks, P

[AmyMikhail]

Hi @pbkeating

Do you mean this:

Also note that with R you can download polygons of buildings directly from Open Street Maps, using the osmdata package. We have not done this in the case study but a chapter is being added to the epiRhandbook which will have a walk-through.

I don't think this has yet been added to the Epi R handbook, but I tagged @aspina7 as he might have more news on this.

In the meantime, there is a nice blog which shows how it can be done here.

The below code seems to work, but would be greate if someone (Joreike?) could check that it makes sense:

1. Add osmdata to the packages to load
2. After importing the shapefile for Kario camp, get the boundary box for it using sf::st_bbox()
3. Query open street map with opq() to get the geometries of the buildings that fall within the boundary box
4. The feature to fetch building geometries is: add_osm_feature(key = "building")
5. Convert the query output to an osmdata_sf object with osmdata_sf()
6. Note that the osm_points slot in this object are the points used to create the building polygons - they are NOT the building centroids.
7. To do further work on the building geometries, the osm_polygons slot of the query output needs to be converted to a simple features object with st_as_sf()
8. Note that some of the buildings returned are within the boundary box but outside the polygon that demarcates the borders of the camp. The geometries for these buildings can be removed with sf::st_intersection() specifying the building polygons as x (the data to be trimmed) and Kario camp polygon as y (the boundaries to use for the trimming).
9. Now we can use st_centroid() to get the centroid for each building polygon.
10. The building centroids can then be plotted with addCircleMarkers() onto the leaflet map.

# Load packages -----------------------------------------------------------

# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman", quiet = TRUE)

# Install and/or load required packages
pacman::p_load(
  leaflet, 
  mapedit, 
  sf,          
  osmdata,  
  here,         
  dplyr,        
  rio     
)

# Import shape files ------------------------------------------------------

# Import the polygon for Kario camp:
kario_polygon <- sf::read_sf(here::here("data", 
                                        "SDN_Kario_17Feb2019_polygon", 
                                        "SDN_Kario_17Feb2019_polygon.shp"))

# Get buildings from OSM --------------------------------------------------

# Get a boundary box for Kario camp from the shapefile:
kario_bbox <- sf::st_bbox(kario_polygon)

# Get the builing points within the boundary box:
households <- osmdata::opq(bbox = kario_bbox) %>%

  # Get the buildings in this boundary box:
  add_osm_feature(key = "building") %>%

  # Convert the building points to a simple feature object:
  osmdata::osmdata_sf() 

# Remove households that are in the boundary box but outside the polygon:
households_trimmed <- households$osm_polygons %>% 

  # Convert to a simple features object:
  sf::st_as_sf() %>% 

  # Trim to points which are inside the Kario camp polygon:
  sf::st_intersection(y = kario_polygon) %>% 

  # Get centroids (points) for each building polygon:
  mutate(points = st_centroid(geometry)) %>% 

  # Remove the empty ID column:
  select(-id)

# Create leaflet map ------------------------------------------------------

# Create a leaflet map of Kario camp:
kario_map <- leaflet() %>% 

  # Add the background tiles of your choice
  addProviderTiles("OpenStreetMap.HOT") %>% 

  # Plot the polygon of Kario camp ontop of the background map tiles:
  addPolygons(
    data  = kario_polygon, 
    color = "red",         
    fillOpacity = 0 
    ) %>% 

  # Plot the points representing buildings ontop:
  addCircleMarkers(
    data = households_trimmed$points, 
    radius = 5 
    )

# View the map:
kario_map

[AmyMikhail]

@aspina7 and @pbkeating note this process gives 2609 households, is that the expected number?