2  Data Acquisition

2.1 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)
Spherical geometry (s2) switched off

To crate 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.2 Online databases

2.3 OpenStreetMap

OpenStreetMap (OSM) is a participatory mapping project that aim s to buil 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.3.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.3.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.3.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.3.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.3.3 Import OSM data

2.3.3.1 osmdata

The package osmdata (Padgham et al. 2017) 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 fast food-food place. To obtain a coherent point layer we can use the centroids of the polygons.

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.3.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.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 × 17
  address        latit…¹ longi…² resul…³ resul…⁴ resul…⁵ resul…⁶ resul…⁷ resul…⁸
  <chr>            <dbl>   <dbl> <chr>     <dbl> <chr>   <chr>   <chr>   <chr>  
1 19 rue Michel…    48.6   -3.82 19 Rue…    0.81 housen… 29151_… 19      Rue Mi…
2 2 Allee Emile…    48.8    2.24 2 Allé…    0.83 housen… 92012_… 2       Allée …
# … with 8 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>, 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.