Rasterio: rio warp paints a wide black band on part of the output

[rasmuse]

Hi, and thanks for great software.

I'm running Ubuntu 16.04. rasterio is installed with conda create -n myvirtualenv python=2 and then conda install rasterio. So it's the conda-packaged GDAL behind it.

$ rio --version
0.36.0
$ gdalinfo --version
GDAL 2.1.1, released 2016/07/07

I am trying to use rio warp to resample (shrink) a GeoTIFF raster, but it's producing a strange artifact.

I'm using a natural earth raster (md5 sum just below).

$ rio --version
0.36.0
$ rio warp --res 1 NE1_50M_SR_W.tif warped.tif
$ md5sum NE1_50M_SR_W.tif warped.tif 
0a436f7b943aef85977017d9fff02f9d  NE1_50M_SR_W.tif
b25a314e9fa573661edb4f71078d9f5b  warped.tif

Screenshots:

Source:

Result:

The output file is (as expected) 360 x 180 px large. The 18 rightmost columns (= 18 rightmost degrees of longitude) are painted black.

$ rio info NE1_50M_SR_W.tif 
{"count": 3, "crs": "EPSG:4326", "colorinterp": ["red", "green", "blue"], "interleave": "pixel", "dtype": "uint8", "driver": "GTiff", "transform": [0.03333333333333, 0.0, -179.99999999999997, 0.0, -0.03333333333333, 90.0], "lnglat": [-1.7962520360015333e-11, 8.995471034722868e-12], "height": 5400, "width": 10800, "shape": [5400, 10800], "tiled": false, "res": [0.03333333333333, 0.03333333333333], "nodata": null, "bounds": [-179.99999999999997, -89.99999999998201, 179.99999999996405, 90.0]}
$ rio info warped.tif 
{"count": 3, "crs": "EPSG:4326", "colorinterp": ["red", "green", "blue"], "interleave": "pixel", "dtype": "uint8", "driver": "GTiff", "transform": [1.0, 0.0, -179.99999999999997, 0.0, -1.0, 90.0], "lnglat": [2.842170943040401e-14, 0.0], "height": 180, "width": 360, "shape": [180, 360], "tiled": false, "res": [1.0, 1.0], "nodata": null, "bounds": [-179.99999999999997, -90.0, 180.00000000000003, 90.0]}
$ rio sample warped.tif
[161, 0]
[162, 0]
[116, 168, 202]
[0, 0, 0]

The image looks the same (but with higher resolution) if I warp at higher destination resolution.

$ rio warp --res .1 NE1_50M_SR_W.tif warped-fine.tif
$ rio sample warped-fine.tif
[161.9, 0]
[162, 0]
[115, 165, 201]
[0, 0, 0]

So here it's still the 18 rightmost degrees of longitude that fail.

Similar result with another input file

Another Natural Earth dataset (medium size) gives similar results, but here it's about the rightmost 9 degrees longitude that are black.

$ rio info HYP_LR_SR_OB_DR.tif 
{"count": 3, "crs": "EPSG:4326", "colorinterp": ["red", "green", "blue"], "interleave": "pixel", "dtype": "uint8", "driver": "GTiff", "transform": [0.02222222222222, 0.0, -180.0, 0.0, -0.02222222222222, 90.0], "lnglat": [-1.801936377887614e-11, 9.00968188943807e-12], "height": 8100, "width": 16200, "shape": [8100, 16200], "tiled": false, "res": [0.02222222222222, 0.02222222222222], "nodata": null, "bounds": [-180.0, -89.99999999998198, 179.99999999996396, 90.0]}
$ rio warp --res 1 HYP_LR_SR_OB_DR.tif warped2.tif
$ rio info warped2.tif 
{"count": 3, "crs": "EPSG:4326", "colorinterp": ["red", "green", "blue"], "interleave": "pixel", "dtype": "uint8", "driver": "GTiff", "transform": [1.0, 0.0, -180.0, 0.0, -1.0, 90.0], "lnglat": [0.0, 0.0], "height": 180, "width": 360, "shape": [180, 360], "tiled": false, "res": [1.0, 1.0], "nodata": null, "bounds": [-180.0, -90.0, 180.0, 90.0]}
$ md5sum HYP_LR_SR_OB_DR.tif warped2.tif 
6bd010d9895d134a7c23570886621289  HYP_LR_SR_OB_DR.tif
cf1c189dd1f180e83aedad2b23462cb4  warped2.tif
$ rio sample warped2.tif 
[170, 0]
[171, 0]
[109, 168, 212]
[0, 0, 0]
$ rio warp --res 0.1 HYP_LR_SR_OB_DR.tif warped2-fine.tif
$ rio info warped2-fine.tif --shape
1800 3600
$ rio sample warped2-fine.tif 
[170.9, 0]
[171, 0]
[102, 166, 207]
[0, 0, 0]

Hope this helps!

[rasmuse]

Oh, and I should add that gdalwarp works as expected:

$ gdalwarp -tr 1 1 NE1_50M_SR_W.tif warped_gdal.tif
Creating output file that is 360P x 180L.
Processing input file NE1_50M_SR_W.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.
$ md5sum warped_gdal.tif 
46e9ca6f5e9c4a2ecb16938dddf060a3  warped_gdal.tif
$ rio info warped_gdal.tif 
{"count": 3, "crs": "EPSG:4326", "colorinterp": ["red", "green", "blue"], "interleave": "pixel", "dtype": "uint8", "driver": "GTiff", "transform": [1.0, 0.0, -179.99999999999997, 0.0, -1.0, 90.0], "lnglat": [2.842170943040401e-14, 0.0], "height": 180, "width": 360, "shape": [180, 360], "tiled": false, "res": [1.0, 1.0], "nodata": null, "bounds": [-179.99999999999997, -90.0, 180.00000000000003, 90.0]}

[rasmuse]

Sorry 'bout the mess. Did not intend to close the issue :)

[brendan-ward]

@rasmuse this should have been resolved in #581, but we were specifically testing reprojection to spherical mercator. Looks like we might be bumping into a case we didn't cover.

[brendan-ward]

I can reproduce on master. Definitely a bug.

[perrygeo]

I can reproduce as well, rasterio@master and gdal 2.1.1 on OSX. But I get a slightly different result - missing _4_ columns of pixels on the eastern edge.

Clearly something is awry. We've seen something similar to this bug before: https://github.com/mapbox/rasterio/issues/541 so I'll revisit that issue and see if we someone introduced a regression in recent versions.

[rasmuse]

@perrygeo Interesting. Exactly the same input file I suppose?

[perrygeo]

Yeah, same MD5. Not sure exactly what's happening but I'll be digging in shortly. Thanks for the great bug report, btw!

[perrygeo]

So I was able to reproduce this with our world.byte.tif test fixture (failing test here at warp-loss-redux branch)

rio warp --res 100000 tests/data/world.byte.tif /tmp/test_rio.tif --dst-crs EPSG:3857

19 columns of data lost.

If I run it with --res 10000, no such loss occurs

gdalwarp (v 2.1.1) appears to suffer from a similar problem (in this case it's missing 20 columns compared to rio warp's 19)

gdalwarp -t_srs EPSG:3857 -tr 100000 100000 tests/data/world.byte.tif /tmp/test_gdal.tif

Conclusions Nothing solid at this point but we can say that

• the bug is affected by resolution, particularly high -> low resolution warping
• the bug exists in both gdalwarp and rio warp indicating that the bug _might_ reside in the core GDAL warping function.

[sgillies]

@perrygeo I'm going to check with @rouault on this one and report back.

[perrygeo]

Per chat, gdalwarp's -tap option can usually solve these things

[perrygeo]

@rasmuse finally picking this up again, thanks for your patience 🙇

So it looks like gdalwarp's -tap option does the trick:

gdalwarp -tap -t_srs EPSG:3857 -tr 100000 100000 tests/data/world.byte.tif /tmp/test_gdal_tap.tif

I've implemented a similar feature, --target-aligned-pixels, in https://github.com/mapbox/rasterio/pull/941

rio warp --target-aligned-pixels --no-check-invert-proj  --dst-crs EPSG:3857 \
    --res 100000 tests/data/world.byte.tif /tmp/test_rio_tap.tif

Note that we need to include --no-check-invert-proj; with the default --check-invert-proj we still get strange results. This time there are blank strips on the left and right. This is due, I think, to over-eager application of projection math checks, something we'll address in https://github.com/mapbox/rasterio/issues/940

rio warp --target-aligned-pixels --dst-crs EPSG:3857 \
    --res 100000 tests/data/world.byte.tif /tmp/test_rio_tap_check.tif

[rasmuse]

Thank you for getting back!

It is great that rasterio can now produce the expected result. But at the same time, it is not entirely satisfying that two a seemingly unrelated flags would be needed to make such a basic operation work. I mean, at least in the original bug report I was not even reprojecting anything. (But perhaps it does take a reprojection roundtrip even if the source CRS equals the destination CRS?)

Edit:
One more thing. I now see that #940 might fix the choice of --no-check-invert-proj. But this would still require me to explicitly say --target-aligned-pixels in order to safely make a simple resampling operation?

[perrygeo]

I'd say that, even in the case of straight resampling, a --target-aligned-pixel option is directly related and necessary. Agreed that --no-check-invert-proj should be unnecessary, hence #940.

The root of the problem is that, in both gdalwarp and rio warp cases, we're using the warping machinery to do resampling. There is no projection roundtripping but the warping functions are used under the hood and all of the baggage comes with it. Ultimately we should have a resampling command that just does straight up/down sampling with numpy without touching reproject. See https://github.com/mapbox/rasterio/issues/755 for a related discussion.

[rasmuse]

Yes, I see now that the check-invert flag could be fixed with #940. You were just faster to respond than I was to edit my first question. Thanks anyway for the clarification!

Thanks also for mentioning #755. I had not seen that. And I would be very thankful for a separate resampling machinery.

[vodp]

I meet a similar bug, by calling reproject directly from Python API. The black band is not on the right edge, but at the top. As pointed out by @perrygeo the bug happens as I warp a high res img to lower res.

The correct one should be like

As I noticed the fixed bugs in rio command line does not apply for internal reproject function.

I am using rasterio version 1.0a6.