the siconc example

[mdsumner]

library(vapour)
file <- "/perm_storage/home/mdsumner/siconc_SIday_MPI-ESM1-2-HR_ssp126_r1i1p1f1_gn_20150101-20191231.nc"
dsn <- sprintf("vrt://NETCDF:%s:siconc?bands=150", file)
opts <- c("-to", "SRC_METHOD=NO_GEOTRANSFORM")
library(ximage)
ximage(x <- gdal_raster_data(dsn, options = opts), col = hcl.colors(12))
ximage(volcano, col = hcl.colors(12), breaks = 1:13)

library(palr)
image_pal(volcano[1:4, 1:5])

ll <- sprintf("NETCDF:%s:%s", file, c("longitude", "latitude"))
info <- vapour_raster_info(dsn)
lon0 <- gdal_raster_data(ll[1L], options = opts)
range(lon0[[1]])

wraplon <- function(x) ((x - -180) %% 360) - 180
lon <- matrix(wraplon(lon0[[1]]), attr(lon0, "dimension")[2], byrow = TRUE)
range(lon)

mk_geoloc <- function(x, bands = 1L) {
  vapour_vrt(x, geolocation = c(dsn::mem(t(lon)), ll[2]), bands = bands)
}

## resolved grid
d1 <- gdal_raster_data(mk_geoloc(dsn, bands = 1L), target_ext = c(-1, 1, -1, 1) * 5e6, target_crs = "+proj=laea +lon_0=83 +lat_0=-90")
ximage::ximage(d1, asp = 1)
range(d1[[1]], na.rm = T)

[mdsumner]

Note that I'm using a very recent GDAL, so it's probably not possible to share this for reuse yet.

See how the long/lat coords are concentrated in the north in this curvlinear grid.

library(vapour)
file <- "~/siconc_SIday_MPI-ESM1-2-HR_ssp126_r1i1p1f1_gn_20150101-20191231.nc"
ll <- sprintf("NETCDF:%s:%s", file, c("longitude", "latitude"))

opts <- c("-to", "SRC_METHOD=NO_GEOTRANSFORM")

lon <- gdal_raster_data(sprintf(ll[1], file), options = opts)
lat <- gdal_raster_data(sprintf(ll[2], file), options = opts)

plot(lon[[1]], lat[[1]], pch =  ".")
maps::map("world2", add = TRUE, col = "hotpink", )

This is a triipolar grid - they put the nodes over land so that the grid is "nice" for the ocean and we don't care about those bits inside continental land area

.

## now the data, lets look at in raw form
dsn <- sprintf("vrt://NETCDF:%s:siconc?bands=1", file)

library(ximage)
ximage(x <- gdal_raster_data(dsn, options = opts))

totally understandable, but note that the longitude latitude coordinates while on a grid are actually not regular at all

and they are pretty odd with a strange edge

 terra::plot(terra::rast(ll[1]))
 terra::plot(terra::rast(ll[2]))

Now it gets a bit hectic, we can use the GDAL warper but we have to do this weird trick to wrap longitude (I want GDAL to do this internally ... and this problem manifests in many different ways downstream in R and Python so it's a valuable thing to get fixed).

## wrap longitude
wraplon <- function(x) ((x - -180) %% 360) - 180

## helper to overwrite the geolocation array metadata in a source
mk_geoloc <- function(x, bands = 1L) {
  vapour_vrt(x, geolocation = c(dsn::mem(t(lon)), ll[2]), bands = bands)
}

info <- vapour_raster_info(dsn)

## I'm going to read this copy, modify it in memory and then have GDAL reference 
## that updated version for the geolocation array
lon0 <- gdal_raster_data(ll[1L], options = opts)
range(lon0[[1]])
#[1]   0.00591556 359.99893487

## modified copy
lon <- matrix(wraplon(lon0[[1]]), attr(lon0, "dimension")[2], byrow = TRUE)
range(lon)
#[1] -179.9941  179.9989

par(bg = "darkgrey")
## resolved grid
d1 <- gdal_raster_data(mk_geoloc(dsn, bands = 1L),  target_crs = "+proj=laea +lat_0=90")
ximage::ximage(d1, asp = 1)
range(d1[[1]], na.rm = T)
#[1]   0 100

Here's our nearly-respectable gridded version of the grid (for northern use at least, but I'm not sure why it's truncated like that, will still explore ), we could make a longlat grid, or transverse mercator isn't bad for both poles ... but it really depends on the required use, and we can define a nice grid exactly and check that it's appopriate - or it might be better still just to treat the points as points and use fast indexed point lookup for them ...

This whole narrative has been beyond my grasp for a while, but it ties a lot of stuff together and I'm really interested in it.