3d

Welcome to #30DayMapChallenge 2023 day 23

3d Italy Population Density Map: this map shows the population density of Italy in 3D spikes.
3D
Italy
Population Density
Published

November 23, 2023

Overview

This is the 3D Italy’s Population Density Map made with the instruction provided by Milos Agathon: Milos Makes Maps tutorial:

3d Italy Population Density Map
library("tidyverse")
library("R.utils")
library("httr")
library("sf")
library("stars")
library("rayshader")

Download and unzip kontur data:

url <-
    "https://geodata-eu-central-1-kontur-public.s3.amazonaws.com/kontur_datasets/kontur_population_IT_20220630.gpkg.gz"
file_name <- "italy-population.gpkg.gz"

get_population_data <- function() {
    res <- httr::GET(
        url,
        write_disk(file_name),
        progress()
    )

    R.utils::gunzip(file_name, remove = F)
}

get_population_data()
load_file_name <- gsub(".gz", "", file_name)
crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
get_population_data <- function() {
    pop_df <- sf::st_read(
        load_file_name
    ) |>
        sf::st_transform(crs = crsLONGLAT)
}

pop_sf <- get_population_data()

head(pop_sf)

First raw image:

ggplot() +
    geom_sf(
        data = pop_sf,
        color = "grey10", fill = "grey10"
    )
# ggsave("raw.png")

Make it a raster

bb <- sf::st_bbox(pop_sf)

get_raster_size <- function() {
    height <- sf::st_distance(
        sf::st_point(c(bb[["xmin"]], bb[["ymin"]])),
        sf::st_point(c(bb[["xmin"]], bb[["ymax"]]))
    )
    width <- sf::st_distance(
        sf::st_point(c(bb[["xmin"]], bb[["ymin"]])),
        sf::st_point(c(bb[["xmax"]], bb[["ymin"]]))
    )

    if (height > width) {
        height_ratio <- 1
        width_ratio <- width / height
    } else {
        width_ratio <- 1
        height_ratio <- height / width
    }

    return(list(width_ratio, height_ratio))
}
width_ratio <- get_raster_size()[[1]]
height_ratio <- get_raster_size()[[2]]

size <- 3000
width <- round((size * width_ratio), 0)
height <- round((size * height_ratio), 0)

get_population_raster <- function() {
    pop_rast <- stars::st_rasterize(
        pop_sf |>
            dplyr::select(population, geom),
        nx = width, ny = height
    )

    return(pop_rast)
}

pop_rast <- get_population_raster()

Second raw image this time as a raster:

plot(pop_rast)
# ggsave("raw2.png")
pop_mat <- pop_rast |>
    as("Raster") |>
    rayshader::raster_to_matrix()
cols <- rev(c(
    "#0b1354", "#283680",
    "#6853a9", "#c863b3"
))

texture <- grDevices::colorRampPalette(cols)(256)

Create the initial 3D object

pop_mat |>
    rayshader::height_shade(texture = texture) |>
    rayshader::plot_3d(
        heightmap = pop_mat,
        solid = F,
        soliddepth = 0,
        zscale = 15,
        shadowdepth = 0,
        shadow_darkness = .95,
        windowsize = c(800, 800),
        phi = 65,
        zoom = .65,
        theta = -30,
        background = "white"
    )

Adjust the view after building the window object

rayshader::render_camera(phi = 75, zoom = .7, theta = 0)

rayshader::render_highquality(
    filename = "italy_population_2022.png",
    preview = T,
    light = T,
    lightdirection = 225,
    lightaltitude = 60,
    lightintensity = 400,
    interactive = F,
    width = width, height = height
)
Back to top