gdalwarp not accurate

[mwooten3]

I'm trying to aggregate data up from 0.05 to a coarser 0.25 resolution and am having trouble getting gdalwarp to work correctly. From comparing the original input data with the aggregated data and manually calculating the averages, I can see that the result I'm getting is not what I would expect. Many (not all) aggregated pixel values are just wrong (FYI, I have specifically been looking where the water meets land on the horn of Africa). Also, it's my understanding that gdal is supposed to ignore NoData pixels in the average calculation, but in some areas, the 25 pixels underneath the aggregated pixels were all NoData, yet the aggregated pixel reported a number, when it should have been NoData.

What's strange is that it works for another dataset I'm trying to run the same process on, and I can't figure out what's different about the input data to cause this. I'm going to attach both the inputs and outputs of each dataset.

Doesn't work:

gdalwarp -tr .25 .25 -r average -srcnodata -15000 -dstnodata -15000 /path/to/CHIRPS_201512.tif /path/to/CHIRPS_201512_ave.tif

Works:

gdalwarp -tr .25 .25 -r average -srcnodata -15000 -dstnodata -15000 /path/to/NDVI_201512.tif /path/to/NDVI_201512_ave.tif

If anyone has any clue as to what is going on and can help solve it, I have a deadline coming up soon so I need solutions ASAP.

Thank you.

[mwooten3]

The files are too big to attach here, so I've uploaded them to dropbox. You can find the four files here: ​https://www.dropbox.com/sh/dk0c41638jg8fir/AAD9z9IXeEfJ2PyEUVnjG4aya?dl=0

[Even Rouault]

I had a look and couldn't really find a difference between your "doesn't work" and "works" cases. In both cases, from what I could see, if at least one pixel value in a 5x5 source block was not no-data, a not no-data pixel value was output. And if all pixel values in the 5x5 source blocks where no-data, a no-data pixel value was output. But I may have missed the pixels you wre looking to.

[mwooten3]

Replying to rouault:

I had a look and couldn't really find a difference between your "doesn't work" and "works" cases. In both cases, from what I could see, if at least one pixel value in a 5x5 source block was not no-data, a not no-data pixel value was output. And if all pixel values in the 5x5 source blocks where no-data, a no-data pixel value was output. But I may have missed the pixels you wre looking to.

Really? It seems like everywhere I look (I have been looking at coastal areas of Horn of Africa, Florida, the Carribbean, Australia, Japan, but my study area is Africa) and there are several pixels that should be NoData but aren't. Additionally, in areas where one of the original pixels is X, the aggregated average should also be X, but it's not. I'm attaching a screenshot of the Horn of Africa and labeling the pixels (black for the smaller orginal pixels, white for the what the larger aggregated pixels turned out to be, and red for what the aggregate pixels *should* be). I also attached a map (ex2) with a rough extent rectangle (drew it in paint so it's not exact). Let me know what you think.

[mwooten3]

pixels screenshot- zoomed in and labeled

[mwooten3]

zoomed out to see rough extent of first screenshot

[Even Rouault]

OK I now also see the issue on the zone you spotted. I believe it might be linked to a rounding issue in the way source pixel coordinates are computed. If you're interested, I can further investigate to fix that : ​http://www.spatialys.com/en/services/

[mwooten3]

Replying to rouault:

OK I now also see the issue on the zone you spotted. I believe it might be linked to a rounding issue in the way source pixel coordinates are computed. If you're interested, I can further investigate to fix that : ​http://www.spatialys.com/en/services/

That would be great. I'm guessing your services are not free?

[Even Rouault]

You can email me at "contact at spatialys.com" for further details.

[mwooten3]

Replying to rouault:

You can email me at "contact at spatialys.com" for further details.

Done, thanks.

[Even Rouault]

Possible workarounds: As no reprojection is actually involved, if you use GDAL 2.0, you could use "gdal_translate -r average -outsize 20% 20% in.tif out.tif". And with lower versions, you could workaroud by adding overviews "gdaladdo -ro -r average in.tif 5" followed by "gdal_translate -outsize 20% 20% in.tif out.tif". In both cases, you'll need to properly set the nodata value first with a prior gdal_translate -a_nodata.

[ctl101]

I note very similar problems using -r average. I wondered when there might be a fix for this in gdalwarp? I note that -r mode seems to produce particularly spurious output although fortunately -r near seems unaffected.

Might the usage option to 'ignore nodata' in resampling calculations be implemented (as per the arcgis functionality in similar tools) in gdalwarp - particularly for -r mode? A specific option to avoid the production of nodata cells in output where -r mode finds two or more values occur equally often would also be extremely useful (perhaps with a rule option such as selecting the lower or higher value in such instances).

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.