Skip to content

Commit

Permalink
gdal_rasterize: fix inverse rasterization of polygon nested inside an…
Browse files Browse the repository at this point in the history 
…other one. Requires GEOS enabled build (fixes OSGeo#8689)
  • Loading branch information
@rouault
rouault committed Nov 10, 2023
1 parent 170bc25 commit 2e3130c
Show file tree
Hide file tree
Showing 3 changed files with 209 additions and 40 deletions.
113 changes: 73 additions & 40 deletions apps/gdal_rasterize_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ static void InvertGeometries(GDALDatasetH hDstDS,
std::vector<OGRGeometryH> &ahGeometries)

{
OGRGeometryH hInvertMultiPolygon = OGR_G_CreateGeometry(wkbMultiPolygon);
OGRMultiPolygon *poInvertMP = new OGRMultiPolygon();

/* -------------------------------------------------------------------- */
/* Create a ring that is a bit outside the raster dataset. */
Expand All @@ -82,46 +82,83 @@ static void InvertGeometries(GDALDatasetH hDstDS,
double adfGeoTransform[6] = {};
GDALGetGeoTransform(hDstDS, adfGeoTransform);

OGRGeometryH hUniverseRing = OGR_G_CreateGeometry(wkbLinearRing);
OGRLinearRing *poUniverseRing = new OGRLinearRing();

OGR_G_AddPoint_2D(
hUniverseRing,
poUniverseRing->addPoint(
adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);

OGR_G_AddPoint_2D(hUniverseRing,
adfGeoTransform[0] + brx * adfGeoTransform[1] +
-2 * adfGeoTransform[2],
adfGeoTransform[3] + brx * adfGeoTransform[4] +
-2 * adfGeoTransform[5]);

OGR_G_AddPoint_2D(hUniverseRing,
adfGeoTransform[0] + brx * adfGeoTransform[1] +
bry * adfGeoTransform[2],
adfGeoTransform[3] + brx * adfGeoTransform[4] +
bry * adfGeoTransform[5]);

OGR_G_AddPoint_2D(hUniverseRing,
adfGeoTransform[0] + -2 * adfGeoTransform[1] +
bry * adfGeoTransform[2],
adfGeoTransform[3] + -2 * adfGeoTransform[4] +
bry * adfGeoTransform[5]);

OGR_G_AddPoint_2D(
hUniverseRing,
poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
-2 * adfGeoTransform[2],
adfGeoTransform[3] + brx * adfGeoTransform[4] +
-2 * adfGeoTransform[5]);

poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
bry * adfGeoTransform[2],
adfGeoTransform[3] + brx * adfGeoTransform[4] +
bry * adfGeoTransform[5]);

poUniverseRing->addPoint(adfGeoTransform[0] + -2 * adfGeoTransform[1] +
bry * adfGeoTransform[2],
adfGeoTransform[3] + -2 * adfGeoTransform[4] +
bry * adfGeoTransform[5]);

poUniverseRing->addPoint(
adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);

OGRGeometryH hUniversePoly = OGR_G_CreateGeometry(wkbPolygon);
OGR_G_AddGeometryDirectly(hUniversePoly, hUniverseRing);
OGRPolygon *poUniversePoly = new OGRPolygon();
poUniversePoly->addRingDirectly(poUniverseRing);
poInvertMP->addGeometryDirectly(poUniversePoly);

OGR_G_AddGeometryDirectly(hInvertMultiPolygon, hUniversePoly);
bool bFoundNonPoly = false;
// If we have GEOS, use it to "subtract" each polygon from the universe
// multipolygon
if (OGRGeometryFactory::haveGEOS())
{
OGRGeometry *poInvertMPAsGeom = poInvertMP;
poInvertMP = nullptr;
CPL_IGNORE_RET_VAL(poInvertMP);
for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
{
auto poGeom = OGRGeometry::FromHandle(ahGeometries[iGeom]);
const auto eGType = OGR_GT_Flatten(poGeom->getGeometryType());
if (eGType != wkbPolygon && eGType != wkbMultiPolygon)
{
if (!bFoundNonPoly)
{
bFoundNonPoly = true;
CPLError(CE_Warning, CPLE_AppDefined,
"Ignoring non-polygon geometries in -i mode");
}
}
else
{
auto poNewGeom = poInvertMPAsGeom->Difference(poGeom);
if (poNewGeom)
{
delete poInvertMPAsGeom;
poInvertMPAsGeom = poNewGeom;
}
}

delete poGeom;
}

ahGeometries.resize(1);
ahGeometries[0] = OGRGeometry::ToHandle(poInvertMPAsGeom);
return;
}

/* -------------------------------------------------------------------- */
/* Add outer rings of polygons as inner rings of hUniversePoly */
/* and inner rings as sub-polygons. */
/* If we don't have GEOS, add outer rings of polygons as inner */
/* rings of poUniversePoly and inner rings as sub-polygons. Note */
/* that this only works properly if the polygons are disjoint, in */
/* the sense that the outer ring of any polygon is not inside the */
/* outer ring of another one. So the scenario of */
/* https://github.com/OSGeo/gdal/issues/8689 with an "island" in */
/* the middle of a hole will not work properly. */
/* -------------------------------------------------------------------- */
bool bFoundNonPoly = false;
for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
{
const auto eGType =
Expand All @@ -139,19 +176,15 @@ static void InvertGeometries(GDALDatasetH hDstDS,
}

const auto ProcessPoly =
[hUniversePoly, hInvertMultiPolygon](OGRPolygon *poPoly)
[poUniversePoly, poInvertMP](OGRPolygon *poPoly)
{
for (int i = poPoly->getNumInteriorRings() - 1; i >= 0; --i)
{
auto poNewPoly = new OGRPolygon();
poNewPoly->addRingDirectly(poPoly->stealInteriorRing(i));
OGRGeometry::FromHandle(hInvertMultiPolygon)
->toMultiPolygon()
->addGeometryDirectly(poNewPoly);
poInvertMP->addGeometryDirectly(poNewPoly);
}
OGRGeometry::FromHandle(hUniversePoly)
->toPolygon()
->addRingDirectly(poPoly->stealExteriorRing());
poUniversePoly->addRingDirectly(poPoly->stealExteriorRing());
};

if (eGType == wkbPolygon)
Expand All @@ -165,16 +198,16 @@ static void InvertGeometries(GDALDatasetH hDstDS,
{
auto poMulti =
OGRGeometry::FromHandle(ahGeometries[iGeom])->toMultiPolygon();
for (int i = 0; i < poMulti->getNumGeometries(); i++)
for (auto *poPoly : *poMulti)
{
ProcessPoly(poMulti->getGeometryRef(i));
ProcessPoly(poPoly);
}
delete poMulti;
}
}

ahGeometries.resize(1);
ahGeometries[0] = hInvertMultiPolygon;
ahGeometries[0] = OGRGeometry::ToHandle(poInvertMP);
}

/************************************************************************/
Expand Down
130 changes: 130 additions & 0 deletions autotest/utilities/test_gdal_rasterize_lib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -563,6 +563,136 @@ def test_gdal_rasterize_lib_inverse():
assert struct.unpack("B" * 121, target_ds.ReadRaster()) == expected


###############################################################################


@pytest.mark.require_geos
def test_gdal_rasterize_lib_inverse_nested_polygons():

# Create a memory raster to rasterize into.
target_ds = gdal.GetDriverByName("MEM").Create("", 10, 10, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((0, 1, 0, 10, 0, -1))

# Create a memory layer to rasterize from.
vector_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
rast_mem_lyr = vector_ds.CreateLayer("poly")

# Add a multi polygon with nested polygons
feat = ogr.Feature(rast_mem_lyr.GetLayerDefn())
wkt = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0),(2 2,2 8,8 8,8 2,2 2)),((4 4,4 6,6 6,6 4,4 4)))"
feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt))

rast_mem_lyr.CreateFeature(feat)

gdal.Rasterize(target_ds, vector_ds, burnValues=[1], inverse=True)

expected = (
0,
0,
0,
0,
0,
0,
0,
0,
0,
0, ################################
0,
0,
0,
0,
0,
0,
0,
0,
0,
0, ################################
0,
0,
1,
1,
1,
1,
1,
1,
0,
0, ################################
0,
0,
1,
1,
1,
1,
1,
1,
0,
0, ################################
0,
0,
1,
1,
0,
0,
1,
1,
0,
0, ################################
0,
0,
1,
1,
0,
0,
1,
1,
0,
0, ################################
0,
0,
1,
1,
1,
1,
1,
1,
0,
0, ################################
0,
0,
1,
1,
1,
1,
1,
1,
0,
0, ################################
0,
0,
0,
0,
0,
0,
0,
0,
0,
0, ################################
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
)

got = struct.unpack("B" * (10 * 10), target_ds.ReadRaster())
assert got == expected


###############################################################################
# Test rasterization of a 64 bit integer attribute

Expand Down
6 changes: 6 additions & 0 deletions doc/source/programs/gdal_rasterize.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,12 @@ raster data is only supported since GDAL 2.1.0.
with the first feature into all parts of the image *not* inside the
provided polygon.

.. note::

When the vector features contain a polygon nested within another polygon
(like an island in a lake), GDAL must be built against GEOS to get
correct results.

.. option:: -at

Enables the ALL_TOUCHED rasterization option so that all pixels touched
Expand Down

0 comments on commit 2e3130c

Please sign in to comment.