Rasterio: Unexpected behavior by rasterio.mask.mask with all_touched=True (too many pixels are cropped out)

[atorch]

Expected behavior and actual behavior.

Based on the documentation at https://rasterio.readthedocs.io/en/latest/api/rasterio.mask.html, I expect rasterio.mask.mask with all_touched=True to always include all pixels that touch any of the shapes that I pass in.

Steps to reproduce the problem.

import pprint
from tempfile import NamedTemporaryFile

import numpy as np
import rasterio

from rasterio import Affine
from rasterio.crs import CRS
from rasterio.mask import mask
from shapely.geometry import mapping, Polygon

ORIGIN_X = 408225.0
ORIGIN_Y = 1985985.0

CDL_RESOLUTION = 30.0

AFFINE = Affine(
    CDL_RESOLUTION,
    0.0,
    ORIGIN_X,
    0.0,
    -CDL_RESOLUTION,
    ORIGIN_Y,
)

CDL_CRS = CRS(
    {
        'proj': 'aea',
        'lat_1': 29.5,
        'lat_2': 45.5,
        'lat_0': 23,
        'lon_0': -96,
        'x_0': 0,
        'y_0': 0,
        'ellps':
        'GRS80',
        'towgs84': '0,0,0,0,0,0,0',
        'units': 'm',
        'no_defs': True,
    }
)


def save_toy_raster(origin_x, origin_y, outfile, height=8, width=10, dtype='uint8'):
    raster_meta = {
        'count': 1,
        'crs': CDL_CRS,
        'driver': 'GTiff',
        'dtype': dtype,
        'height': height,
        'width': width,
        'nodata': None,
    }

    raster_meta['transform'] = AFFINE

    raster_values = np.zeros((height, width), dtype=dtype)  # Y, x

    raster_values[0, 0] = 1
    raster_values[0, 1] = 2
    raster_values[0, 3] = 7
    raster_values[4, 0] = 9

    with rasterio.open(outfile, 'w', **raster_meta) as outfile:
        outfile.write(raster_values, indexes=1)

    return


def get_mask_and_affine(raster, shape):
    masked_array, affine = mask(
        raster,
        [mapping(shape)],
        filled=False,
        crop=True,
        all_touched=True,
    )

    return masked_array, affine


def get_example_geometry():
    top_right_coord = [408318.19810377056, 1985962.5740539758]

    shape = Polygon(
        [
            [408240.36536500533, 1985846.4332173148],
            [408234.6170334989, 1985958.2779089767],
            top_right_coord,
            [408323.94761217793, 1985850.7294228084],
            [408240.36536500533, 1985846.4332173148],
        ],
    )

    return shape, top_right_coord


def main():
    pp = pprint.PrettyPrinter(indent=4)

    raster_file = NamedTemporaryFile().name
    save_toy_raster(
        origin_x=ORIGIN_X,
        origin_y=ORIGIN_Y,
        outfile=raster_file,
    )

    raster = rasterio.open(raster_file)
    print('raster values:')
    pp.pprint(raster.read())

    shape, top_right_coord = get_example_geometry()

    masked_array, _affine = get_mask_and_affine(raster, shape)

    value_at_top_right = list(raster.sample([top_right_coord]))[0]

    print('masked_array:')
    pp.pprint(masked_array)

    # masked_array:
    # masked_array(
    # data=[[[1, 2, 0],
    #         [0, 0, 0],
    #         [0, 0, 0],
    #         [0, 0, 0]]],
    # mask=[[[False, False, False],
    #         [False, False, False],
    #         [False, False, False],
    #         [False, False, False]]],
    # fill_value=999999,
    # dtype=uint8)

    print('value_at_top_right:')
    pp.pprint(value_at_top_right)

    # value_at_top_right:
    # array([7], dtype=uint8)


if __name__ == '__main__':
    main()

In this example, we define a toy raster (whose projection and resolution mimic https://nassgeodata.gmu.edu/CropScape/, although that doesn't matter for this bug report). Notice that the raster's [0, 3] value is 7, and that our example geometry touches the [0, 3] pixel. That pixel is missing from masked_array, i.e. it appears to have been incorrectly cropped out.

Qgis screenshot: https://imgur.com/a/wM8sY8Y. The bright pixel at the top right intersects the polygon (in red), but that pixel appears to be cropped out by rasterio.mask.mask even with all_touched=True. This appears to contradict the docs, which say "all_touched (bool (opt)) – Include a pixel in the mask if it touches any of the shapes."

Operating system

Ubuntu 17.10

Rasterio version and provenance

I noticed this behavior on rasterio 1.0a12, and also saw it on rasterio 1.0.2.

[atorch]

Here is another piece of code that reproduces the issue I am describing:

import numpy as np
import rasterio
from tempfile import NamedTemporaryFile

from rasterio import Affine
from rasterio.crs import CRS
from rasterio.mask import mask
from shapely.geometry import mapping, Polygon

import numpy as np

AFFINE = Affine(
    1.0,
    0.0,
    -1.0,  # X origin
    0.0,
    -1.0,
    1.0,  # Y origin
)

CRS = CRS.from_epsg(4326)

def save_toy_raster(outfile, height=3, width=3, dtype='uint8'):
    raster_meta = {
        'count': 1,
        'crs': CRS,
        'driver': 'GTiff',
        'dtype': dtype,
        'height': height,
        'width': width,
        'nodata': None,
        'transform': AFFINE
    }

    raster_values = np.zeros(height * width, dtype=dtype)
    for idx, value in enumerate(raster_values):
        raster_values[idx] = idx
    raster_values = np.reshape(raster_values, (height, width))

    with rasterio.open(outfile, 'w', **raster_meta) as outfile:
        outfile.write(raster_values, indexes=1)

def get_mask_and_affine(raster, shape):
    masked_array, affine = mask(
        raster,
        [mapping(shape)],
        filled=False,
        crop=True,
        pad=True,
        all_touched=True
    )

    return masked_array, affine

def get_example_geometry(shift_x=0, shift_y=0):

    shape = Polygon([
        [0. + shift_x, 0. + shift_y],
        [1. + shift_x, 0. + shift_y],
        [1. + shift_x, -1. + shift_y],
        [0. + shift_x, -1. + shift_y],
        [0. + shift_x, 0. + shift_y],
    ])

    return shape

def test_rasterio_mask(raster_file):

    raster = rasterio.open(raster_file)
    print('Raw raster: ', raster.read(1))

    for dy in np.arange(-0.1, 0.1, 0.1):
        for dx in np.arange(-0.1, 0.1, 0.1):
            geometry = get_example_geometry(dx, dy)
            masked_array, affine = get_mask_and_affine(
                raster,
                geometry
            )
            print(f'(dx, dy): ({dx}, {dy})')
            coords_x, coords_y = geometry.exterior.coords.xy
            print('geometry (x, y): ', list(zip(coords_x.tolist(), coords_y.tolist())))
            print('pixels: ', masked_array)
            print('affine: ', affine.to_gdal())
            print()

def main():

    raster_file = NamedTemporaryFile().name
    save_toy_raster(outfile=raster_file)
    test_rasterio_mask(raster_file)

if __name__ == '__main__':
    main()

When (dx, dy): (-0.1, -0.1), the geometry should touch four pixels in the raster (since the geometry is 1-by-1 pixels and aligned/on a grid with (0, 0)). Instead, the results from rasterio.mask.mask say that the geometry only touches two pixels.

[sgillies]

Yes, looks like a bug to me. Thanks for the clear report and code @atorch !

[sgillies]

@atorch I think I found a bug in geometry_window() that can leave the extents one row or column short. Here are the results I am getting now using your code:

Raw raster:  [[0 1 2]
 [3 4 5]
 [6 7 8]]
(dx, dy): (-0.1, -0.1)
geometry (x, y):  [(-0.1, -0.1), (0.9, -0.1), (0.9, -1.1), (-0.1, -1.1), (-0.1, -0.1)]
pixels:  [[[-- -- --]
  [3 4 --]
  [6 7 --]]]
affine:  (-1.0, 1.0, 0.0, 1.0, 0.0, -1.0)

(dx, dy): (0.0, -0.1)
geometry (x, y):  [(0.0, -0.1), (1.0, -0.1), (1.0, -1.1), (0.0, -1.1), (0.0, -0.1)]
pixels:  [[[-- -- --]
  [-- 4 5]
  [-- 7 8]]]
affine:  (-1.0, 1.0, 0.0, 1.0, 0.0, -1.0)

(dx, dy): (-0.1, 0.0)
geometry (x, y):  [(-0.1, 0.0), (0.9, 0.0), (0.9, -1.0), (-0.1, -1.0), (-0.1, 0.0)]
pixels:  [[[-- -- --]
  [3 4 --]
  [6 7 --]]]
affine:  (-1.0, 1.0, 0.0, 1.0, 0.0, -1.0)

(dx, dy): (0.0, 0.0)
geometry (x, y):  [(0.0, 0.0), (1.0, 0.0), (1.0, -1.0), (0.0, -1.0), (0.0, 0.0)]
pixels:  [[[-- -- --]
  [-- 4 5]
  [-- 7 8]]]
affine:  (-1.0, 1.0, 0.0, 1.0, 0.0, -1.0)

4 pixels instead of 2 in all cases.