Cleaning invalid polygons

Configuration examples Cleaning invalid polygons

Polygon data from other sources is often not valid: a ring crosses itself, a lake is wound the same way as the shell that contains it, two parcels that should share a boundary are a few millimetres apart, or a ring simply does not close. The polygon operators of the four geometry libraries need valid input and either refuse such geometry or repair it first, which changes it in ways the model cannot see. The dms polygon operators read it directly, under the even-odd rule, and take an optional grid size as the tolerance at which vertices that are too close to be meant apart become one.

This example shows four tasks on that basis. Every number below is checked in the regression configuration that ships as examples\testcases\oper_dms_overlay.dms.

the geometry used here 

container CleaningInvalidPolygons
{
	unit<uint32> Solo : nrofrows = 1;
	unit<uint32> RPt  : nrofrows = 5;
	attribute<Solo>   seq (RPt) : [ 0, 0, 0, 0, 0 ];
	attribute<uint32> ord (RPt) : [ 0, 1, 2, 3, 4 ];

	// a square [0,100] x [0,100], clockwise and closed: valid geometry
	attribute<int32>  sq_x (RPt) : [ 0,   0, 100, 100, 0 ];
	attribute<int32>  sq_y (RPt) : [ 0, 100, 100,   0, 0 ];
	attribute<ipoint> square (Solo, polygon) := points2polygon(point_xy(sq_x, sq_y, ipoint), seq, ord);

	// the same corners in an order that makes the ring cross itself at (50, 50): a bow tie
	unit<uint32> BowPt : nrofrows = 4;
	attribute<int32>  bow_x (BowPt) : [ 0, 100, 100,   0 ];
	attribute<int32>  bow_y (BowPt) : [ 0, 100,   0, 100 ];
	attribute<Solo>   bow_seq (BowPt) : [ 0, 0, 0, 0 ];
	attribute<uint32> bow_ord (BowPt) : [ 0, 1, 2, 3 ];
	attribute<ipoint> bowtie (Solo, polygon) := points2polygon(point_xy(bow_x, bow_y, ipoint), bow_seq, bow_ord);

	// a neighbour that misses the square's right edge by three units: [103,203] x [0,100]
	attribute<int32>  nb_x (RPt) : [ 103, 103, 203, 203, 103 ];
	attribute<int32>  nb_y (RPt) : [   0, 100, 100,   0,   0 ];
	attribute<ipoint> neighbour (Solo, polygon) := points2polygon(point_xy(nb_x, nb_y, ipoint), seq, ord);

	// an inner square [20,60] x [20,60], to be cut out of the square
	attribute<int32>  in_x (RPt) : [ 20, 20, 60, 60, 20 ];
	attribute<int32>  in_y (RPt) : [ 20, 60, 60, 20, 20 ];
	attribute<ipoint> inner (Solo, polygon) := points2polygon(point_xy(in_x, in_y, ipoint), seq, ord);
}

1. cleaning a source

dms_polygon is the clean-up: the sweep reads the operand, resolves whatever it finds and writes a valid polygon.

attribute<ipoint> cleaned (Solo, polygon) := dms_polygon(geometry);

Applied to the bow tie this is the even-odd reading: the ring crosses itself in the middle, so the result is the two triangles, area 2 x 2500 = 5000, and not the 10000 of the square that its four corners span.

attribute<ipoint>  cleaned_bowtie (Solo, polygon) := dms_polygon(bowtie);
parameter<float64> a_bowtie := sum(float64(area(cleaned_bowtie)));   // 5000

To clean and dissolve a whole layer in one step, dms_union_polygon reads every element the same way before it merges them.

The same call on geometry that is already valid returns it unchanged apart from vertices that lie on a straight line between their neighbours, which are dropped.

2. using an invalid polygon as an operand directly

A cleaning step is not required: the invalid ring may be an operand of the operation itself. The left half plane cuts the left triangle out of the bow tie.

attribute<int32>  half_x (RPt) : [ 0,   0, 50, 50, 0 ];
attribute<int32>  half_y (RPt) : [ 0, 100, 100,  0, 0 ];
attribute<ipoint> left_half (Solo, polygon) := points2polygon(point_xy(half_x, half_y, ipoint), seq, ord);

attribute<ipoint>  left_triangle (Solo, polygon) := dms_intersect(bowtie, left_half);
parameter<float64> a_left := sum(float64(area(left_triangle)));      // 2500

3. closing slivers with a grid size

Neighbouring parcels from different sources rarely share their boundary exactly. Without a tolerance the union of the square and its neighbour keeps the three-unit gap between them, so the result is two shells; with a grid of 10 the neighbour’s edge snaps onto the square’s and the two become one rectangle. The sequence_element_count tells the two cases apart: 5 + 5 + 1 points for two polygons strung together, 5 points for one closed ring.

attribute<ipoint> apart  (Solo, polygon) := dms_union(square, neighbour);        // two shells, 11 points
attribute<ipoint> merged (Solo, polygon) := dms_union(square, neighbour, 10.0);  // one ring, 5 points

parameter<uint32>  n_apart  := sum(sequence_element_count(apart));    // 11
parameter<uint32>  n_merged := sum(sequence_element_count(merged));   // 5
parameter<float64> a_merged := sum(float64(area(merged)));            // 20000

Choose the grid from what the coordinates mean, not from the size of the gap: it is the resolution below which the model does not distinguish two positions. On integer coordinates it is a whole number of coordinate units; on float coordinates any positive number. One grid serves the whole attribute, so that two neighbours snap a shared boundary the same way.

4. cutting a hole

A difference whose second argument lies inside the first cuts a hole, written in the corridor layout of a multi polygon: the shell of 5 points, the hole of 5 points, and one point back to the shell’s start.

attribute<ipoint>  holed (Solo, polygon) := dms_difference(square, inner);
parameter<float64> a_holed := sum(float64(area(holed)));              // 8400 = 10000 - 1600
parameter<uint32>  n_holed := sum(sequence_element_count(holed));     // 11

checking a result

Three cheap checks say most of what there is to know about polygon output, and they work on any of the families:

  • Area on the result, compared with the area you can compute by hand, or with the same operation on another family.
  • sum(float64(area(dms_xor(a, b)))) == 0.0 as an exact equality test between two polygon attributes: the symmetric difference of two equal geometries is empty.
  • sequence_element_count for the layout: an unexpected point count means an unexpected number of rings.

further examples

The configurations that pin these operators ship with GeoDMS in examples\testcases and can be run with GeoDmsRun or opened in the GUI:

  • oper_dms_overlay.dms - the areas of the four operations, the identities A * A == A, A + A == A, A ^ A and A - A empty, agreement with the bp_ and geos_ results, the bow tie, an unclosed ring, a hole, and the grid argument.
  • oper_dms_overlay_types.dms - the same operations on all six polygon value types, on negative coordinates and on fractional float coordinates.
  • oper_dms_overlay_topology.dms - the multi polygon layout: several shells, a hole, an island in a lake, rings that touch along an edge and rings that touch in one point, a ring walked twice, and a Parameter as operand.
  • oper_dms_overlay_neg1.dms to oper_dms_overlay_neg5.dms - what is refused: a grid of zero, a fractional grid on integer coordinates, a grid given per element, a grid too fine for the extent, and two operands in different coordinate systems.

see also