Skip to content

Instantly share code, notes, and snippets.

@gedankenstuecke
Created March 7, 2024 19:52
Show Gist options
  • Save gedankenstuecke/d3e9575c76116ecababf5b858ee40c2b to your computer and use it in GitHub Desktop.
Save gedankenstuecke/d3e9575c76116ecababf5b858ee40c2b to your computer and use it in GitHub Desktop.
Making 3D population maps of Argentina
# ADAPTED FROM:
# - Spencer Schien's rayshader+kontur tutorial: https://justjensen.co/making-population-density-maps-with-rayrender-in-r/
# - Erik Jensen's adaptation of it: https://justjensen.co/making-population-density-maps-with-rayrender-in-r/
# INSTALL NECESSARY PACKAGES
install.packages("sf")
install.packages("rayrender")
install.packages("rayshader")
install.packages("R.utils")
install.packages('tidyverse')
install.packages("stars")
install.packages("RColorBrewer")
install.packages("colorspace")
install.packages("OpenStreetMap")
# LOAD PACKAGES
library(RColorBrewer)
library(colorspace)
library(sf)
library(R.utils)
library(tidyverse)
library(rayrender)
library(rayshader)
library(stars)
library(OpenStreetMap)
# LOAD KONTUR DATA: POPULATION & SHAPEFILES BY
# DOWNLOADING FROM https://data.humdata.org/dataset
# load population data for Argentina
df_pop_st <- st_read(gunzip("Downloads/kontur_population_AR_20231101.gpkg.gz", remove=FALSE, skip=TRUE))
# transform to WGS84
df_pop_transformed <- st_transform(df_pop_st,"WGS84")
#### ASPECT RATIO
# define aspect ratio based on bounding box
bb <- st_bbox(df_pop_transformed)
bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |>
st_sfc(crs = st_crs(df_pop_transformed))
bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |>
st_sfc(crs = st_crs(df_pop_transformed))
width <- st_distance(bottom_left, bottom_right)
top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |>
st_sfc(crs = st_crs(df_pop_transformed))
height <- st_distance(bottom_left, top_left)
# handle conditions of width or height being the longer side
if (width > height) {
w_ratio <- 1
h_ratio <- as.numeric(height / width)
} else {
h_ratio <- 1
w_ratio <- as.numeric(width / height)
}
# convert to raster so we can then convert to matrix
# set size to something larger after initial testing, but will increase render time/resources
size <- 1000
rast_arg <- st_rasterize(subset(df_pop_transformed,
df_pop_transformed$population > 10),
nx = floor(size * w_ratio),
ny = floor(size * h_ratio))
mat_arg <- matrix(rast_arg$population,
nrow = floor(size * w_ratio),
ncol = floor(size * h_ratio))
# set color palette
colors = brewer.pal(n=9, name = "PuRd")
texture <- grDevices::colorRampPalette(colors, bias = 3)(256)
swatchplot(texture)
rgl::close3d() # Close any existing renders
# make temp plot
mat_arg |>
height_shade(texture = texture) |>
plot_3d(heightmap = mat_arg,
zscale = 20,
solid = FALSE,
shadowdepth = 10,
windowsize = c(500,750),
background = 'grey80')
# set camera
render_camera(theta = -15, phi = 45, zoom = 0.5)
# render high-quality output, increase samples for better quality on same x/y resolution image
render_highquality(
filename = "argentina.png",
interactive = FALSE,
light = TRUE,
width = 2000,
height = 3000,
samples = 256
)
### RENDER ONLY PROVINCE, SANTA FE IN THIS CASE
# READ SHAPE FILE
shape <- read_sf("Downloads/arg_adm_unhcr2017_shp/arg_admbnda_adm1_unhcr2017.shp")
# SUBSET TO SANTA FE & TRANSFORM TO WGS84
shape_sf <- subset(shape,shape$ADM1_ES == "Santa Fe")
shape_sf <- st_transform(shape_sf,"WGS84")
df_pop_st <- st_transform(df_pop_st,"WGS84")
# INTERSECT TO LIMIT DATA TO SHAPE
stfe_df <- st_intersection(df_pop_st,shape_sf)
# make sure intersected data remains in wgs84
stfe_df_transformed <- st_transform(stfe_df,"WGS84")
#### ASPECT RATIO
# define aspect ratio based on bounding box
bb <- st_bbox(stfe_df_transformed)
bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |>
st_sfc(crs = st_crs(stfe_df_transformed))
bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |>
st_sfc(crs = st_crs(stfe_df_transformed))
width <- st_distance(bottom_left, bottom_right)
top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |>
st_sfc(crs = st_crs(stfe_df_transformed))
height <- st_distance(bottom_left, top_left)
# handle conditions of width or height being the longer side
if (width > height) {
w_ratio <- 1
h_ratio <- as.numeric(height / width)
} else {
h_ratio <- 1
w_ratio <- as.numeric(width / height)
}
# convert to raster so we can then convert to matrix
size <- 1000 # again, increase for real plot at end
rast_stfe <- st_rasterize(subset(stfe_df_transformed,
stfe_df_transformed$population > 10),
nx = floor(size * w_ratio),
ny = floor(size * h_ratio))
mat_stfe <- matrix(rast_stfe$population,
nrow = floor(size * w_ratio),
ncol = floor(size * h_ratio))
rgl::close3d() # Close exisitng windows
mat_stfe |>
height_shade(texture = texture) |>
plot_3d(heightmap = mat_stfe,
zscale = 25,
solid = FALSE,
shadowdepth = 10,
windowsize = c(400,600),
background = 'grey80')
render_camera(theta = -15, phi = 45, zoom = .5)
# render output
render_highquality(
filename = "santa_fe.png",
interactive = FALSE,
light = TRUE,
width = 2000,
height = 3000,
samples = 256
)
### MAP OF DEPARTMENT (GENERAL LOPEZ)
# load pop data again
df_pop_st <- st_read(gunzip("Downloads/kontur_population_AR_20231101.gpkg.gz", remove=FALSE, skip=TRUE))
# load shape
shape <- read_sf("Downloads/arg_adm_unhcr2017_shp/arg_admbnda_adm2_unhcr2017.shp")
shape_gl <- subset(shape,shape$ADM2_REF == "General Lopez")
shape_gl <- st_transform(shape_gl,"WGS84")
# combine
df_pop_st <- st_transform(df_pop_st,"WGS84")
gl_dt <- st_intersection(df_pop_st,shape_gl)
gl_dt_transformed <- st_transform(gl_dt,"WGS84")
#### ASPECT RATIO
# define aspect ratio based on bounding box
bb <- st_bbox(gl_dt_transformed)
bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |>
st_sfc(crs = st_crs(gl_dt_transformed))
bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |>
st_sfc(crs = st_crs(gl_dt_transformed))
width <- st_distance(bottom_left, bottom_right)
top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |>
st_sfc(crs = st_crs(gl_dt_transformed))
height <- st_distance(bottom_left, top_left)
if (width > height) {
w_ratio <- 1
h_ratio <- as.numeric(height / width)
} else {
h_ratio <- 1
w_ratio <- as.numeric(width / height)
}
# convert to raster 4 matrix
size <- 1000 # don't forget to set higher
rast_gl <- st_rasterize(subset(gl_dt_transformed,
gl_dt_transformed$population > 10),
nx = floor(size * w_ratio),
ny = floor(size * h_ratio))
mat_gl <- matrix(rast_gl$population,
nrow = floor(size * w_ratio),
ncol = floor(size * h_ratio))
rgl::close3d()
mat_gl |>
height_shade(texture = texture) |>
plot_3d(heightmap = mat_gl,
zscale = 25,
solid = FALSE,
shadowdepth = 10,
windowsize = c(1000,700),
background = 'grey80')
render_camera(theta = -15, phi = 40, zoom = .5)
render_highquality(
filename = "test_plot2.png",
interactive = FALSE,
light = TRUE,
width = 3000,
height = 2100,
samples = 256
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment