vapour_vrt(): geolocation is only working when 'bands' specified is length 1

[mdsumner]

Only when we have 1 band chosen does it work, when multiple or NULL (all bands) it doesn't get written into the GEOLOCATION metadata domain.

sds <- vapour::vapour_sds_names("CMIP_extract_forgdal.nc")

for (i in list(1, 1:2, NULL)) {
## 1 is intpp variable (x, y, time), 2,3 is lon(x), lat(y)
vrt1 <- vapour::vapour_vrt(sds[1], geolocation = sds[3:2], bands = i)
tx <- unlist(strsplit(vrt1, "\n"))
writeLines(tx[grep("GEOLOCATION", tx) + c(0:10)])
print(i)
}


Warning 1: dimension #2 (x) is not a Longitude/X dimension.
Warning 1: dimension #1 (y) is not a Latitude/Y dimension.
  <Metadata domain="GEOLOCATION">
    <MDI key="LINE_OFFSET">0</MDI>
    <MDI key="LINE_STEP">1</MDI>
    <MDI key="PIXEL_OFFSET">0</MDI>
    <MDI key="PIXEL_STEP">1</MDI>
    <MDI key="SRS">OGC:CRS84</MDI>
    <MDI key="X_BAND">1</MDI>
    <MDI key="X_DATASET">NETCDF:"CMIP_extract_forgdal.nc":lon</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="Y_DATASET">NETCDF:"CMIP_extract_forgdal.nc":lat</MDI>
  </Metadata>
[1] 1
Warning 1: dimension #2 (x) is not a Longitude/X dimension.
Warning 1: dimension #1 (y) is not a Latitude/Y dimension.
[1] 1 2
Warning 1: dimension #2 (x) is not a Longitude/X dimension.
Warning 1: dimension #1 (y) is not a Latitude/Y dimension.
NULL

CMIP_extract_forgdal.nc.zip

file created with

# data/CMIP6_data/ACCESS-ESM1-5/ssp126/CMIP6.ScenarioMIP.CSIRO.ACCESS-ESM1-5.ssp126.r1i1p1f1.Omon.intpp.gn.v20210318.2015-2024_2091-2100.nc
# ncks -v intpp,lon,lat -d time,60280.5 $f CMIP_extract_forgdal.nc


nc <- RNetCDF::open.nc("CMIP_extract_forgdal.nc", write = TRUE)
lon <- RNetCDF::var.get.nc(nc, "lon")
lon[lon > 180] <- lon[lon > 180] - 360

RNetCDF::var.def.nc(nc, "lon180", vartype = "NC_DOUBLE", c("x", "y", "time"))
RNetCDF::var.put.nc(nc, "lon180", data = array(lon, c(dim(lon), 1L)))
RNetCDF::close.nc(nc)

[mdsumner]

fixed the band thing

same deal as with #168, and now also works with terra:

x <- rast(vrt1)
set.crs(x, "OGC:CRS84")
plot(project(x, rast()))