library(sf)
library(leaflet)
<- st_read("data_cambodia/cambodia.gpkg", layer = "district", quiet = TRUE)
district <- st_read("data_cambodia/cambodia.gpkg", layer = "hospital", quiet = TRUE)
hospital
<- district[district$ADM2_PCODE == "KH0201", ] #Select one district (Banan district: KH0201)
banan <- hospital[hospital$DCODE == "201", ] #Select Health centers in Banan
health_banan
<- st_transform(banan, 4326) #Transform coordinate system to WGS84
banan <- st_transform(health_banan, 4326)
health_banan
<- leaflet(banan) %>% #Create interactive map
banan_map addTiles() %>%
addPolygons() %>%
addMarkers(data = health_banan)
banan_map
2 Data Acquisition
2.1 Online databases
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.
rnaturalearth
(South 2017): retrieves Natural Earth map data.
gadmr
(Guevarra 2021): retrieves data from the GADM (national and sub-national administrative divisions of all countries in the world).
rgeoboundaries
(Dicko 2021) : R client for the geoBoundaries API, providing political administrative boundaries of countries.cshapes
(Weidmann, Schvitz, and Girardin 2021): makes available national boundaries, from 1886 to present.
osmextract
(Gilardi and Lovelace 2021): allows importing OpenStreetMap data.osmdata
(Padgham et al. 2017b): to download and use OpenStreetMap data.
maptiles
(Giraud 2021) : This package downloads, composes and displays tiles from a large number of providers (OpenStreetMap, Stamen, Esri, CARTO or Thunderforest).
geonames
(Rowlingson 2019) : allows you to query the geonames DB, which provides locations in particular.wbstats
(wbstats2020?) andWDI
(R-WDI?): provide access to World Bank data and statistics.
sen2r
(R-sen2r?): allows automatic download and preprocessing of Sentinel-2 satellite data.
MODIStsp
(MODIStsp2016?): find, download and process MODIS images.
geodata
(R-geodata?): provides access to data on climate, elevation, soil, species occurrence and administrative boundaries.
elevatr
(R-elevatr?): provides access to elevation data made available by Amazon Web Services Terrain Tiles, the Open Topography Global Datasets API and the USGS Elevation Point Query Service.rgee
(R-rgee?): allows use of the Google Earth Engine API, a public data catalog and computational infrastructure for satellite images.
nasapower
(nasapower2018?): NASA client API (global energy resource forecasting, meteorology, surface solar energy, and climatology).geoknife
(geoknife2015?): allows processing (online) of large raster data from the Geo Data Portal of the U.S. Geological Survey.
wopr
(R-wopr?): provides API access to the WorldPop Open Population Repository database.
rdhs
(rdhs2019?) : Demographic and Health Survey (DHS) client API and data managements.
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:
- khm_admbnda_adm0_gov_20181004.zip
- khm_admbnda_adm1_gov_20181004.zip
- khm_admbnda_adm2_gov_20181004.zip
- khm_admbnda_adm3_gov_20181004.zip
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.
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)
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)
<- st_read("data_cambodia/cambodia.gpkg", layer = "district", quiet = TRUE)
district <- st_transform(district, 3857)
district <- get_tiles(x = district, zoom = 10, crop = TRUE)
osm_tiles 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)
<- st_read("data_cambodia/cambodia.gpkg", layer = "country", quiet = TRUE)
country <- opq(bbox = st_bbox(st_transform(country, 4326))) #Define the bounding box
ext <- add_osm_feature(opq = ext, key = 'amenity', value = "hospital") #Health Center Extraction
query <- osmdata_sf(query)
hospital <- unique_osmdata(hospital) #Result reduction (points composing polygon are detected) hospital
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$osm_points
hospital_point <- hospital$osm_polygons #Extracting centroids of polygons
hospital_poly <- st_centroid(hospital_poly)
hospital_poly_centroid
<- intersect(names(hospital_point), names(hospital_poly_centroid)) #Identify fields in Cambodia boundary
cambodia_point <- rbind(hospital_point[, cambodia_point], hospital_poly_centroid[, cambodia_point]) #Gather the 2 objects hospitals
Result display
library(mapview)
mapview(country) + mapview(hospitals)
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.
For administrative boundaries, check here the administrative levels by country:
library(osmextract)
library(mapsf)
<- oe_get(
province 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)
<- oe_get(
roads 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)
<- st_as_sf(read.csv("data_cambodia/adress.csv"), coords = c("long", "lat"), crs = 4326)
place_sf 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)
<- st_as_sf(data.frame(x = 0.5, y = 45.5), coords = c("x", "y"), crs = 4326)
test_point 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)
<- get_tiles(x = place_sf, zoom = 12)
osm 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)
<- data.frame(
test_adresses address = c("Phnom Penh International Airport, Phnom Penh, Cambodia",
"Khmer Soviet Friendship Hospital, Phnom Penh, Cambodia"))
<- geocode(test_adresses, address)
places1 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.
2.4.2 banR
(Base Adresse Nationale)
# remotes::install_github("joelgombin/banR")
library(banR)
<- data.frame(
mes_adresses address = c("19 rue Michel Bakounine, 29600 Morlaix, France",
"2 Allee Emile Pouget, 920128 Boulogne-Billancourt")
)<- geocode_tbl(tbl = mes_adresses, adresse = address)
places2 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
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.
mapedit
website