gdalwarp: transformation method GEOLOC_ARRAY from longitude 0,360

mdsumner commented 1 year ago

When using 'SRC_METHOD=GEOLOC_ARRAY' it seems that longitudes > 180 are not meaningfully wrapped.

Using two VRT files 'CMIP_pacific.vrt' and 'CMIP_atlantic.vrt' that reference GEOLOCATION arrays in the 'CMIP_extract_forgdal.nc' file, these are zipped into this link (and wget used in example below)

The pacific version references arrays lon,lat - longitude is in the range 0,360 (as per the original source).

The atlantic version references arrays lon180,lat - longitude is a modified version of the original 'lon' converted to -180,180.

Expected behavior and actual behavior.

I expect the geolocation array method to work with longitudes in the range 0,360 - but not sure this is a deficiency vs a feature request, or perhaps there's a workaround I am missing.

This figure shows the result of gdalwarp on the two VRT variants, with data missing from the western hemisphere (when the GEOLOCATION longitudes are > 360), but data complete when the GEOLOCATION longitudes are -180,180 range).

Steps to reproduce the problem.

Download and unzip this file containing

'CMIP_pacific.vrt'
'CMIP_atlantic.vrt'
'CMIP_extract_forgdal.nc'

wget https://github.com/hypertidy/vapour/files/9867052/geoloc_array_pacific_orientation.zip
unzip geoloc_array_pacific_orientation.zip

## create the pacific form (with geolocation longitudes in the range 0,360)
gdalwarp CMIP_pacific.vrt pacific.tif

## create the atlantic form (with geolocation longitudes in the range -180,180)
gdalwarp CMIP_atlantic.vrt atlantic.tif

Operating system

Ubuntu 20.04.5 LTS

GDAL version and provenance

GDAL 3.6.0dev-83983f521d-dirty, released 2022/10/25 (debug build)

mdsumner commented 1 year ago

It is in the clamping, this psTransform->bGeographicSRSWithMinus180Plus180LongRange (in several places) tests true for the Pacific case (longitudes 0,360) based on the IsGeographicSRS() test - the geolocation domain has 'SRS=OGC:CRS84'.

https://github.com/OSGeo/gdal/blob/83983f521dbe065452aa2418081678bf46e7f2ff/alg/gdalgeoloc.cpp#L302

So, a partial workaround is remove the SRS item from the geolocation metadata, but then gdalwarp provides output on -358.895, 361.2911 -90,90 with data in 0,360 - but if you specify the '-te -180 -90 180 90' then the data is still missing. I feel like this could be made to work, seamlessly handle data in 180,360 but I need to explore the overall logic some more before I can attempt a PR 🙏

mdsumner commented 10 months ago

see an in-memory workaround inspired by rasterio discussion: https://gist.github.com/mdsumner/46471f4b82d224fb755ae9bf96d069de

mdsumner commented 10 months ago

here's another example, with an extra complication for the nodata value (I just infill that and it seems to be fine - but presumably we can set the nodata for the DATAPOINTER workaround, but better to have the warper do the wrap) - I'm just glad to have a woarkound and will follow up for a PR deeper down.

https://gist.github.com/mdsumner/f9331f909d692bc9968251dd2b8bd678?permalink_comment_id=4757875#gistcomment-4757875

OSGeo / gdal