[gdal-dev] gdalwarp across the antimeridian with a source extent >180E

[Even Rouault]

Michael,

https://github.com/OSGeo/gdal/pull/10108 will fix it.

You can also workaround the issue by overriding the extent of the source 
dataset to be exactly -179.995,89.995,180.005,-89.995 with

  "vrt://NETCDF:20240102090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326&a_ullr=-179.995,89.995,180.005,-89.995"

The issue was twofolds:

- the dataset uses lon/lat arrays with single precision floating-point 
numbers. Hence -179.99 gets actually read a a read as 
-179.99000549316406 as a double-precision number, and thus the maxx - 
minxx difference was slightly over 360 degrees. I've added a heuristics 
to try to round -179.99-single-precision as -179.99-double-precision

- which defeated a heuristics in the warper to identify the center 
longitude of a dataset registered in a geographic CRS. This center 
longitude is just (min_lon + max_lon) / 2 if the dataset doesn't cover 
more than 360 degrees of longitude. This value, when computed, is passed 
as a hint OGRCoordinateTransformation so that it can post-correct 
longitudes to apply a +/- 360 degree offset, to be in the range of the 
source dataset. I've relaxed the sanity check to allow slighly more than 
360 degrees

Even


Le 01/06/2024 à 08:19, Michael Sumner via gdal-dev a écrit :
> This data source has an odd georeferencing, it's a 36000x17999 raster 
> in the extent -179.995 -89.995 180.005 89.995.
>
> vrt://NETCDF:/vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20240101090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326
>
> (Why is it 17999: well there's one missing row - you'd expect 18000 at 
> 0.01 resolution and that manifests as half a cell short at the top and 
> bottom).
>
> But also the x registration is half a cell to the east from what you 
> would expect.
>
> It's correctly registered in that frame though, if we shift it west by 
> the half cell I've verified visually that it "looks wrong" compared to 
> coastlines etc.
>
> The point of this report is, gdalwarp misses the western window of 
> data for a request across the antimeridian (there's no data in the 
> output east of 180).
>
> gdalwarp $dsn -te 173.00 -19.85 185.00 -15.00 out1.tif -co COMPRESS=LZW
>
> If we "fix" the extent to be more longitudinally pleasing, we get the 
> data east of the antimeridian, the map is properly filled either side 
> of the antimeridian.
>
>
> vrt://NETCDF:/vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20240101090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326&a_ullr=-180,89.995,180,-89.995 
> <https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20240101090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc:analysed_sst?a_srs=EPSG:4326&a_ullr=-180,89.995,180,-89.995>
>
> gdalwarp $dsnfix -te 173.00 -19.85 185.00 -15.00 out2.tif -co COMPRESS=LZW
>
>
> Which is good, but as mentioned above the workaround means we are now 
> half a cell out. It works fine for projected windows over the 
> antimeridian in either situation.
>
> The compressed file sizes reflect the missing data in the second one, 
> these are 625K and 1.1Mb.
>
> I'd like the warper to be able to correctly fill the data from the 
> western window in the original case, we have determined that we can't 
> feasibly "fix" the extent.
>
> Or maybe something else I'm missing? I don't *think* SAMPLE_GRID or 
> SOURCE_EXTRA helps here. I know that I could craft a meridian-crossing 
> combination in VRT but I'd rather like this workflow to work as-is.  
> (Chasing down the weird grid referencing is another project, but it's 
> a really massive dataset so I'm pretty unconfident of that occurring).
>
> Cheers, Mike
>
> -- 
> Michael Sumner
> Research Software Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsumner at gmail.com
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240602/d3435d8b/attachment.htm>