vapour_warp_raster: problem with GEOLOC to projected grid

[mdsumner]

(but does do atlantic and pacific, see #167 for anotherissue creating these VRT)

v <- vapour_warp_raster("CMIP_atlantic.vrt", extent = c(-1, 1, -1, 1) * 1e7, dimension  = c(512, 512), projection = "+proj=laea")
ERROR 1: Too many points (528 out of 529) failed to transform, unable to compute output bounds.
Warning 1: Unable to compute source region for output window 0,0,512,512, skipping.

range(v[[1]])
# [1] 0 0 

geoloc_array_pacific_orientation.zip

[mdsumner]

the problem is really that there is a CRS, but we think there is not

so, the logic that checks this should look for GCP and other geolocation sources that the warper might use ...

for workaround, just set the 'source_projection = "OGC:CRS84"'

[mdsumner]

this also a problem in the general case, but here it's clearly my fault

dsn <- "NETCDF:\"https://thredds.aodn.org.au/thredds/dodsC/IMOS/SRS/Surface-Waves/SAR_Wind/DELAYED/SENTINEL-1A/2018/04/14/IMOS_SRS-Surface-Waves_M_20180414_Coastal-Wind-Sentinel-1A_FV01_DM00-021467-024F96-0798.nc\":WSPD"
gdal_raster_data(dsn, target_crs = "+proj=laea +lon_0=145 +lat_0=-42", target_dim = c(256, 256))
# Warning message:
#   In warp_general_cpp(dsn, target_crs, as.numeric(target_ext), as.integer(target_dim),  :
#                         no source crs, target crs is ignored
                      

[mdsumner]

here

https://github.com/hypertidy/vapour/blob/main/inst/include/gdalwarpgeneral/gdalwarpgeneral.h#L127

and

https://github.com/hypertidy/vapour/blob/main/inst/include/gdalwarpmem/gdalwarpmem.h#L90

we need to allow it if there are geolocation arrays (or just allow anyway?)

[mdsumner]

I just turned off the stop when no source_crs in gdalwarpgeneral, seems fine (but in GDAL 3.7 ... might be a version thing)

[mdsumner]

Here is R code to do it all from scratch, with osgeo.gdal

library(reticulate)
gdal <- import("osgeo.gdal")
gdal$UseExceptions()
cmip_dsn <- "/vsizip//vsicurl/https://github.com/hypertidy/vapour/files/9867052/geoloc_array_pacific_orientation.zip"

gdal$ReadDirRecursive(cmip_dsn)
#> [1] "CMIP_extract_forgdal.nc" "CMIP_pacific.vrt"       
#> [3] "CMIP_atlantic.vrt"
ds0 <- gdal$Open(file.path(cmip_dsn, "CMIP_extract_forgdal.nc"))
sds = ds0$GetMetadata_List("SUBDATASETS")
sds <- sds[seq(1, length(sds), by = 2)]

sds <- gsub("^SUBDATASET_.*_NAME=", "", sds)

md <- list("LINE_OFFSET" = "0",
        "LINE_STEP" = "1",
        "PIXEL_OFFSET" = "0",
        "PIXEL_STEP" = "1",
        "X_DATASET" = grep(":lon180$", sds, value = T),
        "X_BAND" = "1",
        "Y_DATASET" = grep(":lat$", sds, value = T),
        "Y_BAND" = "1"
)

## now we want the actual dataset
ds <- gdal$Open(grep(":intpp$", sds, value = T))
ds$SetMetadata(md, "GEOLOCATION")
#> [1] 0

gdal$Warp("/vsimem/geo.tif", ds)
#> <osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7f68cbd129f0> >

library(vapour)
library(ximage)
ximage(gdal_raster_data("/vsimem/geo.tif", target_dim = c(256, 0)))

^{Created on 2023-09-09 with reprex v2.0.2}

now, need to modify the actual lon sds without hacking it upfront