Compare city sizes by putting a city onto another

CompareCities.R

library(tidyverse)
library(jsonlite)
library(glue)
library(sf)
library(leaflet)

make_url <- function(cityname) {
  cityname <- URLencode(cityname)
  glue("https://nominatim.openstreetmap.org/search?city={cityname}&format=geojson&polygon_geojson=1")
}

get_feature <- function(feature) {
  
  # Display Name and coordinates
  display_name <- feature$properties$display_name
  message(display_name)
  coords <- feature$geometry$coordinates
  message(feature$geometry$type)
  
  # If Point
  if(feature$geometry$type == "Point") {
    long <- coords[[1]]
    lat <- coords[[2]]
    g <- st_sfc(st_point(cbind(long, lat)))
    sdf <-  st_sf(display_name = display_name, g) %>% st_set_crs(4326)
  }
  
  # If Polygon
  if(feature$geometry$type == "Polygon") {
    out <- vector(mode = "list")
    coords <- coords[[1]]
    for(i in 1:length(coords)) {
      elt <- coords[[i]]
      long <- elt[[1]] %>% as.numeric
      lat <- elt[[2]] %>% as.numeric
      out[[i]] <- cbind(long, lat)
    }
    res <- do.call(rbind, out)
    res <- rbind(res, head(res, 1))   
    g <- st_polygon(list(res)) %>% st_sfc
    sdf <- st_sf(display_name = display_name, g) %>% st_set_crs(4326)
  }
  
  # if MultiPolygon
  if(feature$geometry$type == "MultiPolygon") {
    out <- vector(mode = "list")
    for(i in 1:length(coords)) {
      elt <- coords[[i]]
      out2 <- vector(mode="list")
      for(j in 1:length(elt)) {
        elt2 <- elt[[j]]
        out3 <- vector(mode="list")
        for(k in 1:length(elt2)) {
          long <- elt2[[k]][[1]] %>% as.numeric
          lat <- elt2[[k]][[2]] %>% as.numeric
          out3[[k]] <- cbind(long, lat)
        }
        res2 <- do.call(rbind, out3)
        res2 <- rbind(res2, head(res2, 1))
        out2[[j]] <- res2
      }
      out[[i]] <- st_polygon(out2)
    }
    
    g <- st_sfc(out) %>% st_combine
    sdf <- st_sf(display_name = display_name, g) %>% st_set_crs(4326)
  }
  
  return(sdf)
}

search_city <- function(cityname) {
  # Request Nominatim
  url <- make_url(cityname)
  json <- read_json(url)
  features <- json$features
  
  if(length(features) == 0) return()
  
  out <- vector(mode = "list")
  for(i in 1:length(features)) {
    print(i)
    feature <- json$features[[i]]
    out[[i]] <- feature %>% get_feature
  }
  
  res <- do.call(rbind, out)
  
  # Deduplicate
  w <- which(!duplicated(res$display_name))
  res <- res[w, ]
  
  # Make valid
  res <- st_make_valid(res)
  
  return(res)
}

shift_feature <- function(feat_source, feat_target) {
  
  # If source = Point, return NULL
  if(st_geometry_type(feat_source) == "POINT") return()
  
  # Calculate difference of source coords and target coords
  # Get source coordinates
  coords_source <- st_centroid(feat_source) %>% st_coordinates
  
  # Coordinates of target
  if(st_geometry_type(feat_target) == "POINT") {
    coords_target <- feat_target %>% st_coordinates
  } else {
    coords_target <- st_centroid(feat_target) %>% st_coordinates
  }
  
  # Centroid will be used for shifting
  cntrd <- st_sfc(st_point(coords_target - coords_source))
  
  # Shift
  feat_shifted <- st_geometry(feat_source) + cntrd
  feat_shifted <- feat_shifted %>% st_set_crs(4326)
  
  return(feat_shifted)
}

shift_city <- function(city_source, city_target) {
  feat_source <- search_city(city_source)[1, ]
  feat_target <- search_city(city_target)[1, ]
  feat_shifted <- shift_feature(feat_source, feat_target)
  return(feat_shifted)
}

# Process ----
# Translate Paris to Houston
feat_shifted <- shift_city("Paris", "Houston")

# Visualize ----
# Get Berlin
feat_target <- search_city("Berlin") %>% head(1)

# Show Berlin and Paris shifted to Berlin
leaflet() %>% addTiles() %>% 
  addPolygons(data = feat_target) %>% 
  addPolygons(data = feat_shifted, color = "red")