rspatial
FEATURE REQUEST: simplify need to set.crs for warper geolocation
This is a feature request, to somehow allow use of the warper with geolocation arrays to be more discoverable (?)
I'm not really sure what the UX would be, but perhaps the status of the raw swath data, understood by GDAL could be conveyed - and the crs of the geolocation arrays auto-applied at the terra level. (I know this is poorly formed but I think it's worth sharing for longer term exploration).
At any rate, we can open the data set and then project() it easily to target rasters, but we need to set.crs() to the geolocation arrays, but afaik we can't discover this situation from terra atm.
dsn <- 'NETCDF:"/vsicurl/https://dapds00.nci.org.au/thredds/fileServer/fj7/Copernicus/Sentinel-5/TROPOMI/L2__AER_AI/2023/2023-05/2023-05-27/S5P_OFFL_L2__AER_AI_20230527T213922_20230527T232052_29120_03_020500_20230529T112412.nc":/PRODUCT/aerosol_index_354_388'
library(terra)
#> terra 1.7.3
## we don't have georeferencing
r0 <- rast(dsn)
#> Error in R_nc4_open: No such file or directory
#> Warning: [rast] GDAL did not find an extent. Cells not equally spaced?
## but, GDAL knows about the geolocation arrays
writeLines(system(sprintf("gdalinfo %s | grep Geolocation -A 10", dsn), intern = TRUE))
#> Geolocation:
#> LINE_OFFSET=0
#> LINE_STEP=1
#> PIXEL_OFFSET=0
#> PIXEL_STEP=1
#> 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_BAND=1
#> X_DATASET=NETCDF:"/vsicurl/https://dapds00.nci.org.au/thredds/fileServer/fj7/Copernicus/Sentinel-5/TROPOMI/L2__AER_AI/2023/2023-05/2023-05-27/S5P_OFFL_L2__AER_AI_20230527T213922_20230527T232052_29120_03_020500_20230529T112412.nc":/PRODUCT/longitude
#> Y_BAND=1
#> Y_DATASET=NETCDF:"/vsicurl/https://dapds00.nci.org.au/thredds/fileServer/fj7/Copernicus/Sentinel-5/TROPOMI/L2__AER_AI/2023/2023-05/2023-05-27/S5P_OFFL_L2__AER_AI_20230527T213922_20230527T232052_29120_03_020500_20230529T112412.nc":/PRODUCT/latitude
#> Corner Coordinates:
## so we invoke the warper
try(project(r0, rast()))
#> Error : [project] input raster CRS not set
## but it needs a source CRS (the geolocation arrays have this)
set.crs(r0, "OGC:CRS84") ## (there's no geotransform so we fall back to the arrays, or force with -geoloc)
## now it works
plot(project(r0, rast())); maps::map(add = TRUE)
prj <- "+proj=laea +lat_0=56 +lon_0=180"
plot(project(r0, rast(ext(-1.3e7, .8e7, -.5e7, 1e7), ncols = 512, nrows = 512, crs = prj)))
lines(project(do.call(cbind, maps::map(plot = F)[1:2]), to = prj, from = "OGC:CRS84"))
#> Warning: [project] 1972 failed transformations
Created on 2023-05-30 by the reprex package (v2.0.1)
This also is kinda related to #107, because simply by nominating the arrays with VRT (if GDAL hasn't already discovered them) the general curvlinear case can be immediately warped to any target raster.
I think this implies a new property on SpatRasterSource to indicate if the GEOLOCATION domain is populated in the gdal metadata, (with potentially multiple sources in a SpatRaster being a situation that shouldn't be allowed). I'm not super clear on how a source with geoloc arrays might present as a SpatRaster but maybe replace the "extent not determined" indicator with something about this being a placeholder until after warp.
I should be done to cover other cases like GCPs, and actually expose the crs of these explicitly. I'll explore this 🙏