Rasterio: Nodata mask doesn't match fill values with non-integer valued windows

[drewbo]

Expected behavior and actual behavior.

When doing a windowed read with masking, the mask should match places where the data is equal to the fill value. This occurs when the window is integer valued but not otherwise. (the existence of the window methods round_shape and round_offsets make me think this is a known issue/expected behavior but I couldn't find any comments related to it)

Steps to reproduce the problem

import rasterio
from rasterio.windows import Window

window = Window(col_off=21195.542639606836, row_off=10993.82477100648, width=37.73706566452529, height=46.59517549596967)

with rasterio.open('s3://mytif.tif') as src:
    data = src.read(window=window, masked=True)
    print('original window -- ', f'mask pixels: {data.mask.sum()},', f'fill value pixels: {(data.data == data.fill_value).sum()}')

    data = src.read(window=window.round_shape(), masked=True)
    print('round shape -- ', f'mask pixels: {data.mask.sum()},', f'fill value pixels: {(data.data == data.fill_value).sum()}')

    data = src.read(window=window.round_offsets(), masked=True)
    print('round offsets -- ', f'mask pixels: {data.mask.sum()},', f'fill value pixels: {(data.data == data.fill_value).sum()}')

    data = src.read(window=window.round_shape().round_offsets(), masked=True)
    print('round all -- ', f'mask pixels: {data.mask.sum()},', f'fill value pixels: {(data.data == data.fill_value).sum()}')

>>> original window --  mask pixels: 787, fill value pixels: 794
>>> round shape --  mask pixels: 787, fill value pixels: 794
>>> round offsets --  mask pixels: 821, fill value pixels: 794
>>> round all --  mask pixels: 794, fill value pixels: 794

The tif I discovered this on is 13G (hence the windowed read) but I'm happy to create a more minimal reproduction if this seems like a valid bug.

Operating system

For example: Mac OS X 10.15.4.

Rasterio version and provenance

rasterio==1.1.7 from pip 20.2.2

[sgillies]

@drewbo can you tell me more about your dataset? If I understand correctly, it has no mask, but has a nodata value? I should be able to reproduce this with rasterio's RGB.byte.tif file.

[drewbo]

@sgillies yes, we didn't create a mask directly but the data does have a nodata value. I can also recreate on the test file

$ rio insp tests/data/RGB.byte.tif

from rasterio.windows import Window
window = Window(20.2, 130.8, 240.3, 450.7)
data = src.read(window=window, masked=True)
(data.mask.sum(), (data.data == data.fill_value).sum())
>>> (72381, 72705)

[sgillies]

@drewbo confirmed with rasterio 1.1.7 and GDAL 3.1.3, thanks! This report and the examples are super helpful. My first hypothesis was that the problem would be around the edges of the window (off-by-one error always being a prime suspect), but in plotting ((data.data[0] == data.fill_value) == data.mask[0]) it seems that is not the case. Yellow is True, Purple is False in the figure below. The discrepancy is around the edge of 0 values in the dataset.

One major difference between reading using integer and float windows is that resampling has to be done in the latter case. In #2023 I've added tests and found that the problem occurs with nearest neighbor resampling (the default), but not with bilinear or bicubic resampling (where I count 70735 masked pixels and 70735 pixels with the value == to the nodata/fill value).

I do not know why the mask doesn't expand as the 0 values in dataset expand under nearest neighbor resampling, while the mask contracts as the 0 values contract under the other resampling algorithms. I think we need to ask upstream. Meanwhile, I think we'll have to work around this one.

[sgillies]

BTW, this issue persists in GDAL 3.1.4RC2.

[sgillies]

https://github.com/OSGeo/gdal/issues/3102 fixes the GDAL bug but won't be in a release until 3.2.0. I've made a patch for GDAL 3.1.4 that also applies to 2.4.4 and am testing it at https://github.com/rasterio/rasterio-wheels/pull/60. If this works, I'll make a 1.1.8 post-release with the GDAL patch.