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
}
service <- file.path(
"https://s3-us-west-2.amazonaws.com",
"mrlc"
)
layers <- paste0(
service,
"/Annual_NLCD_LndCov_",
years,
"_CU_C1V0.tif"
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)
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
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"),
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)
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)
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'})
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)
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
)
NewerOlder