vapour_warp_raster: problem with GEOLOC to projected grid

mdsumner commented 2 years ago

(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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

hypertidy / vapour

vapour_warp_raster: problem with GEOLOC to projected grid #168