mdim_oisst.md

1. List all of OISST netcdf files from object storage and derive date from path
2. Remove duplicates that are 'preliminary' files replaced by final, and sort by date
3. Stack into MDIM VRT with xml and templating from the first few files

### 1. 
Sys.setenv("AWS_NO_SIGN_REQUEST" = "YES")
root <- "/vsis3/noaa-cdr-sea-surface-temp-optimum-interpolation-pds/data/v2.1/avhrr"
## root can be replaced with
## /vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/
f <- gdalraster::vsi_read_dir(root, recursive = TRUE)
d <- as.Date(stringr::str_extract(f, "[0-9]{8}"), "%Y%m%d")

### 2. 
paths <- tibble::tibble(date = d, source = f)
## drop any preliminary that also have a final
paths <- paths |> dplyr::filter_out(grepl("preliminary", source) & duplicated(date)) |> 
  dplyr::arrange(date) |> dplyr::mutate(source = sprintf("%s/%s", root, source))

range(diff(paths$date))
#Time differences in days
#[1] 1 1
diff(range(paths$date))
#Time difference of 16371 days
nrow(paths)
#[1] 16372

### 3. 
## remotes::install_github("hypertidy/vrtstack")
## create one VRT from input files, template on the date-in-path
vrt <- vrtstack::vrtstack(paths$source, concat = "(\\d{8})", parse_format = "%Y%m%d",
                origin = "1978-01-01", unit = "days",   # must match the template
                template = TRUE, concat_dim = "time")


gdalraster::mdim_info(vrt)
# {
#   "type": "group",
#   "driver": "VRT",
#   "name": "/",
#   "dimensions": [
#     {
#       "name": "lat",
#       "full_name": "/lat",
#       "size": 720,
#       "type": "HORIZONTAL_Y",
#       "direction": "NORTH",
#       "indexing_variable": "/lat"
#     },
#     {
#       "name": "lon",
#       "full_name": "/lon",
#       "size": 1440,
#       "type": "HORIZONTAL_X",
#       "direction": "EAST",
#       "indexing_variable": "/lon"
#     },
#     {
#       "name": "time",
#       "full_name": "/time",
#       "size": 16372,
#       "type": "TEMPORAL",
#       "indexing_variable": "/time"
#     },
#     {
#       "name": "zlev",
#       "full_name": "/zlev",
#       "size": 1,
#       "type": "VERTICAL",
#       "direction": "DOWN",
#       "indexing_variable": "/zlev"
#     }
#   ],
#   "arrays": {
#     "time": {
#       "full_name": "/time",
#       "datatype": "Float64",
#       "dimensions": [
#         "/time"
#       ],
#       "dimension_size": [
#         16372
#       ],
#       "attributes": {
#         "long_name": "Center time of the day"
#       }
#     },
#     "zlev": {
#       "full_name": "/zlev",
#       "datatype": "Float64",
#       "dimensions": [
#         "/zlev"
#       ],
#       "dimension_size": [
#         1
#       ],
#       "attributes": {
#         "actual_range": "0, 0",
#         "long_name": "Sea surface height",
#         "positive": "down"
#       }
#     },
#     "lat": {
#       "full_name": "/lat",
#       "datatype": "Float64",
#       "dimensions": [
#         "/lat"
#       ],
#       "dimension_size": [
#         720
#       ],
#       "attributes": {
#         "grids": "Uniform grid from -89.875 to 89.875 by 0.25",
#         "long_name": "Latitude"
#       }
#     },
#     "lon": {
#       "full_name": "/lon",
#       "datatype": "Float64",
#       "dimensions": [
#         "/lon"
#       ],
#       "dimension_size": [
#         1440
#       ],
#       "attributes": {
#         "grids": "Uniform grid from 0.125 to 359.875 by 0.25",
#         "long_name": "Longitude"
#       }
#     },
#     "anom": {
#       "full_name": "/anom",
#       "datatype": "Int16",
#       "dimensions": [
#         "/time",
#         "/zlev",
#         "/lat",
#         "/lon"
#       ],
#       "dimension_size": [
#         16372,
#         1,
#         720,
#         1440
#       ],
#       "block_size": [
#         1,
#         1,
#         720,
#         1440
#       ],
#       "attributes": {
#         "long_name": "Daily sea surface temperature anomalies",
#         "valid_max": 1200,
#         "valid_min": -1200
#       },
#       "unit": "Celsius",
#       "nodata_value": -999,
#       "offset": 0,
#       "scale": 0.0099999997764825821
#     },
#     "err": {
#       "full_name": "/err",
#       "datatype": "Int16",
#       "dimensions": [
#         "/time",
#         "/zlev",
#         "/lat",
#         "/lon"
#       ],
#       "dimension_size": [
#         16372,
#         1,
#         720,
#         1440
#       ],
#       "block_size": [
#         1,
#         1,
#         720,
#         1440
#       ],
#       "attributes": {
#         "long_name": "Estimated error standard deviation of analysed_sst",
#         "valid_max": 1000,
#         "valid_min": 0
#       },
#       "unit": "Celsius",
#       "nodata_value": -999,
#       "offset": 0,
#       "scale": 0.0099999997764825821
#     },
#     "ice": {
#       "full_name": "/ice",
#       "datatype": "Int16",
#       "dimensions": [
#         "/time",
#         "/zlev",
#         "/lat",
#         "/lon"
#       ],
#       "dimension_size": [
#         16372,
#         1,
#         720,
#         1440
#       ],
#       "block_size": [
#         1,
#         1,
#         720,
#         1440
#       ],
#       "attributes": {
#         "long_name": "Sea ice concentration",
#         "valid_max": 100,
#         "valid_min": 0
#       },
#       "unit": "%",
#       "nodata_value": -999,
#       "offset": 0,
#       "scale": 0.0099999997764825821
#     },
#     "sst": {
#       "full_name": "/sst",
#       "datatype": "Int16",
#       "dimensions": [
#         "/time",
#         "/zlev",
#         "/lat",
#         "/lon"
#       ],
#       "dimension_size": [
#         16372,
#         1,
#         720,
#         1440
#       ],
#       "block_size": [
#         1,
#         1,
#         720,
#         1440
#       ],
#       "attributes": {
#         "long_name": "Daily sea surface temperature",
#         "valid_max": 4500,
#         "valid_min": -300
#       },
#       "unit": "Celsius",
#       "nodata_value": -999,
#       "offset": 0,
#       "scale": 0.0099999997764825821
#     }
#   }
# }

get the latest one

gdalraster::plot_raster(x <- gdalraster::mdim_as_classic(vrt, "sst", 1L, 0L, view_expr = "[-1,0,:,:]"))