Pre-process

A slightly longer title

Author

James Goldie

Published

May 2, 2023

Observed and interpolated subsistence rates are available for each city as geoTIFFs from NTU Dataverse.

Although some Dataverse instances have an API, I couldn’t verify that this one does. Instead, you can download the dataset manually and place it in /data as dataverse_files.zip.

Code

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

Code

library(terra)

terra 1.7.29

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract

Code

library(tidyterra)


Attaching package: 'tidyterra'

The following object is masked from 'package:stats':

    filter

Code

library(here)

here() starts at /workspaces/report-sinking-cities

Let’s first unzip the dataverse files:

Code

# unzip the downloaded dataverse zip
unzip(
  here("data", "dataverse_files.zip"),
  exdir = here("data", "src"))

# unzip each city
list.files(here("data", "src"), pattern = glob2rx("*.zip"),
  full.names = TRUE) %>%
  walk(unzip, exdir = here("data", "src"))

# delete the city zip files
unlink(here("data", "src", "*.zip"), expand = TRUE)

Each city has four geoTIFFs:

velocity_InSAR: the subsistence velocity observed using InSAR;
velocity_interpolation: the subsistence velocity interpolated between observations;
velocityStd_InSAR: the standard deviation of the observed velocity; and
velocityStd_interpolation: the standard deviation of the interpolated velocity.

Note

the velocities and standard deviations of them are both in units of m yr^-1. The paper’s figures present them in units of mm yr^-1.

For each city, I want to:

Combine the two velocity maps and the two SD maps;
Calculate a normalised standard deviation
Mask the combined velocity map using bins of the combined SD map.

Code

tibble(
  path = list.files(here("data", "src"), pattern = glob2rx("*.tif"),
    full.names = TRUE),
  fname = str_remove(basename(path), ".tif")) %>%
  # discard the area bmsl files
  filter(str_detect(fname, "areaBMSL", negate = TRUE)) %>%
  # extract filename info
  separate(fname,
    into = c("city_name", "city_code", "measure", "source"),
    sep = "_") ->
city_maps

We’ll need to use {terra} or a similar raster package to do the above two tasks.

First, we need functions to do the two tasks above:

Code

# combine_maps: load two geotiff paths in, merge and write back to disk.
# return the output path
combine_maps <- function(path_a, path_b, city, measure) {
  out_path <- here("data", "1-combined", paste0(city, "_", measure, ".tif"))

  merge(rast(path_a), rast(path_b), filename = out_path, overwrite = TRUE)
  return(out_path)
}

Code

get_norm_sd <- function(velocity_path, sd_path, city) {
  out_path <- here("data", "1-combined",
    paste0(city, "_velocityNormStd", ".tif"))

  norm_sd <- rast(sd_path) / rast(velocity_path)
  writeRaster(norm_sd, out_path, overwrite = TRUE)
  return(out_path)
}

Code

# mask_map: given a velocity raster and a normalised std deviation path,
# clip the velocities to only contain cells with normalised SD within an
# absolute range of low to high. return the path of the result (written to disk)
mask_map <- function(velocity_path, normsd_path, low, high, city, thresh_label) {
  out_path <- here("data", "2-masked",
    paste0(city, "_", thresh_label, "p.tif"))

  vel <- rast(velocity_path)
  normsd <- rast(normsd_path)
  
  # get an absolute value of the normalised sd
  abs_normsd <- ifel(normsd < 0, -normsd, normsd)

  # first mask the normalised sd to within the absolute threshold
  normsd_thresh <- ifel(abs_normsd > high | abs_normsd < low, NA, abs_normsd)

  # then mask the velocity map based on the remaining normalised sd map
  mask(vel, normsd_thresh, filename = out_path, overwrite = TRUE)

  return(out_path)
}

Code

dir.create(here("data", "1-combined"), showWarnings = FALSE)
dir.create(here("data", "2-masked"), showWarnings = FALSE)

city_maps %>%
  group_by(city_name, measure) %>%
  # TODO - assert two rows per group
  summarise(
    path = combine_maps(path[1], 
      path[2], city_name[1], measure[1])) %>%
  ungroup() ->
city_maps_combined

`summarise()` has grouped output by 'city_name'. You can override using the
`.groups` argument.

Code

# calculate the normalised standard deviations
city_maps_combined %>%
  # ensure velocity is first in each group and sd is second
  arrange(measure) %>%
  group_by(city_name) %>%
  summarise(path = get_norm_sd(path[1], path[2], city_name[1])) %>%
  ungroup() %>%
  mutate(measure = "velocityNormStd") ->
city_norm_sds

# combine the normalised standard velocities with the velocities (drop the sds)
city_maps_combined %>%
  filter(measure == "velocity") %>%
  bind_rows(city_norm_sds) %>%
  arrange(city_name, measure) ->
city_maps_all

# now do the masking, in bins of normalised standard deviation:
# 0-25%, 25-50%, 50-1005, > 100%
city_maps_all %>%
  group_by(city_name) %>%
  summarise(
    masked_path_0to25p =
      mask_map(path[1], path[2], 0,    0.25, city_name[1], "0to25p"),
    masked_path_25to50p =
      mask_map(path[1], path[2], 0.25, 0.5,  city_name[1], "25to50p"),
    masked_path_50to100p =
      mask_map(path[1], path[2], 0.5,  1,    city_name[1], "50to100p"),
    masked_path_gt100p =
      mask_map(path[1], path[2], 1,    Inf,  city_name[1], "gt100p")) ->
city_maps_masked

We also need to convert these TIFFs to PNGs for web display, and we need to run off their extents (along with the reference point KMLs) as a data frame to help us fit the PNGs on the map.

But before we do that, we’re going to reproject everything in data/2-masked to EPSG:3857 to see if it helps us reconcile Mapbox and DeckGL:

Code

# reproject_to_3857: read in a raster path and write back out reprojected
reproject_to_3857 <- function(path) {
  out_path <- here("data", "2a-masked-3857", basename(path))
  project(rast(path), "epsg:3857", gdal = FALSE, filename = out_path,
    overwrite = TRUE)
  return(out_path)
}

We’ll use this function in a moment, before we get the map bounds (reprojected).

First, let’s get the reference points:

Code

tibble(
  path = list.files(here("data", "src"), pattern = glob2rx("*.kml"),
    full.names = TRUE),
  fname = basename(path)) %>%
  separate(fname, into = c("city", "code", "ref"), sep = "_", remove = FALSE) %>%
  select(-code, -ref) %>%
  group_by(city) %>%
  # read kml files in
  summarise(point = st_read(path)) %>%
  unpack(point) %>%
  # extract lat and lon
  mutate(
    coord = st_coordinates(geometry),
    ref_lon = coord[, "X"],
    ref_lat = coord[, "Y"]) %>%
  select(city, ref_lon, ref_lat) %>%
  print() ->
reference_points

Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Abidjan_TA118_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -4.04875 ymin: 5.45625 xmax: -4.04875 ymax: 5.45625
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Ahmedabad_TD034_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 72.66875 ymin: 23.06875 xmax: 72.66875 ymax: 23.06875
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Alexandria_TD065_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 29.81792 ymin: 30.99958 xmax: 29.81792 ymax: 30.99958
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Bangkok_TD062_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 100.6546 ymin: 13.76625 xmax: 100.6546 ymax: 13.76625
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Barcelona_TA132_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 2.13375 ymin: 41.41125 xmax: 2.13375 ymax: 41.41125
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/BuenosAires_TD170_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -58.52125 ymin: -34.61208 xmax: -58.52125 ymax: -34.61208
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Chennai_TD092_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 80.11792 ymin: 13.18375 xmax: 80.11792 ymax: 13.18375
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Chittagong_TD077_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 91.83208 ymin: 22.37458 xmax: 91.83208 ymax: 22.37458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Dalian_TA098_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 121.9987 ymin: 39.32875 xmax: 121.9987 ymax: 39.32875
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/DaresSalaam_TD006_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 39.16292 ymin: -6.837917 xmax: 39.16292 ymax: -6.837917
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Dhaka_TD150_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 90.42208 ymin: 23.80792 xmax: 90.42208 ymax: 23.80792
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Dongguan_TA011_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 114.0354 ymin: 22.95792 xmax: 114.0354 ymax: 22.95792
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Foshan_TA011_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 113.0962 ymin: 23.09625 xmax: 113.0962 ymax: 23.09625
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Fukuoka_TD163_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 130.6696 ymin: 33.61875 xmax: 130.6696 ymax: 33.61875
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Guangzhou_TA011_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 113.3004 ymin: 23.23625 xmax: 113.3004 ymax: 23.23625
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Hangzhou_TA069_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 120.2229 ymin: 29.96792 xmax: 120.2229 ymax: 29.96792
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/HoChiMinhCity_TD018_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 106.5921 ymin: 10.81708 xmax: 106.5921 ymax: 10.81708
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/HongKong_TA011_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 114.1729 ymin: 22.51542 xmax: 114.1729 ymax: 22.51542
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Houston_TD143_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -95.27125 ymin: 29.79958 xmax: -95.27125 ymax: 29.79958
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Istanbul_TD138_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 28.87042 ymin: 41.08625 xmax: 28.87042 ymax: 41.08625
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Jakarta_TA098_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 106.8754 ymin: -6.34125 xmax: 106.8754 ymax: -6.34125
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Karachi_TD078_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 67.15625 ymin: 24.93458 xmax: 67.15625 ymax: 24.93458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Kolkata_TD048_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 88.47875 ymin: 22.72208 xmax: 88.47875 ymax: 22.72208
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Lagos_TD095_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 3.30625 ymin: 6.62125 xmax: 3.30625 ymax: 6.62125
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Lima_TD069_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -76.95958 ymin: -12.04958 xmax: -76.95958 ymax: -12.04958
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/London_TD081_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -0.07458351 ymin: 51.54708 xmax: -0.07458351 ymax: 51.54708
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/LosAngeles_TD071_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -118.2438 ymin: 34.04292 xmax: -118.2438 ymax: 34.04292
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Luanda_TA131_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 13.40458 ymin: -9.064583 xmax: 13.40458 ymax: -9.064583
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Manila_TD032_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 120.9929 ymin: 14.70208 xmax: 120.9929 ymax: 14.70208
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Miami_TA048_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -80.29625 ymin: 25.76542 xmax: -80.29625 ymax: 25.76542
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Mumbai_TD034_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 72.84542 ymin: 19.18292 xmax: 72.84542 ymax: 19.18292
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Nagoya_TD119_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 137.0104 ymin: 35.17458 xmax: 137.0104 ymax: 35.17458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Nanjing_TA069_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 119.0237 ymin: 31.64292 xmax: 119.0237 ymax: 31.64292
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/NewYork_TA033_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -73.85958 ymin: 40.72792 xmax: -73.85958 ymax: 40.72792
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Osaka_TD017_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 135.6787 ymin: 34.78958 xmax: 135.6787 ymax: 34.78958
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Philadelphia_TA106_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -75.15958 ymin: 40.06292 xmax: -75.15958 ymax: 40.06292
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Qingdao_TD076_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 120.0229 ymin: 36.75458 xmax: 120.0229 ymax: 36.75458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/RiodeJaneiro_TD155_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -43.34958 ymin: -22.80792 xmax: -43.34958 ymax: -22.80792
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/SaintPetersburg_TD007_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 30.31542 ymin: 59.86292 xmax: 30.31542 ymax: 59.86292
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Seoul_TD134_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 126.9987 ymin: 37.59042 xmax: 126.9987 ymax: 37.59042
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Shanghai_TA171_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 121.2712 ymin: 31.02792 xmax: 121.2712 ymax: 31.02792
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Singapore_TD018_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.8737 ymin: 1.377917 xmax: 103.8737 ymax: 1.377917
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Surat_TD034_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 72.87542 ymin: 21.17458 xmax: 72.87542 ymax: 21.17458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Suzhou_TA171_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 120.5137 ymin: 31.31458 xmax: 120.5137 ymax: 31.31458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Tianjin_TD149_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 117.2537 ymin: 39.75875 xmax: 117.2537 ymax: 39.75875
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Tokyo_TD046_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 139.3146 ymin: 35.77458 xmax: 139.3146 ymax: 35.77458
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Washington,D.C._TA004_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: -77.02292 ymin: 38.90625 xmax: -77.02292 ymax: 38.90625
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Reading layer `referencepoint' from data source 
  `/workspaces/report-sinking-cities/data/src/Yangon_TD033_referencePoint.kml' 
  using driver `LIBKML'
Simple feature collection with 1 feature and 12 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 96.13542 ymin: 16.95708 xmax: 96.13542 ymax: 16.95708
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
# A tibble: 48 × 3
   city        ref_lon ref_lat
   <chr>         <dbl>   <dbl>
 1 Abidjan       -4.05    5.46
 2 Ahmedabad     72.7    23.1 
 3 Alexandria    29.8    31.0 
 4 Bangkok      101.     13.8 
 5 Barcelona      2.13   41.4 
 6 BuenosAires  -58.5   -34.6 
 7 Chennai       80.1    13.2 
 8 Chittagong    91.8    22.4 
 9 Dalian       122.     39.3 
10 DaresSalaam   39.2    -6.84
# ℹ 38 more rows

Now, here’s a function to extract the domain and range of the raster file:

Code

# get_range_bounds: given a raster path, return a list containing:
#   - the `min` and `max` values of the raster band
#   - the 
get_range_bounds <- function(path) {
  
  # open raster and get summary and extent
  vel <- rast(path)
  vel_summary <- summary(vel)

  # extract raster range from summary
  vel_min <-
    vel_summary %>% pluck(1) %>% str_remove("Min.   :") %>% as.numeric()
  vel_max <-
    vel_summary %>% pluck(6) %>% str_remove("Max.   :") %>% as.numeric()

  return(list(
    min = vel_min,
    max = vel_max,
    xmin = xmin(vel),
    xmax = xmax(vel),
    ymin = ymin(vel),
    ymax = ymax(vel)))
}

Let’s bring all that info together:

Code

dir.create(here("data", "2a-masked-3857"), showWarnings = FALSE)

city_maps_masked %>%
  pivot_longer(starts_with("masked"), names_to = "measure",
    values_to = "path") %>%
  mutate(reprojected_path = map_chr(path, reproject_to_3857)) ->
city_maps_reprojected

city_maps_reprojected %>%
  mutate(
    range_bounds_latlon = map(path, get_range_bounds),
    range_bounds_3857 = map(reprojected_path, get_range_bounds),
    measure = str_remove(measure, "masked_path_"),
    filename = basename(reprojected_path)) %>%
  unnest_wider(range_bounds_latlon) %>%
  rename(lon_min = xmin, lon_max = xmax, lat_min = ymin, lat_max = ymax) %>%
  select(-min, -max) %>%
  unnest_wider(range_bounds_3857) %>%
  # finally, merge in reference point data
  left_join(reference_points, by = join_by(city_name == city),
    multiple = "all") ->
city_maps_info

Warning: There were 328 warnings in `mutate()`.
The first warning was:
ℹ In argument: `range_bounds_latlon = map(path, get_range_bounds)`.
Caused by warning:
! [summary] used a sample
ℹ Run `dplyr::last_dplyr_warnings()` to see the 327 remaining warnings.

We want to ensure all the maps are made using a common, symmetrical colour scale. To do that, I’ll get a range that covers all of the cities’ data ranges, and I’ll push it out so that it’s symmetrical around zero:

Code

dir.create(here("data", "3-pngs"), showWarnings = FALSE)

# first, we need to settle on a common range that encompasses all the files
# (and make it symmetrical around 0)
# data_range <- c(
#   min(city_maps_info$min, na.rm = TRUE),
#   max(city_maps_info$max, na.rm = TRUE)) %>%
#   abs() %>%
#   max()

# compute breaks based on that range
# number_colours <- 10
# palette_breaks <- seq(-data_range, data_range, length.out = number_colours + 1)

Here’s our plotting function:

Code

# plot_raster_png: save the raster as a png (using tidyterra to make applying
# a common colour scale across files much easier).
# note: we really need to ensure that the png's aspect ratio reflects the extent
# of the tiff in order to get decent projection on the interactive map. i've
# reprojected the tiff to web mercator (epsg:3857) JUST for apsect ratio
# calculation.
plot_raster_png <- function(path, palette_breaks) {

  vel <- rast(path)
  vel_3857 <- project(vel, "epsg:3857")

  # this might not work well because of polar distortion, but in lieu of
  # knowing the actual panel aspect ratio, i'm going to work it out from the
  # raster extent and see if it's a better fit than these square plots with
  # lots of padding
  vel_ext_3857 <- ext(vel_3857)
  aspect_ratio_3857 <-
    (vel_ext_3857[2] - vel_ext_3857[1]) /
    (vel_ext_3857[4] - vel_ext_3857[3])

  # use zero scale expansion + void theme to remove margins
  vel_plot <- ggplot() +
    geom_spatraster(data = vel) +
    scale_fill_fermenter(
      type = "div",
      palette = "RdYlBu",
      direction = 1,
      na.value = NA,
      # use `breaks` and `limits` to supply our common breaks
      # (and force all of them to be used)
      breaks = palette_breaks,
      limits = range(palette_breaks),
      guide = NULL
      ) +
    # push plot out to edges of the panel
    scale_x_continuous(expand = expansion(0, 0)) +
    scale_y_continuous(expand = expansion(0, 0)) +
    # remove backgrounds and other elements
    theme_void() +
    theme(
      plot.background = element_blank(),
      panel.background = element_blank())


  # write to disk, varying the plot height based on the epsg:3857 aspect ratio
  ggsave(
    here("data", "3-pngs", str_replace(basename(path), ".tif$", ".png")),
    vel_plot,
    width = 800,
    height = 800 / aspect_ratio_3857,
    units = "px")

}

Let’s do the export on each processed TIF now:

Code

number_colours <- 10

city_maps_info %>%
  # compute city-specific ranges for the colour palettes
  group_by(city_name) %>%
  mutate(
    city_min = min(min, na.rm = TRUE),
    city_max = max(max, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    palette_max_abs = pmax(abs(city_min), abs(city_max)),
    city_breaks = map(palette_max_abs,
      ~ seq(-.x, .x, length.out = number_colours + 1)),
    walk2(path, city_breaks, plot_raster_png)
    ) %>%
  # export to csv (dropping paths)
  select(filename, ref_lon, ref_lat, file_min = min, file_max = max,
    palette_extent = palette_max_abs, lon_min, lon_max, lat_min, lat_max, xmin,
    xmax, ymin, ymax) %>%
  write_csv(here("data", "map-bounds-extent.csv"))

SpatRaster resampled to ncells = 500976
SpatRaster resampled to ncells = 500976
SpatRaster resampled to ncells = 500976
SpatRaster resampled to ncells = 500976

SpatRaster resampled to ncells = 501320
SpatRaster resampled to ncells = 501320
SpatRaster resampled to ncells = 501320
SpatRaster resampled to ncells = 501320

SpatRaster resampled to ncells = 500465
SpatRaster resampled to ncells = 500465
SpatRaster resampled to ncells = 500465
SpatRaster resampled to ncells = 500465

SpatRaster resampled to ncells = 500588
SpatRaster resampled to ncells = 500588
SpatRaster resampled to ncells = 500588
SpatRaster resampled to ncells = 500588

SpatRaster resampled to ncells = 500535
SpatRaster resampled to ncells = 500535
SpatRaster resampled to ncells = 500535
SpatRaster resampled to ncells = 500535

SpatRaster resampled to ncells = 500308
SpatRaster resampled to ncells = 500308
SpatRaster resampled to ncells = 500308
SpatRaster resampled to ncells = 500308

SpatRaster resampled to ncells = 501034
SpatRaster resampled to ncells = 501034
SpatRaster resampled to ncells = 501034
SpatRaster resampled to ncells = 501034

SpatRaster resampled to ncells = 500736
SpatRaster resampled to ncells = 500736
SpatRaster resampled to ncells = 500736
SpatRaster resampled to ncells = 500736

SpatRaster resampled to ncells = 500520
SpatRaster resampled to ncells = 500520
SpatRaster resampled to ncells = 500520
SpatRaster resampled to ncells = 500520

SpatRaster resampled to ncells = 500858
SpatRaster resampled to ncells = 500858
SpatRaster resampled to ncells = 500858
SpatRaster resampled to ncells = 500858

SpatRaster resampled to ncells = 500797
SpatRaster resampled to ncells = 500797
SpatRaster resampled to ncells = 500797
SpatRaster resampled to ncells = 500797

SpatRaster resampled to ncells = 501095
SpatRaster resampled to ncells = 501095
SpatRaster resampled to ncells = 501095
SpatRaster resampled to ncells = 501095

SpatRaster resampled to ncells = 500760
SpatRaster resampled to ncells = 500760
SpatRaster resampled to ncells = 500760
SpatRaster resampled to ncells = 500760

SpatRaster resampled to ncells = 501808
SpatRaster resampled to ncells = 501808
SpatRaster resampled to ncells = 501808
SpatRaster resampled to ncells = 501808

SpatRaster resampled to ncells = 500832
SpatRaster resampled to ncells = 500832
SpatRaster resampled to ncells = 500832
SpatRaster resampled to ncells = 500832

SpatRaster resampled to ncells = 501147
SpatRaster resampled to ncells = 501147
SpatRaster resampled to ncells = 501147
SpatRaster resampled to ncells = 501147

SpatRaster resampled to ncells = 500904
SpatRaster resampled to ncells = 500904
SpatRaster resampled to ncells = 500904
SpatRaster resampled to ncells = 500904

SpatRaster resampled to ncells = 500400
SpatRaster resampled to ncells = 500400
SpatRaster resampled to ncells = 500400
SpatRaster resampled to ncells = 500400

SpatRaster resampled to ncells = 500536
SpatRaster resampled to ncells = 500536
SpatRaster resampled to ncells = 500536
SpatRaster resampled to ncells = 500536

SpatRaster resampled to ncells = 500488
SpatRaster resampled to ncells = 500488
SpatRaster resampled to ncells = 500488
SpatRaster resampled to ncells = 500488

SpatRaster resampled to ncells = 501160
SpatRaster resampled to ncells = 501160
SpatRaster resampled to ncells = 501160
SpatRaster resampled to ncells = 501160

SpatRaster resampled to ncells = 500703
SpatRaster resampled to ncells = 500703
SpatRaster resampled to ncells = 500703
SpatRaster resampled to ncells = 500703

SpatRaster resampled to ncells = 500816
SpatRaster resampled to ncells = 500816
SpatRaster resampled to ncells = 500816
SpatRaster resampled to ncells = 500816