3 Using vector data

3.1 Importing and exporting data

The st_read() and st_write() function are used to import and export many types of files. The following lines import the administrative data in district level layer located in the district.shp shapefile file.

library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.1; sf_use_s2() is TRUE

district <- st_read("data_cambodia/district.shp")

Reading layer `district' from data source 
  `/home/lucas/Documents/ForgeIRD/rspatial-for-onehealth/data_cambodia/district.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 197 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 211534.7 ymin: 1149105 xmax: 784612.1 ymax: 1625495
Projected CRS: WGS 84 / UTM zone 48N

Shapefile format limitations

For the multiple limitations of this format (multi-file, limited number of records…) we advise you to prefer another format such as the geopackage *.gpkg. All the good reasons not to use the shapefile are here.

A geopackage is a database, to load a layer, you must know its name

st_layers("data_cambodia/cambodia.gpkg")

Driver: GPKG 
Available layers:
  layer_name     geometry_type features fields              crs_name
1    country     Multi Polygon        1     10 WGS 84 / UTM zone 48N
2   district     Multi Polygon      197     10 WGS 84 / UTM zone 48N
3  education     Multi Polygon       25     19 WGS 84 / UTM zone 48N
4   hospital             Point      956     13 WGS 84 / UTM zone 48N
5      cases       Multi Point      972      2 WGS 84 / UTM zone 48N
6       road Multi Line String        6      9 WGS 84 / UTM zone 48N

road <- st_read("data_cambodia/cambodia.gpkg", layer = "road")

Reading layer `road' from data source 
  `/home/lucas/Documents/ForgeIRD/rspatial-for-onehealth/data_cambodia/cambodia.gpkg' 
  using driver `GPKG'
Simple feature collection with 6 features and 9 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 212377 ymin: 1152214 xmax: 784654.7 ymax: 1625281
Projected CRS: WGS 84 / UTM zone 48N

library(sf)

district = st_read("data_cambodia/cambodia.gpkg", layer = "district")   #import district data

Reading layer `district' from data source 
  `/home/lucas/Documents/ForgeIRD/rspatial-for-onehealth/data_cambodia/cambodia.gpkg' 
  using driver `GPKG'
Simple feature collection with 197 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 211534.7 ymin: 1149105 xmax: 784612.1 ymax: 1625495
Projected CRS: WGS 84 / UTM zone 48N

The following lines export the district object to a data folder in geopackage and shapefile format.

st_write(obj = district, dsn = "data_cambodia/district.gpkg", delete_layer = TRUE)
st_write(obj = district, "data_cambodia/district.shp", layer_options = "ENCODING=UTF-8", delete_layer = TRUE)

3.2 Display

Preview of the variables via the function head() and plot().

head(district)

Simple feature collection with 6 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 300266.9 ymin: 1180566 xmax: 767313.9 ymax: 1563861
Projected CRS: WGS 84 / UTM zone 48N
       ADM2_EN ADM2_PCODE      ADM1_EN ADM1_PCODE  Male Female T_POP Area.Km2.
1    Aek Phnum     KH0205   Battambang       KH02 41500  43916 85416 1067.8638
2 Andoung Meas     KH1601 Ratanak Kiri       KH16  7336   7372 14708  837.7064
3   Angk Snuol     KH0808       Kandal       KH08 45436  47141 92577  183.9050
4 Angkor Borei     KH2101        Takeo       KH21 26306  27168 53474  301.0502
5  Angkor Chey     KH0701       Kampot       KH07 42448  44865 87313  316.7576
6  Angkor Chum     KH1701     Siemreap       KH17 34269  34576 68845  478.6988
    Status      DENs                           geom
1 <4500km2  79.98773 MULTIPOLYGON (((306568.1 14...
2 <4500km2  17.55747 MULTIPOLYGON (((751459.2 15...
3 <4500km2 503.39580 MULTIPOLYGON (((471954.3 12...
4 <4500km2 177.62485 MULTIPOLYGON (((490048.2 12...
5 <4500km2 275.64610 MULTIPOLYGON (((462702.2 12...
6 <4500km2 143.81696 MULTIPOLYGON (((363642.5 15...

plot(district)

for Geometry display only.

plot(st_geometry(district))

3.3 Coordinate systems

3.3.1 Look up the coordinate system of an object

The function st_crs() makes it possible to consult the system of coordinates used and object sf.

st_crs(district)

Coordinate Reference System:
  User input: WGS 84 / UTM zone 48N 
  wkt:
PROJCRS["WGS 84 / UTM zone 48N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 48N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",105,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 102°E and 108°E, northern hemisphere between equator and 84°N, onshore and offshore. Cambodia. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Russian Federation. Singapore. Thailand. Vietnam."],
        BBOX[0,102,84,108]],
    ID["EPSG",32648]]

3.3.2 Changing the coordinate system of an object

The function st_transform() allows to change the coordinate system of an sf object, to re-project it.

plot(st_geometry(district))
title("WGS 84 / UTM zone 48N")

dist_reproj <- st_transform(district, "epsg:4326")
plot(st_geometry(dist_reproj))
title("WGS84")

The Spatial Reference site provides reference for a large number of coordinate systems.

3.4 Selection by attributes

The object sf are data.frame, so you can select their rows and columns in the same way as data.frame.

# row Selection
district[1:2, ]

Simple feature collection with 2 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 300266.9 ymin: 1449408 xmax: 767313.9 ymax: 1563861
Projected CRS: WGS 84 / UTM zone 48N
       ADM2_EN ADM2_PCODE      ADM1_EN ADM1_PCODE  Male Female T_POP Area.Km2.
1    Aek Phnum     KH0205   Battambang       KH02 41500  43916 85416 1067.8638
2 Andoung Meas     KH1601 Ratanak Kiri       KH16  7336   7372 14708  837.7064
    Status     DENs                           geom
1 <4500km2 79.98773 MULTIPOLYGON (((306568.1 14...
2 <4500km2 17.55747 MULTIPOLYGON (((751459.2 15...

district[district$ADM1_EN == "Phnom Penh", ]

Simple feature collection with 12 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 468677.5 ymin: 1262590 xmax: 505351.9 ymax: 1297419
Projected CRS: WGS 84 / UTM zone 48N
First 10 features:
              ADM2_EN ADM2_PCODE    ADM1_EN ADM1_PCODE  Male Female  T_POP
29        Chamkar Mon     KH1201 Phnom Penh       KH12 52278  54478 106756
31        Chbar Ampov     KH1212 Phnom Penh       KH12 64816  68243 133059
43    Chraoy Chongvar     KH1210 Phnom Penh       KH12 30920  31087  62007
48            Dangkao     KH1205 Phnom Penh       KH12 46999  48525  95524
50          Doun Penh     KH1202 Phnom Penh       KH12 33844  36471  70315
93          Mean Chey     KH1206 Phnom Penh       KH12 68381  70366 138747
117        Praek Pnov     KH1211 Phnom Penh       KH12 27566  27698  55264
118 Prampir Meakkakra     KH1203 Phnom Penh       KH12 31091  33687  64778
133       Pur SenChey     KH1209 Phnom Penh       KH12 95050 109297 204347
141        Russey Keo     KH1207 Phnom Penh       KH12 67357  68419 135776
     Area.Km2.   Status       DENs                           geom
29   11.049600 <4500km2  9661.5265 MULTIPOLYGON (((494709.4 12...
31   86.780498 <4500km2  1533.2823 MULTIPOLYGON (((498855.3 12...
43   85.609156 <4500km2   724.3034 MULTIPOLYGON (((491161.3 12...
48  113.774833 <4500km2   839.5881 MULTIPOLYGON (((489191.1 12...
50    7.734808 <4500km2  9090.7234 MULTIPOLYGON (((492447.1 12...
93   28.998026 <4500km2  4784.7051 MULTIPOLYGON (((491068.2 12...
117 115.384300 <4500km2   478.9560 MULTIPOLYGON (((481483.3 12...
118   2.224892 <4500km2 29115.1253 MULTIPOLYGON (((491067.6 12...
133 148.357984 <4500km2  1377.3913 MULTIPOLYGON (((479078.8 12...
141  23.381517 <4500km2  5806.9800 MULTIPOLYGON (((490264.8 12...

# column selection
district[district$ADM1_EN == "Phnom Penh", 1:4]

Simple feature collection with 12 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 468677.5 ymin: 1262590 xmax: 505351.9 ymax: 1297419
Projected CRS: WGS 84 / UTM zone 48N
First 10 features:
              ADM2_EN ADM2_PCODE    ADM1_EN ADM1_PCODE
29        Chamkar Mon     KH1201 Phnom Penh       KH12
31        Chbar Ampov     KH1212 Phnom Penh       KH12
43    Chraoy Chongvar     KH1210 Phnom Penh       KH12
48            Dangkao     KH1205 Phnom Penh       KH12
50          Doun Penh     KH1202 Phnom Penh       KH12
93          Mean Chey     KH1206 Phnom Penh       KH12
117        Praek Pnov     KH1211 Phnom Penh       KH12
118 Prampir Meakkakra     KH1203 Phnom Penh       KH12
133       Pur SenChey     KH1209 Phnom Penh       KH12
141        Russey Keo     KH1207 Phnom Penh       KH12
                              geom
29  MULTIPOLYGON (((494709.4 12...
31  MULTIPOLYGON (((498855.3 12...
43  MULTIPOLYGON (((491161.3 12...
48  MULTIPOLYGON (((489191.1 12...
50  MULTIPOLYGON (((492447.1 12...
93  MULTIPOLYGON (((491068.2 12...
117 MULTIPOLYGON (((481483.3 12...
118 MULTIPOLYGON (((491067.6 12...
133 MULTIPOLYGON (((479078.8 12...
141 MULTIPOLYGON (((490264.8 12...

3.5 Spatial selection

3.5.1 Intersections

Selection of roads that are intersecting dangkao district

road <- st_read("data_cambodia/cambodia.gpkg", layer = "road", quiet = TRUE) %>% st_cast("LINESTRING")
dangkao <-  district[district$ADM2_EN == "Dangkao", ]
inter <- st_intersects(x = road, y = dangkao, sparse = FALSE)
head(inter)

      [,1]
[1,] FALSE
[2,] FALSE
[3,] FALSE
[4,] FALSE
[5,] FALSE
[6,] FALSE

dim(inter)

[1] 108285      1

The inter object is a matrix which indicates for each of element of the road object (6 elements) whether it intersects each elements the dangkao object (1 element). The dimension of the matrix is therefore indeed 6 rows * 1 column. Note the use of the parameter sparse = FALSE here. It is then possible to create a column from this object:

road$intersect_dangkao <- inter
plot(st_geometry(dangkao), col = "lightblue")
plot(st_geometry(road), add = TRUE)
plot(st_geometry(road[road$intersect_dangkao, ]),
      col = "tomato", lwd = 1.5, add = TRUE)

3.5.1.1 Difference between `sparse = TRUE` and `sparse = FALSE`

sparse = TRUE

inter <- st_intersects(x = grid, y = pt, sparse = TRUE)
inter

Sparse geometry binary predicate list of length 4, where the predicate
was `intersects'
 1: (empty)
 2: 6, 7
 3: 1, 4
 4: 2, 3, 5, 8

sparse = FALSE

inter <- st_intersects(x = grid, y = pt, sparse = FALSE)
rownames(inter) <- grid$id
colnames(inter) <- pt$id
inter

      a     b     c     d     e     f     g     h
1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
2 FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE
3  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
4 FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE

3.5.2 Contains / Within

Selection of roads contained in the municipality of Dangkao. The function st_within() works like the function st_intersects()

road$within_dangkao <- st_within(road, dangkao, sparse = FALSE)
plot(st_geometry(dangkao), col = "lightblue")
plot(st_geometry(road), add = TRUE)
plot(st_geometry(road[road$within_dangkao, ]), col = "tomato",
     lwd = 2, add = TRUE)

3.6 Operation of geometries

3.6.1 Extract centroids

dist_c <- st_centroid(district)
plot(st_geometry(district))
plot(st_geometry(dist_c), add = TRUE, cex = 1.2, col = "red", pch = 20)

3.6.2 Aggregate polygons

cambodia_dist <- st_union(district)                 
plot(st_geometry(district), col = "lightblue")
plot(st_geometry(cambodia_dist), add = TRUE, lwd = 2, border = "red")

3.6.3 Aggregate polygons based on a variable

dist_union  <- aggregate(x = district[,c("T_POP")],
                   by = list(STATUT = district$Status),
                   FUN = "sum")
plot(dist_union)

3.6.4 Create a buffer zone

dangkao_buffer <- st_buffer(x = dangkao, dist = 1000)
plot(st_geometry(dangkao_buffer), col = "#E8DAEF", lwd=2, border = "#6C3483")
plot(st_geometry(dangkao), add = TRUE, lwd = 2)

3.6.5 Making an intersection

By using the function st_intersection() we will cut one layer by another.

library(magrittr)
# creation of a buffer zone around the centroid of the municipality of Dangkao district
# using the pipe
zone <- st_geometry(dangkao) %>%
  st_centroid() %>%
  st_buffer(30000)
plot(st_geometry(district))
plot(zone, border = "#F06292", lwd = 2, add = TRUE)

dist_z <- st_intersection(x = district, y = zone)
plot(st_geometry(district))
plot(st_geometry(dist_z), col="#AF7AC5", border="#F9E79F", add=T)

plot(st_geometry(dist_z))

3.6.6 Create regular grid

The function st_make_grid() allows you to create regular grid. The function produce and object sfc, you must then use the function st_sf() to transform the object sfc into and object sf. During this transformation we add here a column of unique identifiers.

grid <- st_make_grid(x = district, cellsize = 10000)
grid <- st_sf(ID = 1:length(grid), geom = grid)

plot(st_geometry(grid), col = "grey", border = "white")
plot(st_geometry(district), border = "grey50", add = TRUE)

3.6.7 Counting points in a polygon (in a grid tile)

# selection of grid tiles that intersect the district

inter <- st_intersects(grid, cambodia_dist, sparse = FALSE)
grid <- grid[inter, ]

case_cambodia <- st_read("data_cambodia/cambodia.gpkg", layer = "cases" , quiet = TRUE)
plot(st_geometry(grid), col = "grey", border = "white")
plot(st_geometry(case_cambodia), pch = 20, col = "red", add = TRUE, cex = 0.8)

inter <- st_intersects(grid, case_cambodia, sparse = TRUE)
length(inter)

[1] 1964

Here we use the argument sparse = TRUE. The inter object is a list the length of the grid and each item in the list contain the index of the object items of cases and grid intersection.

For example grid tile 35th intersect with four cases 97, 138, 189, 522, 624, 696

inter[35]

[[1]]
[1]  97 138 189 522 624 696

plot(st_geometry(grid[35, ]))
plot(st_geometry(case_cambodia), add = T)
plot(st_geometry(case_cambodia[c(97, 138, 189, 522, 624, 696), ]), 
     col = "red", pch = 19, add = TRUE)

To count number of case, simply go to the list and report length of the elements.

grid$nb_case <- sapply(X = inter, FUN = length)   # create 'nb_case' column to store number of health centers in each grid tile 
plot(grid["nb_case"])

3.6.8 Aggregate point values into polygons

In this example we import a csv file that contain data from a population grid. Once import we transform it data.frame into an object sf.

The objective is to aggregate the values id these points (the population contained in the “DENs” field) in the municipalities of the district.

pp_pop_raw <- read.csv("data_cambodia/pp_pop_dens.csv")            # import file
pp_pop_raw$id <- 1:nrow(pp_pop_raw)                                # adding a unique identifier
pp_pop <- st_as_sf(pp_pop_raw, coords = c("X", "Y"), crs = 32648)  # Transform into object sf
pp_pop <- st_transform(pp_pop, st_crs(district))                   # Transform projection
inter <- st_intersection(pp_pop, district)                         # Intersection
inter

Simple feature collection with 1295 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 469177.5 ymin: 1263090 xmax: 505177.5 ymax: 1297090
Projected CRS: WGS 84 / UTM zone 48N
First 10 features:
        DENs  id    ADM2_EN ADM2_PCODE ADM1_EN ADM1_PCODE  Male Female T_POP
149       NA 149 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
150       NA 150 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
151       NA 151 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
186       NA 186 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
187       NA 187 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
188       NA 188 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
223       NA 223 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
224       NA 224 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
225       NA 225 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
226 3.400075 226 Angk Snuol     KH0808  Kandal       KH08 45436  47141 92577
    Area.Km2.   Status   DENs.1                 geometry
149   183.905 <4500km2 503.3958 POINT (469177.5 1267090)
150   183.905 <4500km2 503.3958 POINT (470177.5 1267090)
151   183.905 <4500km2 503.3958 POINT (471177.5 1267090)
186   183.905 <4500km2 503.3958 POINT (469177.5 1268090)
187   183.905 <4500km2 503.3958 POINT (470177.5 1268090)
188   183.905 <4500km2 503.3958 POINT (471177.5 1268090)
223   183.905 <4500km2 503.3958 POINT (469177.5 1269090)
224   183.905 <4500km2 503.3958 POINT (470177.5 1269090)
225   183.905 <4500km2 503.3958 POINT (471177.5 1269090)
226   183.905 <4500km2 503.3958 POINT (472177.5 1269090)

By using the function st_intersection() we add to each point of the grid all the information on the municipality in which it is located.

We can then use the function aggregate() to aggregate the population by municipalities.

resultat <- aggregate(x = list(pop_from_grid = inter$DENs), 
                      by = list(ADM2_EN = inter$ADM2_EN), 
                      FUN = "sum")
head(resultat)

          ADM2_EN pop_from_grid
1      Angk Snuol            NA
2     Chamkar Mon    10492.7159
3     Chbar Ampov     1593.9593
4 Chraoy Chongvar     1434.1785
5         Dangkao      942.3595
6       Doun Penh    10781.8026

We can then create a new object with this result.

dist_result <- merge(district, resultat, by = "ADM2_EN", all.x = TRUE)
dist_result

Simple feature collection with 197 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 211534.7 ymin: 1149105 xmax: 784612.1 ymax: 1625495
Projected CRS: WGS 84 / UTM zone 48N
First 10 features:
         ADM2_EN ADM2_PCODE        ADM1_EN ADM1_PCODE  Male Female T_POP
1      Aek Phnum     KH0205     Battambang       KH02 41500  43916 85416
2   Andoung Meas     KH1601   Ratanak Kiri       KH16  7336   7372 14708
3     Angk Snuol     KH0808         Kandal       KH08 45436  47141 92577
4   Angkor Borei     KH2101          Takeo       KH21 26306  27168 53474
5    Angkor Chey     KH0701         Kampot       KH07 42448  44865 87313
6    Angkor Chum     KH1701       Siemreap       KH17 34269  34576 68845
7    Angkor Thum     KH1702       Siemreap       KH17 13802  14392 28194
8  Anlong Veaeng     KH2201 Oddar Meanchey       KH22 24122  23288 47410
9          Aoral     KH0504   Kampong Speu       KH05 19874  19956 39830
10      Ba Phnum     KH1401      Prey Veng       KH14 46562  49852 96414
   Area.Km2.   Status      DENs pop_from_grid                       geometry
1  1067.8638 <4500km2  79.98773            NA MULTIPOLYGON (((306568.1 14...
2   837.7064 <4500km2  17.55747            NA MULTIPOLYGON (((751459.2 15...
3   183.9050 <4500km2 503.39580            NA MULTIPOLYGON (((471954.3 12...
4   301.0502 <4500km2 177.62485            NA MULTIPOLYGON (((490048.2 12...
5   316.7576 <4500km2 275.64610            NA MULTIPOLYGON (((462702.2 12...
6   478.6988 <4500km2 143.81696            NA MULTIPOLYGON (((363642.5 15...
7   357.8890 <4500km2  78.77862            NA MULTIPOLYGON (((376584.4 15...
8  1533.5702 <4500km2  30.91479            NA MULTIPOLYGON (((404936.4 15...
9  2381.7084 <4500km2  16.72329            NA MULTIPOLYGON (((414000.6 13...
10  342.3439 <4500km2 281.62910            NA MULTIPOLYGON (((545045.4 12...

3.7 Measurements

3.7.1 Create a distance matrix

If the dataset’s projection system is specified, the distance are expressed in the projection measurement unit (most often in meter)

mat <- st_distance(x = dist_c, y = dist_c)
mat[1:5,1:5]

Units: [m]
         [,1]     [,2]      [,3]      [,4]      [,5]
[1,]      0.0 425993.7 232592.12 298254.12 299106.92
[2,] 425993.7      0.0 386367.88 414428.82 452431.87
[3,] 232592.1 386367.9      0.00  67060.05  82853.88
[4,] 298254.1 414428.8  67060.05      0.00  40553.15
[5,] 299106.9 452431.9  82853.88  40553.15      0.00

3.7.2 Calculate routes

The package osrm (R-osrm?) acts as an interface R and the OSRM (luxen-vetter-2011?). This package allows to calculate time and distance matrices, road routes, isochrones. The package uses the OSRM demo server by default. In case of intensive use it is strongly recommended to use your own instance of OSRM (with Docker).

3.7.2.1 Calculate a route

The fonction osrmRoute() allows you to calculate routes.

library(sf)
library(osrm)
library(maptiles)
district <- st_read("data_cambodia/cambodia.gpkg",layer = "district", quiet = TRUE)
district <- st_transform(district, 32648)

odongk  <- district[district$ADM2_PCODE == "KH0505", ] # Itinerary between Odongk district and Toul Kouk
takmau <- district[district$ADM2_PCODE == "KH0811",]
route <- osrmRoute(src = odongk, 
                   dst = takmau, 
                   returnclass = "sf")
osm <- get_tiles(route, crop = TRUE)
plot_tiles(osm)
plot(st_geometry(route), col = "#b23a5f", lwd = 6, add = T)
plot(st_geometry(route), col = "#eee0e5", lwd = 1, add = T)

3.7.2.2 Calculation of a time matrix

The function osrmTable() makes it possible to calculate matrices of distances or times by road.

In this example we calculate a time matrix between 2 addresses and health centers in Phnom Penh on foot.

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

hospital_pp <- hospital[hospital$PCODE == "12", ]     # Selection of health centers in Phnom Penh

adresses <- data.frame(adr = c("Royal Palace Park, Phnom Penh Phnom, Cambodia",
                              "Wat Phnom Daun Penh, Phnom Penh, Cambodia"))  # Geocoding of 2 addresses in Phnom Penh

places <- tidygeocoder::geocode(.tbl = adresses,address = adr)
places

# A tibble: 2 × 3
  adr                                             lat  long
  <chr>                                         <dbl> <dbl>
1 Royal Palace Park, Phnom Penh Phnom, Cambodia  11.6  105.
2 Wat Phnom Daun Penh, Phnom Penh, Cambodia      11.6  105.

# Calculation of the distance matrix between the 2 addresses and the health center in Phnom Penh

cal_mat <- osrmTable(src = places[,c(1,3,2)], 
                 dst = hospital_pp, 
                 osrm.profile = "foot")

cal_mat$durations[1:2, 1:5]

                                               684  685  686  687  691
Royal Palace Park, Phnom Penh Phnom, Cambodia 55.9 71.6 64.4 40.2 76.7
Wat Phnom Daun Penh, Phnom Penh, Cambodia     60.1 80.4 40.1 32.8 53.1

# Which address has better accessibility to health center in Phnom Penh?

boxplot(t(cal_mat$durations[,]), cex.axis = 0.7)

3.1 Importing and exporting data

3.2 Display

3.3 Coordinate systems

3.3.1 Look up the coordinate system of an object

3.3.2 Changing the coordinate system of an object

3.4 Selection by attributes

3.5 Spatial selection

3.5.1 Intersections

3.5.1.1 Difference between sparse = TRUE and sparse = FALSE

3.5.2 Contains / Within

3.6 Operation of geometries

3.6.1 Extract centroids

3.6.2 Aggregate polygons

3.6.3 Aggregate polygons based on a variable

3.6.4 Create a buffer zone

3.6.5 Making an intersection

3.6.6 Create regular grid

3.6.7 Counting points in a polygon (in a grid tile)

3.6.8 Aggregate point values into polygons

3.7 Measurements

3.7.1 Create a distance matrix

3.7.2 Calculate routes

3.7.2.1 Calculate a route

3.7.2.2 Calculation of a time matrix

3.5.1.1 Difference between `sparse = TRUE` and `sparse = FALSE`