Download Ipixuna OpenStreetMap

Load required libraries and data

rm(list = ls())
library(day2day)
library(sf)

#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(purrr)
library(ggplot2)
library(osmplotr)

#> Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright

path_main <- git_path()
path_data <- file.path(path_main, "data")
path_raw <- file.path(path_data, "raw")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")
path_fig <- file.path(path_data, "figures")
path_src <- file.path(path_main, "src")

fidata <- file.path(path_processed, "fi-items-ipixuna-urban.gpkg") |>
    st_read(as_tibble = TRUE)

#> Reading layer `fi-items-ipixuna-urban' from data source 
#>   `/home/rstudio/documents/projects/food-insecurity-mapping/data/processed/fi-items-ipixuna-urban.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 200 features and 36 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -71.70038 ymin: -7.06058 xmax: -71.68109 ymax: -7.03724
#> Geodetic CRS:  WGS 84

Create bbox for OpenStreetMap

# get bounding box
ipi_bbox <- st_bbox(fidata) |>
    st_as_sfc() |>
    st_buffer(500) |>
    st_bbox() |>
    st_as_sfc()

Download OpenStreetMap features

# bbox
bbox <- get_bbox(st_bbox(ipi_bbox))
bbox_3857 <- get_bbox(st_bbox(st_transform(ipi_bbox, st_crs(3857))))

# download some objects
ipi_land <- extract_osm_objects(bbox, key = "landuse", return_type = "polygon") |>
  st_transform(st_crs(3857))

#> Issuing query to Overpass API ...

#> Announced endpoint: none

#> Query complete!

#> converting OSM data to sf format

ipi_highway <- extract_osm_objects(bbox, key = "highway", return_type = "line") |>
  st_transform(st_crs(3857))

#> Issuing query to Overpass API ...

#> Announced endpoint: none

#> Query complete!

#> converting OSM data to sf format

ipi_water_multipoly <- extract_osm_objects(bbox, key = "water", return_type = "multipolygon") |>
  st_transform(st_crs(3857))

#> Issuing query to Overpass API ...

#> Announced endpoint: none

#> Query complete!

#> converting OSM data to sf format

ipi_water_poly <- extract_osm_objects(bbox, key = "water", return_type = "polygon") |>
  st_transform(st_crs(3857))

#> Issuing query to Overpass API ...

#> Announced endpoint: none

#> Query complete!

#> converting OSM data to sf format

ipi_water_line <- extract_osm_objects(bbox, key = "waterway", return_type = "line") |>
  st_transform(st_crs(3857))

#> Issuing query to Overpass API ...

#> Announced endpoint: none

#> Query complete!

#> converting OSM data to sf format

ipi_aeroway <- extract_osm_objects(bbox, key = "aeroway", return_type = "polygon") %>%
  st_transform(st_crs(3857))

#> Issuing query to Overpass API ...

#> Announced endpoint: none

#> Query complete!

#> converting OSM data to sf format

# simple visualization
ggplot(st_transform(ipi_bbox, st_crs(3857))) +
    geom_sf(fill = NA) +
    geom_sf(data = ipi_land, fill = 3) +
    geom_sf(data = st_union(st_buffer(ipi_highway, 5))) +
    geom_sf(data = st_buffer(st_union(ipi_water_line), 15), fill = "lightblue") +
    geom_sf(data = st_union(subset(ipi_water_multipoly, water == "river")), fill = "lightblue") +
    geom_sf(data = subset(ipi_water_multipoly, water != "river"), fill = "lightblue") +
    geom_sf(data = ipi_water_poly, fill = "lightblue") +
    geom_sf(data = ipi_aeroway, fill = "gray50") +
    coord_sf(xlim = bbox_3857[1, ], ylim = bbox_3857[2, ]) +
    theme_bw()

Save Ipixuna OpenStreetMap data

list(box = ipi_bbox, land = ipi_land, highway = ipi_highway,
    water_multipoly = ipi_water_multipoly, water_poly = ipi_water_poly,
    water_line = ipi_water_line, aeroway = ipi_aeroway) |>
    saveRDS(file.path(path_processed, "ipixuna-osm.rds"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>   4.389   1.085  11.228