Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Gdalwarp producing odd results around 0 degree longitude #1182

Open
simonrp84 opened this issue Jan 8, 2019 · 4 comments
Open

Gdalwarp producing odd results around 0 degree longitude #1182

simonrp84 opened this issue Jan 8, 2019 · 4 comments

Comments

@simonrp84
Copy link

Hello,

I am warping ECMWF grib data (full globe, 0.25 degree resolution) to
match the SEVIRI geostationary satellite projection. I use this command:

gdalwarp -overwrite -ot Float32 -t_srs "+proj=geos +a=6378169
+b=6356583.8 +lon_0=0.0 +h=35785831" -te -5570248.832537 5570244.832537
5567248.429179 -5567252.429179 -ts 3712 3712 -wo NUM_THREADS=6 -wo
SOURCE_EXTRA=1000 -q -r cubicspline --config CENTER_LONG 180
/data/ECMWF_T2_20181223T00_000.grib /data/cubicspline.tiff

The output projection looks good, and if I try this with the ECMWF
land-sea mask then I can see that it is correct for most of the view.
However, along a vertical centred on zero degrees longitude I see this
artefact: https://i.stack.imgur.com/tQn5F.png

Any ideas what is causing that, and how to resolve this bug? The artefact
appears for all resampling methods except nearest neighbour.

Thanks,
Simon

@rouault
Copy link
Member

rouault commented Jan 11, 2019

Can you point to a source dataset that can be used to reproduce the issue ?

@rouault rouault added the awaiting_feedback Awaiting feedback from reporter label Jan 11, 2019
@simonrp84
Copy link
Author

Apologies for the delayed reply. You can find an example dataset here:
https://github.com/simonrp84/temp_bugreport_data/blob/master/ECMWF_CAPE_20190121T12_000.grib

I tested this with:
gdalwarp -overwrite -ot Float32 -t_srs "+proj=geos +a=6378169 +b=6356583.8 +lon_0=0.0 +h=35785831" -te -5570248.832537 5570244.832537 5567248.429179 -5567252.429179 -ts 3712 3712 -wo NUM_THREADS=6 -wo SOURCE_EXTRA=1000 -q -r cubicspline --config CENTER_LONG 180 /data/ECMWF_CAPE_20190121T12_000.grib /data/cubicspline.tiff

I also tried several other CENTER_LONG values, and omitting the parameter entirely.

@simonrp84
Copy link
Author

Is there any update on this?

@rouault
Copy link
Member

rouault commented Feb 4, 2019

I managed to reproduce the issue. The fundamental problem is that the resampling kernel is not geography aware, so when interpolating around longitude = 0 or 360 in your source raster, it is has the edge of the raster and fallback to nearest resampling.
A workaround is to "shift" your raster in the more usual [-180,180] range

gdal_translate ECMWF_CAPE_20190121T12_000.grib left.vrt -a_ullr -360.125 90.125 -0.125 -90.125 -of VRT
gdalbuildvrt in.vrt left.vrt ECMWF_CAPE_20190121T12_000.grib
gdal_translate in.vrt in_shifted.vrt -projwin -180.125 90.125 180.125 -90.125 -of VRT
gdalwarp -overwrite -ot Float32 -t_srs "+proj=geos +a=6378169 +b=6356583.8 +lon_0=0.0 +h=35785831" -te -5570248.832537 5570244.832537 5567248.429179 -5567252.429179 -ts 3712 3712 -wo NUM_THREADS=6 -wo SOURCE_EXTRA=1000  -r cubicspline in_shifted.vrt  out.tif 

@rouault rouault removed the awaiting_feedback Awaiting feedback from reporter label Feb 4, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants