Border polygons
configuration-examples Border polygons
introduction
This script presents examples on how to determine the polygons of a polygon data-item A located at the border of polygon data-item B.
An example is to determine the municipalities at the border of the Netherlands:
three approaches
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.
- The inflate_approach starts with 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 between the inflated municipalities and the area outside the Netherlands. Municipalities that have intersections are presented as municipalities at the border (gemeente_border).
- The polygon_connectivity_approach uses the polygon_connectivity function to find the connected municipality polygons of all municipalities 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).
- The outer_points_approach uses the sequence2points function to generate 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 Netherlands, the municipality is presented as municipality at the border (gemeente_border).
Advise: use integer points for these approaches, as most geometric functions require integer coordinates.
The first approach is time/memory consuming, the third approach is the fastest.
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];
}
}
}