OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.94k stars 2.57k forks source link

netCDF Split&Swap #5497

Open zakjan opened 2 years ago

zakjan commented 2 years ago

Expected behavior and actual behavior.

GRIB Split&Swap mode (#4526) would be useful for netCDF driver as well. What approach would be suitable to generalize the feature? Copy it to netCDF? Or expose as a separate gdal_translate/gdalwarp flag, so that it is supported by all drivers?

Steps to reproduce the problem.

Download source file sst_anomaly.nc.zip (comes from CMEMS).

$ gdalinfo NETCDF:sst_anomaly.nc:sst_anomaly
...
Size is 1440, 720
Origin = (0.000000000000000,90.000000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
...
Corner Coordinates:
Upper Left  (   0.0000000,  90.0000000)
Lower Left  (   0.0000000, -90.0000000)
Upper Right (     360.000,      90.000)
Lower Right (     360.000,     -90.000)
Center      ( 180.0000000,   0.0000000)

gdal_translate it to GeoTIFF. Expected result is in -180..180 longitude range, same as GRIB. Actual result is in 0..360 longitude range.

Operating system

macOS 12.1 M1

GDAL version and provenance

3.4.2 from brew

rouault commented 2 years ago

Or expose as a separate gdal_translate/gdalwarp flag

That should reasonably fit in gdal_translate architecture (controlled by an option), by using underneath a VRT assembling 2 sources.. gdalwarp -t_srs "+proj=longlat +datum=WGS84 +lon_wrap=0" could perhaps do the job too.

zakjan commented 2 years ago

Yes, I use +lon_wrap currently, but it requires a separate call to gdalwarp with a non-standard projection string and a temporary file. I like how the new GRIB Split&Swap allows to use a single gdal_translate command. Is there support for this transformation in VRT already?

rouault commented 2 years ago

Is there support for this transformation in VRT already?

no, I gave above the hint how that could possibly be implemented, but work needed

mmomtchev commented 2 years ago

This old stackoverflow by @warmerdam contains the right command: https://gis.stackexchange.com/questions/37790/reprojecting-raster-from-0-360-to-180-180-with-cutting-180-meridian-using-gdalw

gdalwarp -s_srs "+proj=longlat +ellps=WGS84" -t_srs WGS84 ~/0_360.tif 180.tif  -wo SOURCE_EXTRA=1000 \
           --config CENTER_LONG 0

GRIBv2 was a very particular case where the format itself specified that 0-360° longitudes are to be used

zakjan commented 2 years ago

There is my comment as well, CENTER_LONG didn't seem to work well. https://gis.stackexchange.com/questions/37790/reprojecting-raster-from-0-360-to-180-180-with-cutting-180-meridian-using-gdalw#comment665452_223966

gdalwarp -s_srs '+proj=longlat +datum=WGS84 +lon_wrap=180' -t_srs '+proj=longlat +datum=WGS84'

I get it that if we were to support the feature outside of gdalwarp, it should be a VRT transformation. Thanks, I'll take a look at it.

mdsumner commented 1 month ago

You could do

## set this so output has CRS set (or use 'vrt://...&a_srs=EPSG:4326') 
 export GDAL_NETCDF_ASSUME_LONGLAT=YES

gdalwarp -te -180 -90 180 90 \ 
 "vrt:///vsizip//vsicurl/https://github.com/OSGeo/gdal/files/8337193/sst_anomaly.nc.zip/sst_anomaly.nc?a_ullr=-360,90,0,-90" \
 /vsizip//vsicurl/https://github.com/OSGeo/gdal/files/8337193/sst_anomaly.nc.zip/sst_anomaly.nc atlantic.tif

If you only do one input, without adding the modified one you only get the east hemisphere, because the warper doesn't wrap (but it could I think). Potentially related to #6582