Skip to content

Instantly share code, notes, and snippets.

@datagistips
Last active January 13, 2023 20:08
Show Gist options
  • Save datagistips/6cd56402b7fadc1f858d1cdfa717b3e7 to your computer and use it in GitHub Desktop.
Save datagistips/6cd56402b7fadc1f858d1cdfa717b3e7 to your computer and use it in GitHub Desktop.
Compare city sizes by putting a city onto another
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment