Skip to content

Instantly share code, notes, and snippets.

View mdsumner's full-sized avatar

Michael Sumner mdsumner

View GitHub Profile

This local example in the Usage for VirtualiZarr is (still) broken:

https://virtualizarr.readthedocs.io/en/latest/usage.html#__tabbed_1_7

I personally find it impossible to remember all this rigmarole (what's wrong with a file path?), so I need a reference example that actually works:

## local file 'data.nc' in your wd e.g. 
# wget -O data.nc https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198109/oisst-avhrr-v02r01.19810901.nc
+proj=stere +a=6378137.0 +b=6356752.3 +lat_0=-90 +lon_0=64 219383.7 309876.1 2473114 2944041

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)