mdsumner

cellcentre.offset <- function(x) {
  ex <- terra::ext(x)
  res <- terra::res(x)
  halfcell <- res/2
  cco <- c(ex[1], ex[4]) + halfcell * c(1, -1)
  cco
}

mdsumner

service <- file.path(
  "https://s3-us-west-2.amazonaws.com", 
  "mrlc"
)
layers <- paste0(
  service, 
  "/Annual_NLCD_LndCov_", 
  years, 
  "_CU_C1V0.tif"

mdsumner

dsn <- "/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif"

library(terra)
r <- rast(ext(-1, 1, -1, 1) * 6378137 * 1.9, res = 25000, crs = "+proj=spilhaus")
geb <- project(rast(dsn), r, by_util = TRUE)

mdsumner

with stars (uses GDAL MULTIDIM)

## note we have to inner-quote the string (I think ZARR driver needs some attention here)
 dsn <- "ZARR:\"/vsicurl/https://mur-sst.s3.us-west-2.amazonaws.com/zarr-v1\""
library(stars)
read_mdim(dsn, proxy = TRUE, bounds = FALSE, variable = "analysed_sst")
stars_proxy object with 1 attribute in 1 file(s):
$analysed_sst

mdsumner

library(maptiles)
library(sf)
nc = st_read(system.file("shape/nc.shp", package = "sf"))

## let's say we have target raster
rr <-  get_tiles(
    sf::st_transform(nc, "EPSG:3857"),

mdsumner

following on from https://gist.github.com/mdsumner/4ff89897e95d834073a2f49bf648707f we're working on bindings for the GDAL multidimensional api in R

## this code is very WIP and only in-dev

library(gdalraster)  ## a draft branch on mdsumner/gdalraster@multidimnew
dsn <- "/vsigs/gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3"
Sys.setenv("GS_NO_SIGN_REQUEST" = "YES")

ds <- new(GDALMultiDimRaster, dsn)

mdsumner

sf::gdal_utils("vectortranslate", "OAPIF:https://ogc-api.nrw.de/inspire-us-feuerwehr", tf <- tempfile(tmpdir = "/vsimem", fileext = ".fgb"), options = c("-f", "FlatGeobuf", "-where", "name = 'Schwelm'"))
terra::vect(tf)

class       : SpatVector
 geometry    : points
 dimensions  : 1, 22  (geometries, attributes)
 extent      : 381347.8, 381347.8, 5682950, 5682950  (xmin, xmax, ymin, ymax)

mdsumner

import s3fs
bucket = 'podaac-ops-cumulus-protected'

creds = "<the literal json available at https://archive.podaac.earthdata.nasa.gov/s3credentials>"

s3 = s3fs.S3FileSystem(key=creds['accessKeyId'], 
                       secret=creds['secretAccessKey'],
                       token=creds['sessionToken'],
 client_kwargs={'region_name':'us-west-2'})

mdsumner

Crop out NZ from COP30 or GEBCO:

export dsn=/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt
gdalwarp $dsn nz.tif -te 165 -51 180 -33 -ts 1024 0 out.tif

if you want bathymetry as well use GEBCO (will also be a lot quicker, so maybe this and then mask, see gdal_calc.py in linux)

mdsumner

Rely on the geolocation arrays and GDAL warp:

gdalwarp -te -5500000 -5500000 5500000 5500000 -t_srs "+proj=geos +lon_0=140.7 +h=35785863 +x_0=0 +y_0=0 +a=6378137 +b=6356752.3 +units=m +no_defs" vrt://AHI-CHGT_v1r1_h09_s202403040030204_e202403040039398_c202403040047497.nc?sd_name=CldTopTemp test.vrt

or, assign the grid and ignore the geolocation arrays (in this case we also need GDAL_NETCDF_BOTTOMUP=NO)