mdsumner

dsn <- "/vsicurl/https://data.source.coop/vizzuality/hfp-100/hfp_2021_100m_v1-2_cog.tif"
dsn <- basename(dsn)
info <- vapour::vapour_raster_info(dsn)

library(grout)
g <- grout(info$dimension, info$extent, blocksize = c(2048, 2048))

index <- tile_index(g)

mdsumner

This NetCDF has a spherical CRS longlat crs in "crs_wkt", but the NetCDF driver hardcodes WGS84 as the geolocation SRS:

gdalinfo /vsicurl/https://github.com/rgdal-dev/rasterwise/raw/master/data-raw/spherical.nc

...
Geolocation:
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  X
...

mdsumner

#dsn <- dsn::vrtcon(sds::gebco(), ovr = 4)
dsn <- "vrt:///vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif?ovr=4"

scaled.vrt <- "scaled.vrt"
system(glue::glue("gdal_translate {dsn} scaled.vrt -scale -ot Byte"))

## here I had to edit the pal.vrt to have "scaled.vrt" rather than "mem" in the source (...)
system(glue::glue("gdalattachpct.py pal scaled.vrt pal.vrt -of VRT"))

mdsumner

Note that the authorization mechanism requires GDAL 3.6 or later, find your version in terra with terra::gdal().

library(terra)
url <- "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20240318090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"

## make a GDAL dsn for an online file
dsn <- sprintf("/vsicurl/%s", url)

mdsumner

gdal_create -outsize 36000 17999 -ot Int8 -co SPARSE_OK=YES -a_srs EPSG:4326 -a_ullr 0 17999 0 36000 weird.tif
gdal_translate weird.tif cog.tif  -of COG
gdalinfo cog.tif | grep Overviews
#  Overviews: 18000x8999, 9000x4499, 4500x2249, 2250x1124, 1125x562, 562x281, 281x14

import xarray

mdsumner

iris$geom <- as.character(wk::as_wkt(wk::xy(1:150, 150:1)))
write.csv(iris, "iris.csv", row.names = F)
system("ogr2ogr iris.gdb  iris.csv  -oo GEOM_POSSIBLE_NAMES=geom iris -f OpenFileGDB -nlt POINT")

Warning 6: Normalized/laundered field name: 'Sepal.Length' to 'Sepal_Length'
Warning 6: Normalized/laundered field name: 'Sepal.Width' to 'Sepal_Width'

mdsumner

Interesting question! I think there's creative multiplication at the "just one level above xarray" layer because of that fundamental disconnect that occurs between storing every grid with coordinate arrays vs storing it with an offset/scale.

Not every grid is regular of course:

• some are truly curvilinear and some are truly rectilinear
• some are curvilinear or rectilinear by not storing the original regular grid in its native coordinate system but devolving that to points in longitude/latitude
• some are apparently rectilinear by virtue of small numeric differences (and sometimes just lack of explicit-intention)

There's also the disconnect of what the intention is, a regular grid can be specified by six(1) numbers xmin,xmax,ymin,ymax,ncol,nrow but so can the target of an operation, we can subset, subsample, reproject - all of which can be done included in one operation by specifying the output grid, xmin,xmax,ymin,ymax,ncol,nrow.

mdsumner

library(sf)
library(ggplot2)
library(gstat)
library(stars)
crs <- st_crs("EPSG:3857")

ex_data <- tibble::tribble(
    ~latitude,       ~longitude,   ~z,
    -39.0554861111111, 131.017636111111, 5.03,

mdsumner

GHRSST has a bit of a broken grid, the coordinates are degenerate rectilinear (i.e. completely described by the grid in xmin,xmax,ymin,ymax and the dimensions 36000x17999 but they store the longitude and latitude 1D coordinate arrays materalized as 32-bit floating point with resultant noise). It's not clear how the grid should be aligned exactly, and the metadata attributes list -180,180, -90,90 as the valid range but this can't be correct at 0.01 degree resolution for that size grid).

(ncdump -h listed below)

## replace this local file with the one at the new earthdata store commented out (but you need your earthdata creds, or Authorization Bearer token)

mdsumner

it's almost right, but stored in float32 with fuzz - degenerate rectilinear is a terrible way to store a regular grid ... is it regular? did they make a mistake? the coords and the metadata don't agree, it's impossible to be sure and that's normal in netcdf, a regular grid is six numbers (maybe eight) and a crs. sadly xarray has elevated this practice to world standard - I appreciate odc has helpers, but that's python - I R has helpers also of course, I like to keep as our in GDAL as much as possible so it's lang-agnostic