Michael Sumner mdsumner

Bare minimum capability and library alignment for using GDAL R and Python together, using mdim or classic raster:

GDAL R and Python alignment

Single GDAL, single PROJ, single GEOS, R packages and Python packages linking the same libraries, reticulate bridging cleanly between them, modern numpy, modern Python, modern R, current GDAL master available or latest release.

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))()
	}