Skip to content

Instantly share code, notes, and snippets.

View mdsumner's full-sized avatar

Michael Sumner mdsumner

  • Integrated Digital East Antarctica, Australian Antarctic Division
  • Hobart, Australia
View GitHub Profile
library(terra)
dsn <- "/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2025.tif"

ex <- c(-1, 1, -1, 1) * pi/2 * 6378137

## generate a template spilhaus 25km grid and project the COG to that
r &lt;- project(rast(dsn), rast(ext(ex), res = 25000, crs = "+proj=spilhaus"), by_util = TRUE)
Sys.setenv("AWS_NO_SIGN_REQUEST" = "YES")

dsn <- "/vsis3/us-west-2.opendata.source.coop/tge-labs/aef/v1/annual/2023/17N/xfbc3s1k0nrye0ub4-0000008192-0000000000.tiff"
r <- rast(sprintf("vrt://%s?ovr=3", dsn))[[1]]
par(mfrow = c(2, 1))
plot(r)

plot(r)
imdsn &lt;- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=World_Imagery"
/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2025.tif

or s3

--endpoint-url=https://projects.pawsey.org.au
s3://idea-gebco-tif/GEBCO_2025.tif

New chunk read in https://github.com/hypertidy/zaro

library(zaro)
store <- zaro("virtualizarr://https://raw.githubusercontent.com/mdsumner/virtualized/refs/heads/main/remote/ocean_temp_2023.parq")
#> [zaro] opening VirtualiZarr Parquet reference store: https://raw.githubusercontent.com/mdsumner/virtualized/refs/heads/main/remote/ocean_temp_2023.parq
meta <- attr(zaro_meta(store), "consolidated")$temp
#> [zaro] found .zmetadata (Zarr V2 consolidated)
#&gt; [zaro] 11 arrays: Time, Time_bnds, average_DT, average_T1, average_T2, nv, st_edges_ocean, st_ocean, temp, xt_ocean, yt_ocean
# library(sedonadb)
# 
# # Works with files
# fgb <- system.file("files/natural-earth_cities.fgb", package = "sedonadb")
# sd_read_sf(fgb)

## worked a while already with GDAL
library(lazysf)
url <- "/vsicurl/https://github.com/apache/sedona-db/raw/refs/heads/main/r/sedonadb/inst/files/natural-earth_cities_geo.parquet"
#' peek - browse XYZ tiles in the RStudio Viewer via magick
#'
#' User provides an EPSG:3857 bounding box, peek picks an appropriate zoom,
#' fetches the tiles, composites them with magick, and displays the result.

# --- tile math (this is what grout does) ---

#' Web Mercator full extent
ORIGIN <- -20037508.342789244

Here get raw tiles and calculate median across time. We're only plotting in a faceted index space, but it's lined up relatively.

library(rustycogs)  ## hypertidy/rustycogs  - we need dev-version of async_tiff which is in the Cargo.toml (#PR 265)
library(ximage)
js <- jsonlite::fromJSON(sds::stacit("55GDN", c("2025-01-01", format(Sys.Date()))))  ## hypertidy/sds relies on recent MGRS grid code rather than extent
hrf <- unlist(lapply(js$features$assets[1:19], function(x) x[["href"]]))
all(grepl("tif$",hrf))
                     
## index all the tifs (like virtualizarr, we get the IFD tables)
# 1. Linear percent clip
# low/high are the natural UX knobs; 2%/98% is a solid default,
# tighter for low-contrast scenes.
s2_stretch_linear <- function(x, low = 0.02, high = 0.98) {
v <- terra::values(x, na.rm = TRUE)
lo <- quantile(v, low, na.rm = TRUE)
hi <- quantile(v, high, na.rm = TRUE)
terra::clamp(x, lo, hi, values = TRUE) |>
(\(r) (r - lo) / (hi - lo))()
}

Myriahedral Projections: Implementations

Myriahedral projections (van Wijk, 2008) are not projections in the PROJ sense — there's no forward/inverse formula. They're a mesh-algorithm-projection pipeline:

  1. Build a fine triangular mesh on the sphere (icosahedral subdivision, graticule grid, or geography-adaptive)
  2. Assign edge weights (land/ocean crossings, graticule alignment, etc.)
  3. Compute a minimum spanning tree of the dual graph to decide where to cut
  4. Unfold the tree of faces flat using per-face gnomonic projections

The result has negligible area and angle distortion at the cost of many interrupts.

GDAL-ready basemap tiles

Access common web basemap tile services via GDAL. Returns either WMTS connection strings (for ESRI services with WMTS support) or TMS minidriver XML for XYZ tile services.

Function

basemap <- function(name = NULL, url = NULL, api_key = NULL, tile_level = 19L, bands = 3L) {
  providers <- list(
    # OpenStreetMap (no key)
    OpenStreetMap = "https://tile.openstreetmap.org/${z}/${x}/${y}.png",