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
reticulate::use_python("/workenv/bin/python3")
reticulate::py_require("numcodecs")
numcodecs <- reticulate::import("numcodecs")

## just found by eye
zlib <- numcodecs$zlib$Zlib(4L)

https://github.com/openlandmap/GEDTM30?tab=readme-ov-file

gdalinfo /vsicurl/https://s3.opengeohub.org/global/edtm/legendtm_rf_30m_m_s_20000101_20231231_go_epsg.4326_v20250130.tif

Driver: GTiff/GeoTIFF
Files: /vsicurl/https://s3.opengeohub.org/global/edtm/legendtm_rf_30m_m_s_20000101_20231231_go_epsg.4326_v20250130.tif
Size is 1440010, 600010
Coordinate System is:
x_from_col <- function(dimension, bbox, col) {
  col[col < 1] <- NA
  col[col > dimension[1L]] <- NA
  xres <- diff(bbox[c(1, 3)]) / dimension[1]
  bbox[1] - xres/2 + col * xres
}
y_from_row <- function(dimension, bbox, row) {
  row[row < 1] <- NA
  row[row > dimension[2]] <- NA

extract rema at points

library(xml2)
library(gdalraster)
dsn <- "/vsicurl/https://raw.githubusercontent.com/mdsumner/rema-ovr/main/REMA-2m_dem_ovr.vrt"
url <- gsub("/vsicurl/", "", dsn)
xml <- read_xml(url)
dst <- xml |> xml_find_all(".//DstRect")
## https://developmentseed.org/obstore/latest/examples/fastapi/

# Example large Parquet file hosted in AWS open data
#store = S3Store("ookla-open-data", region="us-west-2", skip_signature=True)
#path = "parquet/performance/type=fixed/year=2024/quarter=1/2024-01-01_performance_fixed_tiles.parquet"


Sys.setenv("AWS_REGION" = "us-west-2")

128 cpus, 158 seconds

options(parallelly.fork.enable = TRUE, future.rng.onMisuse = "ignore")
library(furrr); plan(multicore)
d <- arrow::read_parquet("https://data.source.coop/ausantarctic/ghrsst-mur-v2/ghrsst-mur-v2.parquet")
dsn <- sprintf("/vsicurl/%s", d$assets$analysed_sst$href)

#(cell <- terra::cellFromXY(terra::rast(dsn[1]), cbind(147, -48)))
# 496796700

Open a virtual dataset, how can I get the dict() of references from .lat?

import virtualizarr
oisst = virtualizarr.open_virtual_dataset("https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198109/oisst-avhrr-v02r01.19810901.nc")

oisst.lat
<xarray.DataArray 'lat' (lat: 720)> Size: 3kB
ManifestArray

Warning, requires very recent GDAL - built on 2025-03-21 source code.

gdal raster pipeline ! \ 
    read /vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif  ! \
    reproject --resolution 1000,1000  --bbox 140,-45,150,-33 --bbox-crs "EPSG:4326"  --dst-crs "https://spatialreference.org/ref/epsg/3111" ! \
     write gebco_tas_vic.gdalg.json
library(terra)
r <- rast(matrix(sample(1:5, 60,  replace = TRUE), 5))

uv <- unique(r)[,1L, drop = TRUE]
l <- setNames(vector("list", length(uv)), sprintf("l=%i", uv))
for (i in seq_along(uv)) {
  l[[i]] <- r == uv[i]
}