Rasterio: Bilinear resampling creates invalid mask values

[vincentsarago]

Ref: https://github.com/cogeotiff/rio-tiler/issues/105

When using one resampling methods for reading the data, it is (IMO) logical to use the same resampling method for mask reading (e.g via dataset_mask())

👇 This is an example when using different resampling methods for read and dataset_mask

Problem (bug) is when you use other resampling than nearest (e.g bilinear) the resulting mask from dataset_mask will contains values between 0-255.

Expected behavior and actual behavior.

Mask should always be binary 0 or 255.

Steps to reproduce the problem.

import rasterio
from rasterio.enums import Resampling

# Create Dataset
imgsize = 6000
arr = np.random.rand(3, imgsize, imgsize) * 100
arr = arr.astype(np.uint8) + 1
mask = np.full((imgsize, imgsize), 255).astype(np.uint8)
arr[:, 0: 1000, 0: 1000] = 0
mask[0: 1000, 0: 1000] = 0        

config = dict(GDAL_TIFF_INTERNAL_MASK=True)
with rasterio.Env(**config):
    params = {
        "driver": "GTiff",
        "blockxsize": 512,
        "blockysize": 512,
        "dtype": np.uint8,
        "count": 3,
        "height": imgsize,
        "width": imgsize,
    }
    with rasterio.open("uint8_mask.tif", "w", **params) as dst:
        dst.write(arr)
        dst.write_mask(mask)

with rasterio.open("uint8_mask.tif") as src_dst:
    tile = src_dst.read(window=((0,2000),(0,2000)), out_shape=(3, 256, 256), resampling=Resampling["bilinear"])
    mask = src_dst.dataset_mask(window=((0,2000),(0,2000)), out_shape=(256, 256), resampling=Resampling["bilinear"])
    mask_from_data = np.all(tile != 0, axis=0).astype(np.uint8) * 255

    assert mask == mask_from_data

Findings

In the code ☝️, dataset_mask is calling (confirmed with logging) https://github.com/mapbox/rasterio/blob/b9f34ee559039239c7c0c97bd911b466701a39cd/rasterio/_io.pyx#L560-L567, and when looking at the code there is a special case that was added (for boundless/window read) https://github.com/mapbox/rasterio/blob/b9f34ee559039239c7c0c97bd911b466701a39cd/rasterio/_io.pyx#L612-L623 to overcome the same behaviour.

Operating system

Mac OS

Rasterio version and provenance

Github master (1.0.24)

Notebook: https://gist.github.com/vincentsarago/551af3049f1a9708e020f23aa85c4803

[sgillies]

@vincentsarago isn't it obvious that resampling a mask will corrupt it? I'm not sure what you are suggesting we do. Should the resampling options be removed from dataset_mask? But then how would we get the resampled mask if we really, really needed it?

[vincentsarago]

Thanks @sgillies,

isn't it obvious that resampling a mask will corrupt it?

I'm not sure, but IMO a mask should always be 0 or 255 not something in between but to be honest I'm not 100%.

What I know is that I can't use bilinear resampling (in dataset_mask) for some data because it creates invalid mask, but then using nearest resampling doesn't create a valid mask either because it misses some part.

Especially when mask is created from nodata value, I guess you should be able to use the same resampling method on both read and dataset_mask.

what I'm proposing is that we apply the same logic from L622 where we cast the array to be either 0 or 255.

[sgillies]

@vincentsarago can you let me know what you think of this b50b1cc ?

[vincentsarago]

looks good to me @sgillies

[geowurster]

@sgillies @vincentsarago This is incorrect, mostly due to non-premultiplied alpha bands and partly due to GDAL weirdness that may or may not have been addressed in more recent versions.

I'm not sure, but IMO a mask should always be 0 or 255 not something in between but to be honest I'm not 100%.

---

I have provided some examples below but do not currently have an up-to-date Rasterio or OSGeo GDAL environment. The Rasterio code is using a wheel for:

$ rio --version && rio --gdal-version
1.1.2
2.4.3

and the OSGeo code is using:

$ gdalinfo --version
GDAL 2.4.1, released 2019/03/15

---

GDAL supports non-premultiplied alpha bands, or values with partial transparency between 0 and 255 independent of the underlying RGB values. For example the code below creates an RGBA TIFF (only green channel is populated) with an alpha band ranging from completely masked to unmasked. The top row is completely masked and the bottom is only partially masked. QGIS properly interprets this band and renders the image with a gradient. This behavior is also mentioned in GDAL RFC $#15, which introduced bit masks. The GeoTIFF driver's ALPHA creation option includes a few different options but I do not fully understand how all impact what data is produced by GDALGetMaskBand().

import numpy as np
import rasterio as rio

height = 255
width = 10

alpha = np.empty((height, width), dtype=np.uint8)

for idx in range(height):
    alpha[idx] = idx


data = np.array([
    # Red
    np.zeros_like(alpha),
    # Green
    np.full_like(alpha, 255),
    # Blue
    np.zeros_like(alpha),
    alpha
])

assert data.dtype == np.uint8


count, height, width = data.shape
profile = {
    'height': height,
    'width': width,
    'count': count,
    'driver': 'GTiff',
    'dtype': rio.uint8,
    'PHOTOMETRIC': 'RGB',
    'ALPHA': 'YES'
}

with rio.open('gradient.tif', 'w', **profile) as dst:
    dst.write(data)

The same behavior is present when stretching the gradient across a uint16 mask band. I had to adjust the stretch slightly given that 255 is a really dim value in uint16.

_Note that the above image has an incorrect filename. It is non-premultiplied, not premultiplied. Blame a bad Wikipedia page..._

---

Additionally, GDAL supports 2 band rasters where band 1 is data and band 2 has its color interpretation set to alpha and thus behaves (somewhat inconsistently) as a PER_DATASET mask. The specific values provided by the mask band are different depending on if the raster is UInt8 or UInt16. The behavior also may differ based on driver but I am most familiar with this happening with the GeoTIFF driver. Arguably an internal mask or a sidecar mask file would be more appropriate in this case, but internal masks are not supported by all image processing packages, and an alpha band still encapsulates the mask in a single file unlike a mask sidecar, so it has some valid uses.

from osgeo import gdal
import numpy as np


height = 10
width = 10

b1 = np.full((height, width), 34, dtype=np.uint16)
alpha = np.zeros_like(b1)
alpha[:5, :5] = 300
print("Alpha min/max:", alpha.min(), alpha.max())


ds = gdal.GetDriverByName('GTiff').Create(
    'uint16.tif', width, height, 2, gdal.GDT_UInt16)
ds.GetRasterBand(1).WriteArray(b1)
alpha_band = ds.GetRasterBand(2)
alpha_band.WriteArray(alpha)
alpha_band.SetColorInterpretation(gdal.GCI_AlphaBand)
band = None
ds = None


ds = gdal.Open('uint16.tif')
b1 = ds.GetRasterBand(1)
print("Band 1:", b1.ReadAsArray())
b2 = ds.GetRasterBand(2)
print("Band 2:", b2.ReadAsArray())
print("Band 1 mask:", b1.GetMaskBand().ReadAsArray())
print("Band 2 mask:", b2.GetMaskBand().ReadAsArray())

Code above produces:

Alpha min/max: 0 300
Band 1: [[34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]
 [34 34 34 34 34 34 34 34 34 34]]
Band 2: [[300 300 300 300 300   0   0   0   0   0]
 [300 300 300 300 300   0   0   0   0   0]
 [300 300 300 300 300   0   0   0   0   0]
 [300 300 300 300 300   0   0   0   0   0]
 [300 300 300 300 300   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]]
Band 1 mask: [[1 1 1 1 1 0 0 0 0 0]
 [1 1 1 1 1 0 0 0 0 0]
 [1 1 1 1 1 0 0 0 0 0]
 [1 1 1 1 1 0 0 0 0 0]
 [1 1 1 1 1 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]]
Band 2 mask: [[255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]
 [255 255 255 255 255 255 255 255 255 255]]

So what does that mean? It is important to remember that GDAL's mask bands are a GDAL construct and abstraction that sometimes sit on top of real raster bands containing real data that can be read outside of GDAL's masking infrastructure. GDALGetMaskBand() has all of the GDAL business logic that knows how to interpret nodata values, alpha bands, and other values (somewhat inconsistently).

GDAL wrote real data to both bands and then set the color interpretation in a manner that makes it later interpret one as a PER_DATASET mask. The first band is entirely populated with the value 34 and the second is mostly 0's aside from the upper left quadrant set to 300. Reading the first and second bands show that GDAL has written the data as instructed, however calling GDALGetBandMask() on the first band shows that GDAL is interpreting 0's as masked and any value >= 1 as unmasked. Oddly getting a mask band for the alpha band produces a totally transparent band even though that band is what supplies the PER_DATSET band, so one would expect it to match band 1's mask, given that band 2 still contains real values regardless of its color interpretation. The behavior may be different when using ds.GetMaskBand().WriteArray() but I did not pursue that avenue. I think this is due GDAL understanding that the raster has a PER_DATASET mask, which is intended to be read by calling GDALGetRasterBand() on the first band, which is something I _think_ I read somewhere but now I cannot find it in the GDAL docs. Maybe that came from trial by fire? Anyway, the GDALGetMaskBand() docs outline how GDAL derives a mask band based on the dataset's mask flags and includes this important (but maybe outdated?) note:

If there is no nodata value, but the dataset has an alpha band that seems to apply to this band (specific rules yet to be determined) and that is of type GDT_Byte then that alpha band will be returned, and the flags GMF_PER_DATASET and GMF_ALPHA will be returned in the flags.

To further complicate things changing the raster's datatype to Byte produces a mask band containing the original values (0's and 200's) insead of 0's and 1's. Presumably this band is still respecting the ALPHA creation option, which is non-premultiplied by default in modern GDAL versions:

Starting with GDAL 1.10, YES is an alias for NON-PREMULTIPLIED alpha, and the other values can be used.

but is probably also somehow affected by GDAL's general assumptions that 3 band Byte images have an RGB color interpretation and 4 band are RGBA. I did not check to see how QGIS interprets this image.

Another complicating factor is that it is possible to set multiple bands with alpha color interpretation, but $ gdalinfo does not report any mask flags for a dataset with this profile and I did not pursue the issue further. A dataset with multiple alpha bands is almost certainly a user error and an edge case way out in the weeds, but if Rasterio adopts smarter mask handling (https://github.com/mapbox/rasterio/issues/1456) then it is a case worth checking for.

GDAL definitely has its oddities here but I don't think Rasterio should necessarily hide them. If Rasterio wants to adopt better mask handling that seems fine, but I do think it is important to provide direct access to whatever GDALGetMaskBand() provides to stay consistent with GDAL – in addition to some docs warning about the pitfalls of either approach.

---

@vincentsarago As far as this goes:

What I know is that I can't use bilinear resampling (in dataset_mask) for some data because it creates invalid mask, but then using nearest resampling doesn't create a valid mask either because it misses some part.

Masks are discrete data, so interpolating with any algorithm other than nearest neighbor runs the risk of synthesizing new values. Interpolating with nearest neighbor does work but seems to mask off a few additional rows, which is just how masks work. It would be reasonable to use bilinear resampling on an alpha band containing transparencies but the edges still wouldn't be perfect. I'm not exactly familiar with how $ gdal_translate handles resampling in a case like this, but it is more aware of a dataset's mask in the context of all bands and may have additional logic to maintain consistence around edge pixels. In this case Rasterio's dataset_mask() mostly just calls GDALGetMaskBand(), at which point GDAL says "Hey here is this band that just happens to contain values that you should interpret as a mask. Hope you do the right thing with it." In this case it just assumes that you have decided that the resampling is appropriate.

---

@sgillies I don't think removing resampling from dataset_mask() is helpful, and IMO it is the user's responsibility to understand when resampling is appropriate and which algorithm to use, but adding a reminder in the docs that resampling with anything other than nearest neighbor may introduce new values can't hurt.

Should the resampling options be removed from dataset_mask? But then how would we get the resampled mask if we really, really needed it?

---

Extra pro tip for future work like https://github.com/mapbox/rasterio/commit/b50b1cca170256f5bcad9e7edc2999604966053a:

Avoid allocating int64 arrays in numpy.where() by providing Numpy scalars instead of Python integers:

>>> import numpy as np
>>> a = np.array([0, 1], dtype=np.uint8)
>>> print(a.dtype)
uint8
>>> a2 = np.where(a == 0, 2, 4)
>>> print(a2.dtype)
int64
>>> a3 = np.where(a == 0, np.uint8(2), np.uint8(4))
>>> print(a3.dtype)
uint8
>>> a4 = np.where(a == 0, a.dtype.type(2), a.dtype.type(4))
>>> print(a4.dtype)
uint8

I am less certain of the correctness of a.dtype.type() in the a4 example. Also if the 2nd and 3rd arguments are Numpy objects (array or scalar) then numpy.where() chooses a more sane default based on the datatypes and min/max values. Not sure why Numpy does not do the same for Python ints, but I'm guessing it is a historical consistency thing at this point.

[vincentsarago]

great write up @geowurster (really appreciated).
Just to make sure, you are saying we shouldn't change what was in place before?

if we keep what was in place I'll just add mask = np.where(mask != 0, 255, 0).astype("uint8") on my application.

[geowurster]

@vincentsarago Depends on the goal. Without actually testing this commit, I am pretty sure that prior to https://github.com/mapbox/rasterio/commit/6f91107f278a4c3e6d67b99c0d56c5f5f2ef56e0 (part of PR https://github.com/mapbox/rasterio/pull/1378) read_masks() (without a boundless read per https://github.com/mapbox/rasterio/issues/1266) produced whatever GDALGetMaskBand() provided. In this state Rasterio had support for non-premultiplied alpha bands, but also exposed users to the inconsistencies of GDAL's masking infrastructure. https://github.com/mapbox/rasterio/commit/6f91107f278a4c3e6d67b99c0d56c5f5f2ef56e0 obfuscated some of GDAL's inconsistencies at the expense of effectively adopting a new data model that assumes non-zero values are totally transparent and completely abandons support for non-premultiplied alpha bands. A non-premultiplied alpha band is useful for all of the reasons a transparency is useful when dealing with visual images, but can also be a measure of pixel quality. A real value can be present in each band, but the alpha band can be a representation of confidence in that value. 255 would be a very high confidence, 0 no confidences, 200 a pretty good confidence, etc.

https://github.com/mapbox/rasterio/commit/b50b1cca170256f5bcad9e7edc2999604966053a certainly matches the behavior introduced by https://github.com/mapbox/rasterio/commit/6f91107f278a4c3e6d67b99c0d56c5f5f2ef56e0, but IMO Rasterio is adopting its own data model that diverges from GDAL's just because GDAL's behavior is poorly defined, or poorly documented, or buggy, or inconsistent between versions. This is likely a problem that should be solved in GDAL and handled by Rasterio via some version checks that in turn provide consistency to users.

In general the consistency that Rasterio has provided is good and an enormous part of its success, but users porting code from the Python SWIG bindings, C/C++ API, older versions of Rasterio, or code from some other language's bindings are probably dealing with whatever GDALGetMaskBand() provides, so IMO it is important that Rasterio make the underlying behavior available, particularly for masks because they're just difficult to get right.

[vincentsarago]

@sgillies from my understanding of @geowurster comments, it's better not to change the actual behavior (meaning keep non binary mask values.

I've patched this in rio-tiler directly https://github.com/cogeotiff/rio-tiler/blob/e798adc1d521dc4b6fb5ae42432f553652287217/rio_tiler/utils.py#L477-L479

[sgillies]

@geowurster thank you, this is the leading candidate for Rasterio Issue Comment of the Year 2020 :smile:

[sgillies]

When I wrote the comments in 6f91107 I had seen some baffling masks in those VRTs. What I need to do is try to reproduce the problem and see if GDAL has or hasn't changed its behavior. I'm a big :+1: on walking back the mask model we've introduced almost by accident to give rasterio users more certainty.

[geowurster]

thank you, this is the leading candidate for Rasterio Issue Comment of the Year 2020

An award I would happily accept.

---

Some other thoughts on on Rasterio handles masks below.

dataset_mask()

Needs some mention of how any value > 0 is converted to 255 on read. I recognize this would be harder to change, but if this method continues to produce a boolean mask then representing it as a Numpy-style mask may provide some clarity for some users. It looks like write_mask() already supports boolean masks, however they are as the inverse of a Numpy mask and the behavior is undocumented.

Based on the method's name I would expect a dataset with GMF_PER_DATASET to produce something identical to read_masks(1). I would also expect a dataset with GMF_PER_DATASET and GMF_ALPHA with a non-premultiplied alpha band containing values other than 0 and 255 to produce those values unaltered.

Writing Per-Band Masks

write_mask() seems to only support writing GMF_PER_DATASET masks, or for configurations that support masks for each band, only the first band.

Baffling VRT Behavior

I think in some cases GMF_PER_DATASET mask bands for UInt16 datasets need to be written as the UInt16 max value (65535), but are then converted to 255 when read from GDALGetMaskBand()? I still haven't managed fully understand the GDAL mask world.

[geowurster]

For the sake of the record, this ticket is currently marked as closed due to https://github.com/mapbox/rasterio/commit/b50b1cca170256f5bcad9e7edc2999604966053a, however that commit contains changes that were reverted in https://github.com/mapbox/rasterio/pull/1871 along with some alterations in behavior.

@sgillies looks like this line in CHANGES.txt needs to be removed/updated:

https://github.com/mapbox/rasterio/blob/65757e96f9914157cfce3deadcf7faa73b85b352/CHANGES.txt#L10