Overlay_polygon

Geometric functions overlay_polygon

Since GeoDMS 18.0.0 the unprefixed name overlay_polygon no longer exists. It was the boost polygon implementation, which lives on as bp_overlay_polygon with the same arguments; a configuration that still uses the old name fails with an unknown-operator error and is repaired by adding the bp_ prefix. The other implementations are bg_overlay_polygon, cgal_overlay_polygon and geos_overlay_polygon; see polygon operators for the choice between them, and for the polygon_connectivity family that finds the same pairs without the intersection geometry.

syntax

• overlay_polygon (first_polygon_data_item, second_polygon_data_item), until GeoDMS 18.0.0; now bp_overlay_polygon for integer coordinates, using the boost polygon library
• bg_overlay_polygon (first_polygon_data_item, second_polygon_data_item) works also on floating point coordinates, using the boost geometry library.

description

overlay_polygon(first_polygon_data_item, second_polygon_data_item) results in a new uint32 domain unit with an entry for each intersecting part of the first_polygon_data_item and the second_polygon_data_item.

The function generates three subitems for the new domain unit:

• geometry: the geometry of the resulting polygons (in the figure the yellow polygons). This attribute has the same values unit as the first_polygon_data_item and second_polygon_data_item attributes.
• first_rel: a relation for the new domain unit towards the domain of the first_polygon_data_item.
• second_rel: a relation for the new domain unit towards the domain of the second_polygon_data_item.

In QGIS, this operation is known as intersect.

It is similar to:

unit<uint32> CartesianProduct := combine(DomainUnit(first_polygon_data_item), DomainUnit(second_polygon_data_item))
{
    attribute<ValuesUnit(first_polygon_data_item)> Geometry(poly) := first_polygon_data_item[first_rel] * second_polygon_data_item[second_rel];
    attribute<bool>   IsIntersecting := area(Geometry, float64) > 0.0;
}
unit<uint32> result := select(CartesianProduct/IsIntersecting)
{
    attribute<DomainUnit(first_polygon_data_item)> first_rel := collect_by_cond(., CartesianProduct/first_rel);
    attribute<DomainUnit(second_polygon_data_item)> second_rel := collect_by_cond(., CartesianProduct/second_rel);
    attribute<ValuesUnit(first_polygon_data_item)> Geometry(poly) := collect_by_cond(., CartesianProduct/Geometry);
}

applies to

attributes first_polygon_data_item and second_polygon_data_item with an ipoint or spoint value type.

conditions

1. The composition type of the first_polygon_data_item and second_polygon_data_item arguments needs to be polygon.
2. The values unit of the first_polygon_data_item and second_polygon_data_item arguments must match.
3. The order of the points in the first_polygon_data_item and second_polygon_data_item needs to be clockwise for exterior bounds and counter clockwise for holes in polygons (right-hand-rule).

This function results in problems for (integer) coordinates larger than 2^25 (after translation where the first point is moved to (0, 0)). If your integer coordinates for instance represent mm, 2^25[mm] = about 33 [km]. The reason is that for calculating intersections, products of coordinates are calculated and casted to float64 with a 53 bits mantissa (in the development/test environment of boost::polygon these were float80 values with a 64 bits mantissa). We advise to keep the size of your integer coordinates for polygons limited and for instance do not use a mm precision for country borders (meter or kilometer might be sufficient).

since version

7.042

performance

O(n * m * log(n + m)) where n = number of first polygons and m = number of second polygons. Uses spatial indexing to prune non-overlapping pairs. Complexity per intersecting pair depends on vertex count. Result size can be large for highly overlapping datasets.

The work is spread over the threads by the polygons of the first argument, so put the argument with the most, and the smallest, polygons first; see Which argument first under polygon operators.

example

unit<uint32> intersect := bp_overlay_polygon(Building/geometry, District/geometry);

see also

• bp_overlay_polygon - the current name of this operator
• polygon_connectivity - which pairs of polygons within one set touch or overlap, as relations only: the question overlay_polygon answers for two sets, without the intersection geometry
• bp_polygon_connectivity, bg_polygon_connectivity, cgal_polygon_connectivity, geos_polygon_connectivity - the same pairs as an overlay finds, as relations only
• mul (polygon intersection) - the intersection of two polygons element by element, for pairs that are already known
• polygon operators - the comparison of the operators that find pairs of polygons
• Border polygons - configuration example