new gdal mdim mosaic, with 180 Bluelink daily ocean_temp netcdf files from Thredds
## obtain a list of Thredds urls (we previously cached, we grep out daily_ocean_temp, 180 NetCDFs)
files <- arrow::read_parquet("https://github.com/mdsumner/dryrun/releases/download/latest/BRAN-netcdf-2023-netcdf.parquet") |>
dplyr::filter(stringr::str_detect(url, ".*daily.*temp.*\\.nc$"))
## generate gdal mdim mosaic command and run, a few minutes
maybe fix for terra::stretch(histeq = TRUE)
stretch_hist <- function(x, ...) {
## stretch as if all the pixels were in the same band (not memory safe)
rv <- terra::stretch(terra::rast(matrix(terra::values(x))), histeq = TRUE)
## set the values to the input, then stretch to 0,255
terra::stretch(terra::setValues(x, c(terra::values(rv))), histeq = FALSE)
}
Just an example that was mentioned, I've been playing with vrtility today looking at other parts of the world too and learning how it works.
A lot to say, this is a pretty big area for Sentinel (~180x220km), which is 5561, 4483 pixels at 40m resolution (so reduce the size of extend_x, extend_y (in degrees), and move the POINT for where you want the centre. Can go as fine as 10m resolution.
You don't really need median calc done if there isn't a lot of cloud, and we don't necessarily need mask either, date can be as fine as a two-day window (but there's only imagery every few days, depending).
library(vrtility)
# La Pampa
x0 <- c(xmin = -65.5, ymin = -39.5, xmax = -62, ymax = -34)
This url has a netcdf:
"https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc"
it has correct georeferencing (bbox 0, -90, 360, 90) but doesn't have the crs set ("EPSG:4326" is appropriate)
it also has multiple variables, which traditionally requires this verbose "subdataset syntax" for GDAL:
terra::vect("/vsicurl/https://github.com/mdsumner/geoboundaries/releases/download/latest/geoBoundariesCGAZ_ADM0.parquet")
class : SpatVector
geometry : polygons
dimensions : 218, 3 (geometries, attributes)
extent : -180, 180, -90, 83.63339 (xmin, xmax, ymin, ymax)
source : geoBoundariesCGAZ_ADM0.parquet
/vsizip//vsicurl/https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.zip
/vsizip//vsicurl/listdata.thelist.tas.gov.au/opendata/data/LIST_PARCELS_HOBART.zip/list_parcels_hobart.shp
/vsizip//vsicurl/https://github.com/AustralianAntarcticDivision/rema.proc/raw/refs/heads/master/01_rock_classification/Medium_resolution_vector_polygons_of_Antarctic_rock_outcrop_-_VERSION_7.3.zip/Medium resolution vector polygons of Antarctic rock outcrop - VERSION 7.3/add_rock_outcrop_medium_res_polygon_v7.3.gpkg
Please, I think I'm doing something wrong here. With open_virtual_mfdataset( , parallel = False, ) this runs in about 2 minutes.
With open_virtual_mfdataset( , parallel = ThreadPoolExecutor, , ) it takes the same length of time (system with 32 cpus).
import xarray as xr
from obstore.store import from_url
import time
from virtualizarr import open_virtual_dataset, open_virtual_mfdataset
I was delighted to learn that flexible coordinates in xarray are fully operational.
Following on from this (unfinished) blog post on unnecessary netcdf longlat coordinates, I now have a neat xarray/netcdf native fix for the problem I was discussing, it is a direct match to the "parse_coordinates": False example in the xarray indexes gallery.
This only a few weeks ago was challenging for me to present to non-R audiences, so I'm very excited that we can now fix these broken grids in a way more accessible to Python communities. Degenerate coordinates represent a huge entropy problem in array metadata, and (I think) we need a coordinated effort to be able to easily assign a fix
(like in GDAL with vrt://{dsn}?a_gt=c,a,b,f,d,e) and have that knowledge also feed back up to the providers to stem the flow of problematic coordinates.
mk_mth <- function(x) {
## replace daily with month
xx <- stringr::str_replace(x, "daily", "month")
## replace _YEAR with mth_YEAR
stringr::str_replace(xx, "(_[0-9]{4})", "_mth\\1")
}