Cell rasterized despite not being covered by polygon

[chau-intl]

What is the bug?

I have a simple polygon I am trying to rasterize and an extra cell is touched even though the polygon doesn't extend into it:

Both the lower vertices are exactly on the cell border with northing = 6114970.

Steps to reproduce the issue

Simple script to reproduce the issue:

from osgeo import gdal, ogr, osr

polygon = ogr.CreateGeometryFromWkt('POLYGON ((872512 6114970, 872512 6114970.4, 872513.3 6114970.4, 872513.1 6114970, 872512 6114970))')

srs = osr.SpatialReference()
srs.ImportFromEPSG(25832)

raster_ds = gdal.GetDriverByName('MEM').Create('', 8, 5, 1, gdal.GDT_Byte)
raster_ds.SetGeoTransform((872511.2, 0.4, 0, 6114971.2, 0, -0.4))
raster_ds.SetProjection(srs.ExportToWkt())

vector_ds = ogr.GetDriverByName('MEM').CreateDataSource('')
layer = vector_ds.CreateLayer('polygon', srs, ogr.wkbPolygon)

feature = ogr.Feature(layer.GetLayerDefn())
feature.SetGeometry(polygon)
layer.CreateFeature(feature)

err = gdal.RasterizeLayer(raster_ds, [1], layer, burn_values=[1], options=['ALL_TOUCHED'])

print(raster_ds.ReadAsArray())

Versions and provenance

OS: Windows 11
Python: 3.12.13 (from the OSGeo4W installer)
GDAL: 3.13.0 (from the OSGeo4W installer)

Additional context

No response

[jratike80]

Off-topic, but considering the name of the option. Shouldn't all_touched by the DE-9IM rules select all these

but the meaning is probably to compare if the feature and the pixel have intersecting interiors, in DE-9IM "Relate T********".

Perhaps there is too little tolerance when computing the relation of the vertex or segment of the polygon and the interpolated point on the edge of the raster cell. If it helps, those geometries are these, in order geometry of the pixel, vertex of the polygon to rasterize, segment of the polygon to rasterize

POLYGON (( 872512.8 6114970, 872512.8 6114969.6, 872513.2 6114969.6, 872513.2 6114970, 872512.8 6114970 ))
POINT (872513.1 6114970)
LINESTRING ( 872513.1 6114970, 872512 6114970 )

[chau-intl]

I think its great that you can polygonize from raster, then rasterize the polygons back using ALL_TOUCHED and still only get the same cells as in the original data set. In this report I'm wondering why the lower right vertex affects the marked cell but the lower left vertex doesn't?

If I understand #10606 correctly then even a microscopic expansion by a vertex into another cell would not affect the rasterization if the distance was below 1e-4? In my case the vertex stays on the border and should thus be covered by the threshold.

Does it have something to do with vertex order, since the first (and last) vertex doesn't expand the rasterization?

[chau-intl]

Does it have something to do with vertex order, since the first (and last) vertex doesn't expand the rasterization?

Apparently not since moving the vertices one place to the right like in the following WKT doesn't seem to affect the result:

POLYGON ((872513.1 6114970, 872512 6114970, 872512 6114970.4, 872513.3 6114970.4, 872513.1 6114970))

[jratike80]

I tried if a smaller pixel size could have an effect. Still POINT (872513.1 6114970) is interpreted to touch (what ever it means) the white pixel in down and right and paint it white. But now the top right vertex is also painting a pixel to the right from it white even there is no "true" intersection as Relate T********.

raster_ds = gdal.GetDriverByName('GTiff').Create('foo.tif', 64, 40, 1, gdal.GDT_Byte)

[jratike80]

Is this then actually a duplicate of #10606?

[rouault]

@dbaston I believe a reliable fix for that will involve reworking substantially the rasterization process and probably require a pre re-winding of polygons in one of CW or CCW order to be able to distinguish which side of a segment is the inner or outer.

[jratike80]

to be able to distinguish which side of a segment is the inner or outer

Does it matter so much when the case is "all-touched"?

[jratike80]

to be able to distinguish which side of a segment is the inner or outer

Does it matter so much when the case is "all-touched"?

Correction: Yes I think it does matter in the corner cases ;) when the geometry is only touching a raster cell. Then there is a difference if the interior of the cell is totally inside or totally outside.

[mdsumner]

just FYI, I've got a package that adapts exactextract to an always-sparse approach (originally adapted from fasterize approx method, and finished with Claude).

hypertidy/controlledburn#10

The approach in burn_scanline() also avoids the epsilon values that ee gives.

There's a long way to go, it does enough of what I want for R - but the core could be stripped out as a C++ lib to share - that's my long term plan.