Improve AggregateRaster for small source geometries

[byrman]

This pull request solves https://github.com/nens/lizard-nxt/issues/5814 by using a different strategy in AggregateRaster.

An alternative solution is to enable the ALL_TOUCHED rasterization option in GDAL (affecting other calculations too): https://gdal.org/programs/gdal_rasterize.html#cmdoption-gdal_rasterize-at.

[caspervdw]

I recall making the # snap the extent to (0, 0) to prevent subpixel shift part to make label computations deterministic.

There was problem that the outcome of an aggregation value depended on other geometries being present. This is because the extent is computed based on all geometries in the request.

I think that effect is back now as x2, y1 are not snapped anymore to (0,0).

[byrman]

Have you seen these results: https://github.com/nens/lizard-nxt/issues/5814#issuecomment-1479321495?

[byrman]

I think that effect is back now as x2, y1 are not snapped anymore to (0,0)

"think"?

[caspervdw]

Have you seen these results: https://github.com/nens/lizard-nxt/issues/5814#issuecomment-1479321495?

I think that effect is back now as x2, y1 are not snapped anymore to (0,0)

"think"?

Yes I read the issue. I am afraid I don't understand what is happening.
Also I can't oversee the effects of this change completely.

Does this test fail now?

dask-geomodeling/dask_geomodeling/tests/test_geometry.py

Line 882 in 111275c

def test_snap_bbox(self):

I just want to be sure that this solution doesn't break other use cases.

[arjanverkerk]

I recall making the # snap the extent to (0, 0) to prevent subpixel shift part to make label computations deterministic.

There was problem that the outcome of an aggregation value depended on other geometries being present. This is because the extent is computed based on all geometries in the request.

I think that effect is back now as x2, y1 are not snapped anymore to (0,0).

Ah, the clash of use cases... What is the subpixel shift we try to avoid? I'd guess that what is meant is that the rasterization of a particular feature does not depend on the envelope of the group of features that are rasterized together. But by anchoring anchoring to (0, 0), the tiniest change in pixel_size may result in a big shift of the bounding box requested to the raster source in the current situation which may have big effects on the rasterization.

So should we add an option here to choose between snapping methods? We also may need an option to enable the ALL_TOUCHED gdal option.

[caspervdw]

I'd guess that what is meant is that the rasterization of a particular feature does not depend on the envelope of the group of features that are rasterized together.

Exactly. Well put!

the tiniest change in pixel_size may result in a big shift of the bounding box requested to the raster source in the current situation which may have big effects on the rasterization.

I don't get that part. The "big shift" is at maximum 1 cell, right? From what perspective is that big? And in any case, the geometry will be enclosed by the raster.

What happens in rasterization to a geometry much smaller than 1 cell?

[byrman]

Does this test fail now?

All pull requests fail, not just mine. I don't know what's wrong with the test setup. I would be happy the run the tests locally in a container (but no Docker files are present).

[caspervdw]

Does this test fail now?

All pull requests fail, not just mine. I don't know what's wrong with the test setup. I would be happy the run the tests locally in a container (but no Docker files are present).

Just make a virtualenv, any python version will do.

[byrman]

Just make a virtualenv, any python version will do.

Doesn't the project require more than that (e.g. gdal libraries)? I prefer to keep my (macOS) host clean, because I work on many different projects.

[caspervdw]

Just make a virtualenv, any python version will do.

Doesn't the project require more than that (e.g. gdal libraries)? I prefer to keep my (macOS) host clean, because I work on many different projects.

Yes you do need to have the GDAL binaries, on Ubuntu that would be apt install libgdal-dev.

[arjanverkerk]

I don't get that part. The "big shift" is at maximum 1 cell, right? From what perspective is that big? And in any case, the geometry will be enclosed by the raster.

From the perspective of the cells near the edge of a feature it means everything, it means to be included or not. We're talking about the rain raster here, which has 1000 m cells. Of course, one should pick a much smaller pixel_size (supersampling) when the geometry's area is in the same order as the pixel area, and that is what I suggested and @byrman tried. But then he noticed how small pixel_size changes affected the result very much.

What happens in rasterization to a geometry much smaller than 1 cell?

Good one. In this use case it is not clearly defined. We have a small diamond shaped geometry. When near or at the center of a cell, it gets rasterized. But when it is over the intersection between 4 cells and away from centers, it is not rasterized (Bresenhams line algorithm, not exactly meant for this use case I guess) unless 'ALL_TOUCHED=TRUE'.