Skip to content

Instantly share code, notes, and snippets.

View mdsumner's full-sized avatar

Michael Sumner mdsumner

598 followers · 1.5k following
  • Integrated Digital East Antarctica, Australian Antarctic Division
  • Hobart, Australia
View GitHub Profile
@mdsumner
mdsumner / Image-servers.md
Created March 4, 2026 10:59

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",
@mdsumner
mdsumner / laea-vs-mercator.md
Last active March 1, 2026 11:33 
library(wk)
library(geos)
library(PROJ)
## flatten all to POLYGON
map <- wk_flatten(as_wkb(rnaturalearth::ne_countries(scale = 50)))
local_laea <- function(x, merctrans) {
  cent <- geos_centroid(x)
  mercxy <- wk_coords(wk_transform(cent, merctrans))[,c("x", "y")]
 xy &lt;- wk_coords(cent)[, c("x", "y")]
@mdsumner
mdsumner / Gdal-caching-range-requests.md
Created February 28, 2026 09:15

Caching Range Requests: What GDAL Already Does (and the Hard Part It Doesn't)

The question of how to cache HTTP range requests — particularly for cloud-native geospatial formats like COG, PMTiles, FlatGeobuf, and cloud-optimized GeoParquet — keeps coming up. Brandon Liu's How many ranges can you fit in one request is a good treatment of the multi-range packing problem. But there's a mature, battle-tested system that already handles much of this at the client level, and its design choices are instructive even for people building entirely different stacks: GDAL's /vsicurl/ virtual filesystem.

GDAL's VSI curl layer and caching

When GDAL reads a cloud-optimized file via /vsicurl/ (or its cloud-specific variants /vsis3/, /vsigs/, /vsiaz/), it performs range request management internally. The behaviour is controlled by a set of configuration options that most users never touch, but that encode a lot of hard-won knowledge about how to read efficie

@mdsumner
mdsumner / Nazis-run-your-town.md
Created February 25, 2026 00:53

this is fine

@mdsumner
mdsumner / gist:675b3cb501ffead0e5fe9b94d3c9946a
Created February 21, 2026 22:53
Nut
@mdsumner
mdsumner / helpers.R
Last active February 25, 2026 22:02
##' ESRI Basemap Services
##'
##' Access ESRI ArcGIS Online basemap services via WMTS (preferred) or TMS XML fallback.
##'
##' @param name Character. Name of the basemap service. Use `esri_list()` to see available options.
##' @param url Character. For custom services, provide the MapServer base URL
##' (e.g., "https://services.arcgisonline.com/arcgis/rest/services/Ocean_Basemap/MapServer").
##' @param tile_level Integer. Maximum tile level for TMS XML (default 19).
##' @param bands Integer. Number of bands - 3 for RGB, 4 for RGBA (default 3).
##'
@mdsumner
mdsumner / test-gdal-zarr-v3-datetime64.sh
Created February 12, 2026 01:43
#!/usr/bin/env bash
# test_gdal_zarr_v3_datetime64.sh
#
# Reproducer for GDAL Zarr V3 numpy.datetime64 extension data type issue.
# Creates a minimal local Zarr V3 store with a datetime64 time coordinate,
# then exercises GDAL classic 2D raster API access patterns to confirm
# which operations are affected by the unsupported extension type.
#
# Requirements:
# - GDAL >= 3.8 (with Zarr V3 support)
@mdsumner
mdsumner / gdal-zarr-v3-numpy_datetime64.md
Last active February 12, 2026 01:52
@mdsumner
mdsumner / crop_to_footprint_b4_warp.md
Last active February 7, 2026 11:52

https://fosstodon.org/@hughagraham/116005729038262380

## generate a tif with a nodata collar
dsn <- "/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif"
library(terra)
#> terra 1.8.96
pr <- rast(ext(c(-1, 1, -1, 1)) * 1e7, res = 4000, crs = "+proj=laea")
pr
@mdsumner
mdsumner / gdalraster_nara.md
Created January 23, 2026 11:51
Native raster plot from a direct GDAL call

Using new nativeRaster support in gdalraster from github, pull a reprojected image from a tile server direct to R's internal image encoding"

library(gdalraster)
#> GDAL 3.13.0dev-0e4f0cadf6 (released 2026-01-20), GEOS 3.12.1, PROJ 9.7.0
dsn <- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml,layer=World_Imagery"

ex <- reproj::reproj_extent(c(143, 149, -44, -39), "EPSG:3031", source = "EPSG:4326")

fit_dims <- function(size = 1024L, wh = c(size, size)) {