Spatial interaction approach about upthat HOT 14 OPEN

# packages ----------------------------------------------------------------

remotes::install_github("zonebuilders/zonebuilder")
#> Skipping install of 'zonebuilder' from a github remote, the SHA1 (239cede3) has not changed since last install.
#>   Use `force = TRUE` to force installation
remotes::install_github("ITSLeeds/od")
#> Skipping install of 'od' from a github remote, the SHA1 (6cb500ea) has not changed since last install.
#>   Use `force = TRUE` to force installation

library(zonebuilder)
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright


# basic parameters --------------------------------------------------------

# idea: fundamental attributes per city

city_name = "accra"
has_od_data = FALSE
has_zone_data = FALSE


# get city data -----------------------------------------------------------

ne_cities = rnaturalearth::ne_download(type = "populated_places", returnclass = "sf", scale = 50)
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/RtmpMBkTHs", layer: "ne_50m_populated_places"
#> with 1249 features
#> It has 119 fields
#> Integer64 fields read as strings:  wof_id ne_id
nrow(ne_cities)
#> [1] 1249
sum(ne_cities$POP_MAX) / 7e9 # around 20% of world's population
#> [1] 0.2108177

# for smaller cities
# ne_cities_10 = rnaturalearth::ne_download(type = "populated_places", returnclass = "sf", scale = 10)
# nrow(ne_cities_10)
# sum(ne_cities_10$POP_MAX, na.rm = TRUE) / 7e9 # around 33% of world's population

cities_match = grepl(pattern = city_name, x = ne_cities$NAME, ignore.case = TRUE)
n_city_name_in_city_data = sum(cities_match)

if(n_city_name_in_city_data == 1) {
  message("The city is in the city database" )
  city = ne_cities[cities_match, ]
}
#> The city is in the city database

region = osmdata::getbb(place_name = city_name, format_out = "sf_polygon")
#> No polygonal boundary for accra
if(is.null(region)) {
  message("Could not find boundary data from OSM. Getting data from another source.")
  if(city_name == "accra") {
    region = sf::read_sf("~/atfutures/who-data/accra/AccraMetroAdmin/AccraMetroAdmin2.shp")
  }
  # mapview::mapview(region)
}
#> Could not find boundary data from OSM. Getting data from another source.

city_centre = tmaptools::geocode_OSM(q = city_name, as.sf = TRUE)
mapview::mapview(region) + mapview::mapview(city_centre)

if(!has_zone_data) {
  city_zones = zonebuilder::zb_zone(x = city_centre, area = region)["label"]
  city_zone_centroids = sf::st_centroid(city_zones)
  zonebuilder::zb_plot(city_zones)
  plot(city_zone_centroids, add = TRUE, col = "grey")
}
#> although coordinates are longitude/latitude, st_contains assumes that they are planar
#> Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
#> although coordinates are longitude/latitude, st_within assumes that they are planar
#> Warning in st_centroid.sf(city_zones): st_centroid assumes attributes are
#> constant over geometries of x
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
#> of_largest_polygon): st_centroid does not give correct centroids for longitude/
#> latitude data
#> Error in `[[<-.data.frame`(`*tmp*`, i, value = character(0)): replacement has 0 rows, data has 40


# generate od data --------------------------------------------------------

od_region = od::points_to_od(city_zone_centroids, interzone_only = TRUE)
nrow(od_region) # 1560
#> [1] 1560
od_region = od::od_oneway(od_region)
nrow(od_region) # 780
#> [1] 780
head(od_region)
#>     o   d
#> 1   A B01
#> 2   A B02
#> 3 B01 B02
#> 4   A B03
#> 5 B01 B03
#> 6 B02 B03
od_sf = od::od_to_sf(x = od_region, z = city_zone_centroids)
plot(od_sf)