2  Data Acquisition

2.1 Online databases

Note

Several data sets are referenced by the ESoR (Environnement, Societies and Health Risk) research group here

2.1.1 On global scale (access with R package)

Since the appearance of the sf package, which has greatly contributed to the popularization of spatial data manipulation with R, many packages for making geographic data (geometries and/or attributes) available have been developed. Most of them are API packages that allow to query data made available on the Web, directly with R. This chapter presents a non-exhaustive list of them.

2.1.2 On regional scale

Several open data portals are particularly useful in the region, including Humanitarian Data Exchange and Open Development.

2.1.3 On national scale (Cambodia)

Cambodia administrative boundaries, divided in 4 levels:

  • level 0: country,
  • level 1: province / khaet and capital / reach thani,
  • level 2: municipality / district
  • level 3: commune / khum quarter / sangkat.

The contours maps (shapefiles) are made available from the Humanitarian Data Exchange (HDX), provided by OCHA (United Nations Offices for the Coordination of Humanitarian Affairs). These maps were originally produced by the Departement of geography of the Ministry of Land Management, Urbanization and Construction in 2008 and unofficially updated in 2014 according to sub-decrees on administrative modifications. They were provided by WFP - VAM unit Cambodia.

You can download these administrative boundaries, as zip folders, here:

Population data:

Population data is available at these different levels from the Humanitarian Data Exchange (HDX) repository. It comes from the Commune database (CDB), provided by the Cambodia Ministry of Planning.

https://data.humdata.org/dataset/cambodia-population-statistics

Health Facility data:

The Humanitarian Data Exchange (HDX) repository provides a dataset on the location of health facilities (Referral Hospitals, Health Centers, Health Posts). These maps were originally produced by the Cambodia Ministry of Health (MoH).

https://data.humdata.org/dataset/cambodia-health

Transportation data:

The roads network is available from Humanitarian Data Exchange (HDX) repository. These maps were originally produced by the Cambodia Department of Geography of the Ministry of Land Management, Urbanization and Construction. They include: National road primary and secondary, Provincial road primary, Provincial and rural roads, Foot path, Cart track, Bridge line.

https://data.humdata.org/dataset/cambodia-roads

Hydrology data:

The hydrological network is available from Humanitarian Data Exchange (HDX) repository. These maps were originally produced by the Cambodia Department of Geography of the Ministry of Land Management, Urbanization and Construction. They include: rivers (“Non-Perenial/Intermittent/Fluctuating” and “Perennial/Permanent”), lakes

https://data.humdata.org/dataset/cambodia-water-courses-0

Digital Elevation Model (DEM):

The SRTM (Shuttle Radar Topography Mission) is a free DEM provided by NASA and NGA (formerly NIMA). Space Shuttle Endeavour (STS-99) collected these altimetry data during an 11-day mission in February 2000 at an altitude of 233 km using radar interferometry. The SRTM covers nearly 80% of the land area from 56° South latitude to 60° North latitude. Spatial resolution is approximately 30 meters on the line of the Equator.

The SRTM data can be downloaded here: http://srtm.csi.cgiar.org

2.2 OpenStreetMap

OpenStreetMap (OSM) is a participatory mapping project that aims to built a free geographic database on a global scale. OpenStreetMap lets you view, edit and use geographic data around the world.

Terms of use

OpenStreetMap is open data : you are free to use it for ant purpose as long as you credit OpenStreetMap and its contributers. If you modify or rely data in any way, you may distribute the result only under the same license. (…)

Contributors

(…) Our contributors incloude enthusiastic mapmakers, GIS professional, engineers running OSM servers, humanitarians mapping disaster-stricken areas and many mmore.(…)

2.2.1 Display and interactive map

The two main packages that allow to display as interactive map based on OSM are leaflet (Cheng, Karambelkar, and Xie 2022) and mapview (Appelhans et al. 2022).

2.2.1.1 leaflet

leaflet uses the javascript library Leaflet (Agafonkin 2015) to create interactive maps.

library(sf)
library(leaflet)

district <- st_read("data_cambodia/cambodia.gpkg", layer = "district", quiet = TRUE)
hospital <- st_read("data_cambodia/cambodia.gpkg", layer = "hospital", quiet = TRUE)


banan <- district[district$ADM2_PCODE == "KH0201", ]     #Select one district (Banan district: KH0201)
health_banan <- hospital[hospital$DCODE == "201", ]      #Select Health centers in Banan

banan <- st_transform(banan, 4326)                       #Transform coordinate system to WGS84
health_banan <- st_transform(health_banan, 4326)

banan_map <- leaflet(banan) %>%                          #Create interactive map
  addTiles() %>%
  addPolygons() %>%
  addMarkers(data = health_banan)
banan_map

Website of leaflet
Leaflet for R

2.2.1.2 mapview

mapview relies on leaflet to create interactive maps, its use is easier and its documentation is a bit dense.

library(mapview)
mapview(banan) + mapview(health_banan)

Website of mapview
mapview

2.2.2 Import basemaps

The package maptiles (Giraud 2021) allows downlaoding and displaying raster basemaps.
The function get_tiles() allow you to download OSM background maps and the function plot_tiles() allows to display them.
Renders are better if the input data used the same coordinate system as the tiles (EPSG:3857).

library(sf)
library(maptiles)
district <- st_read("data_cambodia/cambodia.gpkg", layer = "district", quiet = TRUE)
district <- st_transform(district, 3857)
osm_tiles <- get_tiles(x = district, zoom = 10, crop = TRUE)
plot_tiles(osm_tiles)
plot(st_geometry(district), border = "grey20", lwd = .7, add = TRUE)
mtext(side = 1, line = -2, text = get_credit("OpenStreetMap"), col="tomato")

2.2.3 Import OSM data

2.2.3.1 osmdata

The package osmdata (Padgham et al. 2017a) allows extracting vector data from OSM using the Overpass turbo API.

library(sf)
library(osmdata)
library(sf)

country <- st_read("data_cambodia/cambodia.gpkg", layer = "country", quiet = TRUE)
ext <- opq(bbox = st_bbox(st_transform(country, 4326)))                    #Define the bounding box
query <- add_osm_feature(opq = ext, key = 'amenity', value = "hospital")   #Health Center Extraction
hospital <- osmdata_sf(query)
hospital <- unique_osmdata(hospital)                                       #Result reduction (points composing polygon are detected)

The result contains a point layer and a polygon layer. The polygon layer contains polygons that represent hospitals. To obtain a coherent point layer we can use the centroids of the polygons.

Spherical geometry (s2) switched off
hospital_point <- hospital$osm_points
hospital_poly <- hospital$osm_polygons                                                             #Extracting centroids of polygons
hospital_poly_centroid <- st_centroid(hospital_poly)

cambodia_point <- intersect(names(hospital_point), names(hospital_poly_centroid))                  #Identify fields in Cambodia boundary
hospitals <- rbind(hospital_point[, cambodia_point], hospital_poly_centroid[, cambodia_point])     #Gather the 2 objects

Result display

library(mapview)
mapview(country) + mapview(hospitals)

Website of osmdata
osmdata

2.2.3.2 osmextract

The package osmextract (Gilardi and Lovelace 2021) allows to extract data from an OSM database directly. This package make it possible to work on very large volumes of data.

Website of osmextract
osmextract

For administrative boundaries, check here the administrative levels by country:

library(osmextract)
library(mapsf)
province <- oe_get(
  place = "Cambodia",
  download_directory = "data_cambodia/",
  layer = "multipolygons",
  extra_tags = c("wikidata", "ISO3166-2", "wikipedia", "name:en"),
  vectortranslate_options = c(
    "-t_srs", "EPSG:32648",
    "-nlt", "PROMOTE_TO_MULTI",
    "-where", "type = 'boundary' AND boundary = 'administrative' AND admin_level = '4'"
  ))
0...10...20...30...40...50...60...70...80...90...100 - done.
Reading layer `multipolygons' from data source 
  `/home/lucas/Documents/ForgeIRD/rspatial-for-onehealth/data_cambodia/geofabrik_cambodia-latest.gpkg' 
  using driver `GPKG'
Simple feature collection with 25 features and 29 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 211418.1 ymin: 1047956 xmax: 784614.9 ymax: 1625621
Projected CRS: WGS 84 / UTM zone 48N
mf_map(x = province)

roads <- oe_get(
  place = "Cambodia",
  download_directory = "data_cambodia/",
  layer = "lines",
  extra_tags = c("access", "service", "maxspeed"),
  vectortranslate_options = c(
  "-t_srs", "EPSG:32648",
  "-nlt", "PROMOTE_TO_MULTI",
  "-where", "
    highway IS NOT NULL
    AND
    highway NOT IN (
    'abandonded', 'bus_guideway', 'byway', 'construction', 'corridor', 'elevator',
    'fixme', 'escalator', 'gallop', 'historic', 'no', 'planned', 'platform',
    'proposed', 'cycleway', 'pedestrian', 'bridleway', 'footway',
    'steps', 'path', 'raceway', 'road', 'service', 'track'
    )
    "
),
  boundary = subset(province, name_en == "Phnom Penh"),
  boundary_type = "clipsrc"
)
0...10...20...30...40...50...60...70...80...90...100 - done.
Reading layer `lines' from data source 
  `/home/lucas/Documents/ForgeIRD/rspatial-for-onehealth/data_cambodia/geofabrik_cambodia-latest.gpkg' 
  using driver `GPKG'
Simple feature collection with 18794 features and 12 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 469524.2 ymin: 1263268 xmax: 503494.3 ymax: 1296780
Projected CRS: WGS 84 / UTM zone 48N
mf_map(x = roads)

2.3 Import from lat / long file

The function st_as_sf() makes it possible to transform a data.frame container of geographic coordinates into an object sf. Here we use the data.frame places2 created in the previous point.

library(sf)
place_sf <- st_as_sf(read.csv("data_cambodia/adress.csv"), coords = c("long", "lat"), crs = 4326)
place_sf
Simple feature collection with 2 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 104.8443 ymin: 11.54366 xmax: 104.9047 ymax: 11.55349
Geodetic CRS:  WGS 84
                                                 address
1 Phnom Penh International Airport, Phnom Penh, Cambodia
2 Khmer Soviet Friendship Hospital, Phnom Penh, Cambodia
                   geometry
1 POINT (104.8443 11.55349)
2 POINT (104.9047 11.54366)

To create a sf POINT type object with only one pair of coordinate (WGS84, longitude=0.5, latitude = 45.5) :

library(sf)
test_point <- st_as_sf(data.frame(x = 0.5, y = 45.5), coords = c("x", "y"), crs = 4326)
test_point
Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0.5 ymin: 45.5 xmax: 0.5 ymax: 45.5
Geodetic CRS:  WGS 84
          geometry
1 POINT (0.5 45.5)

We can display this object sf on an OpenStreetMap basesmap with the package maptiles maptiles (Giraud 2021).

library(maptiles)
osm <- get_tiles(x = place_sf, zoom = 12)
plot_tiles(osm)
plot(st_geometry(place_sf), pch = 2, cex = 2, col = "red", add = TRUE)

2.4 Geocoding

Serveral pakages alow you to geocode addresses. The package tidygeocoder (Cambon et al. 2021) allow the use of a large number of online geocoding sevices. The package banR (Gombin and Chevalier 2022), which is based on the National Address Base, is the particularly suitable for geocoding addresses in France.

2.4.1 tidygeocoder

library(tidygeocoder)
test_adresses <- data.frame(
  address = c("Phnom Penh International Airport, Phnom Penh, Cambodia",
              "Khmer Soviet Friendship Hospital, Phnom Penh, Cambodia"))
places1 <- geocode(test_adresses, address)
places1
# A tibble: 2 × 3
  address                                                  lat  long
  <chr>                                                  <dbl> <dbl>
1 Phnom Penh International Airport, Phnom Penh, Cambodia  11.6  105.
2 Khmer Soviet Friendship Hospital, Phnom Penh, Cambodia  11.5  105.

Website by tidygeocoder :
tidygeocoder

2.4.2 banR (Base Adresse Nationale)

# remotes::install_github("joelgombin/banR")
library(banR)
mes_adresses <- data.frame(
  address = c("19 rue Michel Bakounine, 29600 Morlaix, France",
              "2 Allee Emile Pouget, 920128 Boulogne-Billancourt")
)
places2 <- geocode_tbl(tbl = mes_adresses, adresse = address)
places2
# A tibble: 2 × 18
  address        latit…¹ longi…² resul…³ resul…⁴ resul…⁵ resul…⁶ resul…⁷ resul…⁸
  <chr>            <dbl>   <dbl> <chr>     <dbl> <chr>   <chr>   <chr>   <chr>  
1 19 rue Michel…    48.6   -4.32 19 All…    0.17 housen… 29124_… 19      Allée …
2 2 Allee Emile…    48.7    2.31 2 Squa…    0.2  housen… 92002_… 2       Square…
# … with 9 more variables: result_street <chr>, result_postcode <chr>,
#   result_city <chr>, result_context <chr>, result_citycode <chr>,
#   result_oldcitycode <chr>, result_oldcity <chr>, result_district <chr>,
#   result_status <chr>, and abbreviated variable names ¹​latitude, ²​longitude,
#   ³​result_label, ⁴​result_score, ⁵​result_type, ⁶​result_id,
#   ⁷​result_housenumber, ⁸​result_name

Website of banR :
An R client for the BAN API

2.5 Digitization

The package mapedit (Appelhans, Russell, and Busetto 2020) allows you to digitize base map directly in R. Although it can be practical in some cases, in package cannot replace the functionalities of a GIS for important digitization tasks.

Gif taken from mapedit website