More edge behaviour issues with RasterizeLayer and ALL_TOUCHED

mathause commented 4 weeks ago

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

rasterize_alltouched_edge

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. ```python 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 commented 4 weeks ago

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 commented 4 weeks ago

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.

OSGeo / gdal