More edge behaviour issues with RasterizeLayer and ALL_TOUCHED

[mathause]

What is the bug?

Unfortunately there are some more "edge-cases" left after #8918 (thanks for fixing that so quickly!). Additional grid points are found as 'touched' if the coordinates of the other side of the polygon are not at the raster edges.

In the example below I think there should be only 64 cells.

Steps to reproduce the issue

from osgeo import gdal, ogr, osr

polygons = {
    "a": "POLYGON ((9 1.5, 9 8.5, 0.9999 8.5, 0.9999 1.5, 9 1.5))",
    "b": "POLYGON ((9.0001 1.5, 9.0001 8.5, 1 8.5, 1 1.5, 9.0001 1.5))",
    "c": "POLYGON ((8.5 0.9999, 8.5 9, 1.5 9, 1.5 0.9999, 8.5 0.9999))",
    "d": "POLYGON ((8.5 1, 8.5 9.0001, 1.5 9.0001, 1.5 1, 8.5 1))",
}

for i, (key, wkt_geom) in enumerate(polygons.items()):


    # Setup working spatial reference
    sr_wkt = 'LOCAL_CS["arbitrary"]'
    sr = osr.SpatialReference(sr_wkt)

    # Create a memory raster to rasterize into.

    target_ds = gdal.GetDriverByName("MEM").Create("", 10, 10, 3, gdal.GDT_Byte)
    target_ds.SetGeoTransform((0, 1, 0, 0, 0, 1))

    target_ds.SetProjection(sr_wkt)

    # Create a memory layer to rasterize from.

    rast_ogr_ds = ogr.GetDriverByName("Memory").CreateDataSource("wrk")
    rast_mem_lyr = rast_ogr_ds.CreateLayer("poly", srs=sr)


    feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())
    feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom))

    rast_mem_lyr.CreateFeature(feat)

    # Run the algorithm.

    err = gdal.RasterizeLayer(
        target_ds,
        [2, 1],
        rast_mem_lyr,
        burn_values=[80, -1],
        options=["ALL_TOUCHED"],
    )

    assert err == 0, "got non-zero result code from RasterizeLayer"

    gdal.GetDriverByName("GTiff").CreateCopy(f"rasterize_{i}.tif", target_ds)

    n_cells = target_ds.GetRasterBand(2).ReadAsArray().sum() / 80

    print(f"Polygon {key}: {n_cells}")

Versions and provenance

GDAL 3.9.1, released 2024/06/22

installed via conda-forge

Additional context

• for "a" miny must be between 1 + 1e-4 and 2 - 1e-4", otherwise it does not happen.
• "a" and "b" (and "c" and "d") don't exactly behave the same:
  • for "a" minx must between 1 - 1e-4 and 1 - eps
  • for "b" maxx must be exactly 9.+1.e-4 (+- np.spacing(9) or so)

I used rasterio to create the mask for the plot as I am more used to it's interface.

import affine
import matplotlib.pyplot as plt
import rasterio.features
import shapely
import shapely.plotting

transform = affine.Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)

polygons = {
    "a": shapely.geometry.box(1-1e-4, 1.5, 9., 8.5),
    "b": shapely.geometry.box(1, 1.5, 9.+1e-4, 8.5),
    "c": shapely.geometry.box(1.5, 1-1e-4, 8.5, 9),
    "d": shapely.geometry.box(1.5, 1, 8.5, 9+1e-4),
}

f, axs = plt.subplots(2, 2, sharex=True, sharey=True)
f.set_size_inches(6, 6)
axs = axs.flatten()
# =============================

for i, (key, poly) in enumerate(polygons.items()):

    print(f'"{key}": "{shapely.to_wkt(poly)}",')

    raster = rasterio.features.rasterize(
        ((poly, 1),),
        out_shape=(10, 10),
        fill=0,
        all_touched=True,
        transform=transform,
    )

    ax = axs[i]
    shapely.plotting.plot_polygon(poly, ax, facecolor="none", edgecolor="r", zorder=5, add_points=False)
    ax.pcolormesh(raster, ec="0.1", lw=0.125)
    ax.set_aspect("equal")
    ax.set_title(key)

[jratike80]

Isn't it rather so that too few cells get selected? For example I would expect that the polygon "a" would select the whole leftmost column because those cells have a non-zero intersection area. This is based on an assumption that "ALL_TOUCHED" really means "ALL_WITH_AREAL_INTERSECTION" that may be wrong.

[rouault]

because those cells have a non-zero intersection area.

There is a threshold of 1e-4 (applied in parts of the code paths) to take into account for numerical imprecision so that very tiny intersections are considered to be non significant. So in theory I'd agree with the diagnostic of the OP. Fixing the code to achieve that in all situations is another thing.