Border polygons

Configuration examples Border polygons

introduction

This script presents examples on how to determine the polygons of a polygon data item A that are located at the border of polygon data item B. An example is to determine the municipalities at the border of the Netherlands, see picture:

img

three approahes

To achieve this goal, the configuration example below presents three approaches. All approaches use a municpality base layer (of 2022) of the Netherlands. From this base map a layer of the Netherlands is made with the Union_polygon (dissolve) function.

1. The inflate_approach first uses the layer for the Netherlands and a large square around this layer to make a layer of the area around the Netherlands (outside_nl). Next the municipality base layer is inflated with 10 meter, using the bg_buffer_multi_polygon function. With the overlay_polygon function the intersections are calculated of the inflated municipalities with the area outside the Netherlands. Municipalities that have intersections are presented as municipalities at the border (gemeente_border).
2. The polygon_connectivity_approach uses the polygon_connectivity function to find the connected municipality polygons of all muncipalities in the Netherlands with the area around the Netherlands. Municipalities that have connections with the area around the Netherlands are presented as municipalities at the border (gemeente_border).
3. The outer_points_approach uses the sequence2points function to genearte the point set of the municipalities and of the Netherlands. If at least one of the points of a municipality also occurs in the point set of the Nethelands, the municipalty is presented as municipality at the border (gemeente_border). Advise: use integer points for this approach to prevent rounding issues.

The first approach is time/memory consuming, the last approach is the fastest. In specific cases there can be issues with this third approach, in these cases we advice the second approach.

example

container polygon_border
{
	container units
	{
		unit<float32> m := baseunit('m', float32), label = "meter";
	}

	container Geography: Using = "units"
	{
		#include<wmts_layer.dms>

		unit<fpoint>  point_rd_base : //baseunit('m', fpoint)
			DialogData = "wmts_layer" 
		,	Format     = "EPSG:28992";
		unit<fpoint> point_rd := range(point_rd_base, point(0[m],300000[m]), point(280000[m],625000[m]));
	}

	unit<uint32> cbs_gemeente_2022_with_water:
		StorageName     = "%SourceDataDir%/CBS/2022/gemeente_2022_v1.shp"
	,	StorageReadOnly = "True"
	,	StorageType     = "gdal.vect"
	{
		attribute<Geography/point_rd> geometry (polygon);
	}

	unit<uint32> cbs_gemeente_2022 := select_with_attr_by_org_rel(cbs_gemeente_2022_with_water, cbs_gemeente_2022_with_water/H2O == 'NEE')
	{
		attribute<ipoint> geometry_ipoint (polygon) := ipolygon(geometry);
	}

	parameter<ipoint>             nl_ipoint (polygon) := union_polygon(cbs_gemeente_2022/geometry_ipoint);
	parameter<Geography/point_rd> nl        (polygon) := nl_ipoint[Geography/point_rd];

	container inflate_approach
	{
		unit<uint32> outside_nl : nrofrows = 1
		{
			container square 
			{
				attribute<Geography/point_rd> geometry (polygon, outside_nl) := points2sequence(pointset/coord, const(0, pointset, outside_nl), id(pointset));

				parameter<units/m> left   := 0[units/m];
				parameter<units/m> right  := 350000[units/m];
				parameter<units/m> top    := 700000[units/m];
				parameter<units/m> bottom := 250000[units/m];

				parameter<Geography/point_rd> left_top     := point_xy( left,    top, Geography/point_rd);
				parameter<Geography/point_rd> right_top    := point_xy(right,    top, Geography/point_rd);
				parameter<Geography/point_rd> right_bottom := point_xy(right, bottom, Geography/point_rd);
				parameter<Geography/point_rd> left_bottom  := point_xy( left, bottom, Geography/point_rd);

				unit<uint32> pointset : nrofrows = 5
				{
					attribute<Geography/point_rd> coord := union_data(.,left_top, right_top, right_bottom, left_bottom, left_top);
				}
			}
			attribute<ipoint>             geometry_ipoint (polygon) := (ipolygon(square/geometry) - union_data(., nl_ipoint)); 
			attribute<Geography/point_rd> geometry        (polygon) := geometry_ipoint[Geography/point_rd]; 
		}

		attribute<Geography/point_rd> cbs_gemeente_2022_inflated (polygon, cbs_gemeente_2022) := bg_buffer_multi_polygon(cbs_gemeente_2022/geometry, 10.0, 16b);

		unit<uint32> intersect := overlay_polygon(ipolygon(cbs_gemeente_2022_inflated), ipolygon(outside_nl/geometry));

		unit<uint32> gemeente_border := unique(intersect/first_rel)
		{
			attribute<string>             GM_NAAM            := cbs_gemeente_2022/GM_NAAM[values];
			attribute<Geography/point_rd> geometry (polygon) := cbs_gemeente_2022/geometry[values];
		}
	}

	container polygon_connectivity_approach
	{
		unit<uint32> buiten_nl_plus_gemeente := union_unit(void, cbs_gemeente_2022)
		{
			attribute<string> GM_NAAM                := union_data(., 'buiten NL', cbs_gemeente_2022/GM_NAAM);
			attribute<ipoint> geometry_ipoint (poly) := union_data(., inflate_approach/outside_nl/geometry_ipoint ,cbs_gemeente_2022/geometry_ipoint);
		}

		unit<uint32> points_nl := polygon_connectivity(buiten_nl_plus_gemeente/geometry_ipoint);

		unit<uint32> connected_to_dutch_order := select_with_attr_by_org_rel(points_nl, points_nl/F1 == 0);

		unit<uint32> gemeente_border:= unique(connected_to_dutch_order/F2)
		{
			attribute<string> GM_NAAM                        := buiten_nl_plus_gemeente/GM_NAAM[values];
			attribute<ipoint> geometry_ipoint      (polygon) := buiten_nl_plus_gemeente/geometry_ipoint[values];
			attribute<Geography/point_rd> geometry (polygon) := geometry_ipoint[Geography/point_rd];
		}
	}

	container outer_points_approach
	{
		unit<uint32> points_gemeente := sequence2points(cbs_gemeente_2022/geometry_ipoint)
		{
			attribute<bool> isBorderPoint := isDefined(rlookup(point, points_nl/point));
		}

		unit<uint32> points_nl := sequence2points(nl_ipoint);

		attribute<bool> is_gemeente_at_border (cbs_gemeente_2022) := any(points_gemeente/isBorderPoint, points_gemeente/sequenceNr);

		unit<uint32> gemeente_border := select_with_org_rel(is_gemeente_at_border)
		{
			attribute<string>             GM_NAAM            := cbs_gemeente_2022/GM_NAAM[org_rel];
			attribute<Geography/point_rd> geometry (polygon) := cbs_gemeente_2022/geometry[org_rel];
		}
	}
}