-
Notifications
You must be signed in to change notification settings - Fork 1
Poly to grid to (simplified) polygon
Jip Claassens edited this page Nov 4, 2025
·
14 revisions
Configuration examples Poly to grid to (simplified) polygon
This script can be used for two use cases:
-
To simplify a polygon dataset while maintaining the topology. This means linestrings used across multiple polygons are simplified in the same manner, resulting in a polygon set with the same topology as the source dataset (the Bg_simplify_polygon and bg_simplify_multi_polygon functions do not guarantee maintaining the topology).
-
To make a set of (simplified) polygons of adjacent grid cells with the same values.
poly_2_grid_2_simplified_poly.zip
In the examples container of this script, you find different options:
- simplify_gridsize: simplifies the source data with the size of the grid. In this container you find options with grid sizes of 500 and 100 meters as spoint, and 10 meters as ipoint.
- simplify_m10: simplifies the source data with as simplify factor 10 meter. In this container you find options with grid sizes of 500 and 100 meters as spoint, and 10 meters as ipoint.
- simplify_m0: simplifies the source data with as simplify factor 0 meter. This means only intermediate points in straight lines are removed, as they do not affect the shape of the feature. In this container you find options with grid sizes of 500 and 100 meters as spoint, and 10 meters as ipoint.
The resulting polygons in the rijksdriehoek coordinates can be requested with the item:
examples/simplify_size_/d_size_i_s_point/SimplifiedPolygons/Results/geometry.
- poly2grid
- points2polygon
- impedance_table
- connected_parts
- point_xy
- collect_by_cond
- union_data
- select
- combine_data
// This is the template that uses a grid attribute to make a (simplified) polygonset of connected grid cells with the same attribute value
// In this template such a set of connected cells is called a zone;
template grid2poly_ipoint
{
// begin case parameters
unit<ipoint> domain_grid;
unit<uint32> Zone;
attribute<Zone> zone_rel (domain_grid);
parameter<float64> simplifyFactor;
// end case parameters
parameter<Zone> NoZone := #Zone + lowerbound(Zone); // outside study area indicator
// Segments := horizontale raster-randjes met ongelijke boven en onder zones + verticale raster-randjes met ongelijke linker en rechter zones
unit<int32> rows := range(int32, pointrow(lowerbound(domain_grid)), pointrow(upperbound(domain_grid)));
unit<int32> cols := range(int32, pointcol(lowerbound(domain_grid)), pointcol(upperbound(domain_grid)));
unit<int32> HorizontalLine := range(int32, lowerbound(rows), upperbound(rows)+1[rows]);
unit<int32> VerticalLine := range(int32, lowerbound(cols), upperbound(cols)+1[cols]);
unit<uint32> Point := combine(VerticalLine, HorizontalLine)
{
attribute<geography/point_rd> geometry := domain_grid_rel[geography/point_rd];
attribute<domain_grid> domain_grid_rel := point_xy(first_rel, second_rel, domain_grid); // the RasterCell NorthEast of this Point
attribute<bool> isBoundaryPoint := has_any(BoundarySegment/F1) || has_any(BoundarySegment/F2);
}
unit<uint32> BoundaryPoint := select_with_org_rel(Point/isBoundaryPoint)
{
attribute<Point> Point_rel := org_rel;
attribute<geography/point_rd> geometry := org_rel->geometry;
attribute<uint32> NrBoundarySegments := pcount(BoundarySegment/bp_F1) + pcount(BoundarySegment/bp_F2);
attribute<bool> isJunction := NrBoundarySegments != 2; // Junction := raster-hoekpunt waarop meer dan 2 segmenten uit komen (1 kan niet voor komen)
attribute<BoundarySegment> FirstLink := MakeDefined(first(ID(BoundarySegment), BoundarySegment/bp_F1), first(ID(BoundarySegment), BoundarySegment/bp_F2));
attribute<BoundarySegment> Last_Link := MakeDefined(last (ID(BoundarySegment), BoundarySegment/bp_F2), last (ID(BoundarySegment), BoundarySegment/bp_F1));
}
unit<uint32> HorizontalSegment := combine(cols, HorizontalLine)
{
attribute<domain_grid> North := point_xy(first_rel, second_rel-1[HorizontalLine], domain_grid);
attribute<domain_grid> South := point_xy(first_rel, second_rel , domain_grid);
attribute<bool> isBoundarySegment := MakeDefined(south->zone_rel, NoZone) != MakeDefined(North->zone_rel, NoZone);
}
unit<uint32> VerticalSegment := combine(VerticalLine, rows)
{
attribute<domain_grid> West := point_xy(first_rel-1[VerticalLine], second_rel, domain_grid);
attribute<domain_grid> East := point_xy(first_rel , second_rel, domain_grid);
attribute<bool> isBoundarySegment := MakeDefined(West->zone_rel, max(zone_rel)+1[zone]) != MakeDefined(East->zone_rel, max(zone_rel)+1[zone]);
}
unit<uint32> HorizontalBoundarySegment := select(HorizontalSegment/isBoundarySegment)
{
attribute<cols> first_rel := collect_by_cond(., HorizontalSegment/first_rel);
attribute<HorizontalLine> second_rel := collect_by_cond(., HorizontalSegment/second_rel);
attribute<domain_grid> North := collect_by_cond(., HorizontalSegment/North);
attribute<domain_grid> South := collect_by_cond(., HorizontalSegment/South);
attribute<Point> F1 := combine_data(Point, first_rel [VerticalLine], second_rel);
attribute<Point> F2 := combine_data(Point, (first_rel+1[cols])[VerticalLine], second_rel);
}
unit<uint32> VerticalBoundarySegment := select(verticalSegment/isBoundarySegment)
{
attribute<VerticalLine> first_rel := collect_by_cond(., VerticalSegment/first_rel);
attribute<rows> second_rel := collect_by_cond(., VerticalSegment/second_rel);
attribute<domain_grid> East := collect_by_cond(., VerticalSegment/East);
attribute<domain_grid> West := collect_by_cond(., VerticalSegment/West);
attribute<Point> F1 := combine_data(Point, first_rel, second_rel [HorizontalLine]);
attribute<Point> F2 := combine_data(Point, first_rel, (second_rel+1[rows])[HorizontalLine]);
}
unit<uint32> BoundarySegment := union_unit(HorizontalBoundarySegment, VerticalBoundarySegment)
{
attribute<Direction> Direction_rel := union_data(., const(Direction/V/South, HorizontalBoundarySegment), const(Direction/V/East, VerticalBoundarySegment));
unit<uint32> BS_Point := union_unit(BoundarySegment, BoundarySegment)
{
attribute<geography/point_rd> geometry := union_data(BS_Point, F1->geometry, F2->geometry);
attribute<BoundarySegment> BoundarySegment_rel := union_data(BS_Point, ID(BoundarySegment), ID(BoundarySegment));
}
attribute<geography/point_rd> geometry(arc) := points2sequence(BS_Point/geometry, BS_Point/BoundarySegment_rel);
attribute<uint32> NrIntermediatePoints := pcount(IntermediatePoint/Ra) + pcount(IntermediatePoint/Rb);
attribute<domain_grid> left := union_data(., HorizontalBoundarySegment/North, VerticalBoundarySegment/East);
attribute<domain_grid> right := union_data(., HorizontalBoundarySegment/South, VerticalBoundarySegment/West);
attribute<Direction> Dir_rel := union_data(., const(Direction/V/East ,HorizontalBoundarySegment), const(Direction/V/South ,VerticalBoundarySegment));
attribute<Point> F1 := union_data(., HorizontalBoundarySegment/F1, VerticalBoundarySegment/F1);
attribute<Point> F2 := union_data(., HorizontalBoundarySegment/F2, VerticalBoundarySegment/F2);
attribute<BoundaryPoint> bp_F1 := rlookup(F1, BoundaryPoint/Point_rel);
attribute<BoundaryPoint> bp_F2 := rlookup(F2, BoundaryPoint/Point_rel);
attribute<bool> f1IsJunction := bp_F1->isJunction;
attribute<bool> f2IsJunction := bp_F2->isJunction;
attribute<Zone> left_Zone_rel := zone_rel[left];
attribute<Zone> rightZone_rel := zone_rel[right];
attribute<JFSection> JFSection_rel := JFSection/part_rel;
}
unit<uint32> IntermediatePoint := select_with_org_rel(not(BoundaryPoint/isJunction))
{
attribute<BoundaryPoint> BoundaryPoint_rel := org_rel;
attribute<BoundarySegment> Ra := collect_by_cond(IntermediatePoint, BoundaryPoint/FirstLink);
attribute<BoundarySegment> Rb := collect_by_cond(IntermediatePoint, BoundaryPoint/Last_Link);
attribute<bool> IsOke := IsDefined(Ra) && IsDefined(Rb) && Ra != Rb;
}
unit<uint32> Junction := select_with_org_rel(BoundaryPoint/isJunction)
{
attribute<geography/point_rd> geometry := collect_by_cond(Junction, BoundaryPoint/geometry);
attribute<.> per_node(BoundaryPoint) := invert(org_rel);
}
unit<uint32> FromJunctionLink := select_with_org_rel(BoundarySegment/f1IsJunction)
{
attribute<BoundarySegment> BoundarySegment_rel := org_rel;
attribute<BoundaryPoint> F1 := collect_by_cond(FromJunctionLink, BoundarySegment/bp_F1);
attribute<BoundaryPoint> F2 := collect_by_cond(FromJunctionLink, BoundarySegment/bp_F2);
attribute<Direction> Dir_rel := collect_by_cond(FromJunctionLink, BoundarySegment/Dir_rel);
attribute<Junction> Junction_rel := Junction/per_node[F1];
attribute<Zone> left_Zone_rel := collect_by_cond(FromJunctionLink, BoundarySegment/left_Zone_rel);
attribute<Zone> rightZone_rel := collect_by_cond(FromJunctionLink, BoundarySegment/rightZone_rel);
}
unit<uint32> ToJunctionLink := select_with_org_rel(BoundarySegment/f2IsJunction)
{
attribute<BoundarySegment> BoundarySegment_rel := org_rel;
attribute<BoundaryPoint> F1 := collect_by_cond(ToJunctionLink, BoundarySegment/bp_F1);
attribute<BoundaryPoint> F2 := collect_by_cond(ToJunctionLink, BoundarySegment/bp_F2);
attribute<Direction> Org_dir_rel := collect_by_cond(ToJunctionLink, BoundarySegment/Dir_rel);
attribute<Direction> Dir_rel := Org_dir_rel->Direction/Reverse;
attribute<Junction> Junction_rel := Junction/per_node[F2];
attribute<Zone> left_Zone_rel := collect_by_cond(ToJunctionLink, BoundarySegment/left_Zone_rel);
attribute<Zone> rightZone_rel := collect_by_cond(ToJunctionLink, BoundarySegment/rightZone_rel);
}
unit<uint32> JunctionDirection := combine_unit(Junction, Direction);
unit<uint32> JunctionLink := union_unit(FromJunctionLink, ToJunctionLink)
{
attribute<geography/point_rd> Junction_geometry := Junction_rel->Geometry;
attribute<BoundarySegment> BoundarySegment_rel := union_data(., FromJunctionLink/BoundarySegment_rel, ToJunctionLink/BoundarySegment_rel);
unit<uint32> JL_Point := union_unit(JunctionLink, JunctionLink)
{
attribute<geography/point_rd> geometry := union_data(JL_Point, F1->geometry, F2->geometry);
attribute<JunctionLink> JunctionLink_rel := union_data(JL_Point, ID(JunctionLink), ID(JunctionLink));
}
attribute<geography/point_rd> geometry(arc) := points2sequence(JL_Point/geometry, JL_Point/JunctionLink_rel);
attribute<BoundaryPoint> F1 := union_data(., FromJunctionLink/F1, ToJunctionLink/F2);
attribute<BoundaryPoint> F2 := union_data(., FromJunctionLink/F2, ToJunctionLink/F1);
attribute<Direction> Dir_rel := union_data(., FromJunctionLink/Dir_rel, ToJunctionLink/Dir_rel);
attribute<Direction> NextDir_rel := Dir_rel->Next;
attribute<Direction> NextNextDir_rel := Dir_rel->Reverse;
attribute<Junction> Junction_rel := rlookup(F1, Junction/org_rel);
attribute<Zone> left_Zone_rel := union_data(., FromJunctionLink/left_Zone_rel, ToJunctionLink/rightZone_rel);
attribute<Zone> rightZone_rel := union_data(., FromJunctionLink/rightZone_rel, ToJunctionLink/left_Zone_rel );
attribute<JunctionDirection> JunctionDirection_rel := combine_data(JunctionDirection, Junction_rel, Dir_rel);
attribute<JunctionDirection> NextJunctionDirection_rel := combine_data(JunctionDirection, Junction_rel, NextDir_rel);
attribute<JunctionDirection> NextNextJunctionDirection_rel := combine_data(JunctionDirection, Junction_rel, NextNextDir_rel);
attribute<JunctionLink> n1_JL := rlookup(NextJunctionDirection_rel, JunctionDirection_rel);
attribute<JunctionLink> n2_JL := rlookup(NextNextJunctionDirection_rel, JunctionDirection_rel);
attribute<JunctionLink> next_JL := MakeDefined(n1_JL, n2_JL);
attribute<JFLS2> outJFLS2 := invert(JFLS2/JL_rel);
attribute<.> other_rel := outJFLS2->Other_rel;
attribute<.> next_rel := other_rel->next_JL;
}
// Junction Free Section := clusters van BoundarySegments
unit<uint32> JFSection := connected_parts(intermediatePoint/Ra, intermediatePoint/Rb)
{
attribute<uint32> nr_BoundarySegments := pcount(PartNr);
attribute<uint8> NrEnds := sum_uint8(2 - BoundarySegment/NrIntermediatePoints, PartNr);
attribute<bool> isLineString := NrEnds == 2b; // it could also be a ring without any Junction
attribute<.> per_IntermediatePoint(IntermediatePoint) := MakeDefined(BoundarySegment/JFSection_rel[IntermediatePoint/Ra], BoundarySegment/JFSection_rel[IntermediatePoint/Rb]);
}
unit<uint32> JFRing := select_with_org_rel(not(JFSection/isLineString))
{
attribute<JFSection> JFSection_rel := org_rel;
attribute<BoundarySegment> JFSection_firstBS_rel (JFSection) := min_index(ID(BoundarySegment), JFSection/PartNr);
attribute<BoundarySegment> firstBoundarySegment_rel := lookup(JFSection_rel, JFSection_firstBS_rel);
attribute<BoundaryPoint> firstBoundaryPoint_rel := firstBoundarySegment_rel->bp_F1;
attribute<bool> IsOke := firstBoundarySegment_rel->Direction_rel == Direction/V/South; // follows from the order or Points and Segments
attribute<bool> BoundarySegment_IsNotFirstRingSegment (BoundarySegment) := !IsDefined(invert(firstBoundarySegment_rel));
// abuse the fact that the first segment of a ring is always North to South and therefore right is outside en left is inside.
attribute<Zone> left_Zone_rel := firstBoundarySegment_rel->rightZone_rel;
attribute<Zone> rightZone_rel := firstBoundarySegment_rel->left_Zone_rel;
unit<uint32> ReducedBoundarySegment := select(not(BoundarySegment/JFSection_rel->isLineString) && BoundarySegment_IsNotFirstRingSegment)
{
attribute<BoundaryPoint> F1 := collect_by_cond(., BoundarySegment/bp_F1);
attribute<BoundaryPoint> F2 := collect_by_cond(., BoundarySegment/bp_F2);
}
attribute<.> per_IntermediatePoint (IntermediatePoint) := invert(JFSection_rel)[JFSection/per_IntermediatePoint];
attribute<.> per_BoundaryPoint (BoundaryPoint) := lookup(invert(intermediatePoint/org_rel), per_IntermediatePoint);
// as the first segment has been removed, the remaining segments form a clock-wise path from the first BoundaryPoint to the other side of the removed segment
attribute<uint32> dist_from_first_BP(BoundaryPoint) := impedance_table('bidirectional;startPoint(Node_rel);node:TraceBack'
, const(1, ReducedBoundarySegment), ReducedBoundarySegment/F1, ReducedBoundarySegment/F2, firstBoundaryPoint_rel);
unit<uint32> Point := select(IsDefined(dist_from_first_BP))
{
attribute<geography/point_rd> geometry := collect_by_cond(Point, BoundaryPoint/Geometry);
attribute<JFRing> sequence_rel := collect_by_cond(Point, per_BoundaryPoint);
attribute<uint32> ordinal := collect_by_cond(Point, dist_from_first_BP);
}
attribute<geography/point_rd> geometry (poly) := points2polygon(Point/geometry, Point/sequence_rel, Point/ordinal);
attribute<units/meter2> area := area(geometry, units/meter2);
attribute<geography/point_rd> SimplifiedRing (poly) := bg_simplify_polygon(geometry, simplifyFactor);
attribute<units/meter2> simplified_area := area(SimplifiedRing, units/meter2);
}
unit<uint32> NonEmptyJFRing := select(abs(JFRing/simplified_area) > 100d)
{
attribute<geography/point_rd> geometry (poly) := collect_by_cond(NonEmptyJFRing, JFRing/geometry);
attribute<geography/point_rd> SimplifiedRing (poly) := collect_by_cond(NonEmptyJFRing, JFRing/SimplifiedRing);
attribute<Zone> left_Zone_rel := collect_by_cond(., JFRing/left_Zone_rel);
attribute<Zone> rightZone_rel := collect_by_cond(., JFRing/rightZone_rel);
unit<uint32> Point := sequence2points(SimplifiedRing);
attribute<geography/point_rd> SimplifiedReversedRing(poly) := points2polygon(Point/Point, Point/sequence_rel, (pcount(Point/sequence_rel) - 1)[Point/sequence_rel] - Point/ordinal);
}
unit<uint32> JFLineString := select_with_org_rel(JFSection/isLineString) // JunctionFree LineString
{
attribute<JunctionLink> firstJL_rel := collect_by_cond(., first(ID(JunctionLink), JunctionLink/BoundarySegment_rel->JFSection_rel));
attribute<JunctionLink> last_JL_rel := collect_by_cond(., last (ID(JunctionLink), JunctionLink/BoundarySegment_rel->JFSection_rel));
attribute<BoundaryPoint> firstBP_rel := firstJL_rel->F1;
attribute<BoundaryPoint> secondBP_rel := firstJL_rel->F2;
attribute<BoundaryPoint> last_BP_rel := last_JL_rel->F2; // just before the end of the last segment !!!
// JFSection/Linestring := d.m.v. impedance_single van alle begin JFSection naar eind JFSection routes t.b.v. points2sequence
attribute<BoundarySegment> last_BoundarySegment := last_JL_rel->BoundarySegment_rel;
attribute<BoundarySegment> firstBoundarySegment := firstJL_rel->BoundarySegment_rel;
unit<uint32> ReducedBoundarySegment := select(BoundarySegment/JFSection_rel->isLineString && !IsDefined(invert(last_BoundarySegment)) && !IsDefined(invert(firstBoundarySegment)))
{
attribute<BoundaryPoint> F1 := collect_by_cond(., BoundarySegment/bp_F1);
attribute<BoundaryPoint> F2 := collect_by_cond(., BoundarySegment/bp_F2);
}
attribute<uint32> dist_from_first_BP(BoundaryPoint) := impedance_table('bidirectional;startPoint(Node_rel);node:TraceBack'
, const(1, ReducedBoundarySegment), ReducedBoundarySegment/F1, ReducedBoundarySegment/F2, secondBP_rel);
unit<uint32> Point := select_with_org_rel(IsDefined(dist_from_first_BP))
{
attribute<BoundaryPoint> BoundaryPoint_rel := org_rel;
attribute<IntermediatePoint> IntermediatePoint_rel := invert(IntermediatePoint/BoundaryPoint_rel)[BoundaryPoint_rel];
attribute<BoundarySegment> BoundarySegmentA_rel := IntermediatePoint_rel->Ra;
attribute<geography/point_rd> geometry := collect_by_cond(Point, BoundaryPoint/Geometry);
attribute<JFLineString> sequence_rel := invert(JFLineString/org_rel)[BoundarySegmentA_rel->JFSection_rel];
attribute<uint32> ordinal := collect_by_cond(Point, dist_from_first_BP);
}
attribute<uint32> NrPoints := pcount(Point/sequence_rel) > 0 ? max(Point/ordinal, Point/sequence_rel) + 1 : 0;
unit<uint32> PointwithConnectors := union_unit(JFLineString, Point, JFLineString)
{
attribute<geography/point_rd> geometry := union_data(., firstJL_rel->Junction_geometry, Point/geometry, last_JL_rel->Junction_geometry);
attribute<JFLineString> sequence_rel := union_data(., id(JFLineString), Point/sequence_rel, id(JFLineString));
attribute<uint32> ordinal := union_data(., const(0, JFLineString), Point/ordinal+1, NrPoints+1);
}
attribute<geography/point_rd> LineString (arc) := points2sequence(Point/geometry, Point/sequence_rel, Point/ordinal);
attribute<geography/point_rd> LineStringWithConnectors (arc) := points2sequence(PointwithConnectors/geometry, PointwithConnectors/sequence_rel, PointwithConnectors/ordinal);
attribute<geography/point_rd> SimplifiedLineString (arc) := bg_simplify_linestring(LineStringWithConnectors, simplifyFactor);
unit<uint32> SimplifiedPoints := sequence2points(SimplifiedLineString);
attribute<uint32> SimplifiedPointOrdinalMax := max(SimplifiedPoints/ordinal, SimplifiedPoints/Sequence_rel);
attribute<uint32> NrSimplifiedPoints := SimplifiedPointOrdinalMax+1;
attribute<SimplifiedPoints> SimplifiedPointBase := value(cumulate(NrSimplifiedPoints) - NrSimplifiedPoints, SimplifiedPoints);
}
unit<uint32> JFLS2 := union_unit(JFLineString, JFLineString)
{
attribute<JunctionLink> JL_rel := union_data(., JFLineString/firstJL_rel, JFLineString/last_JL_rel);
attribute<uint32> NrSimplifiedPoints := union_data(., JFLineString/NrSimplifiedPoints, JFLineString/NrSimplifiedPoints);
attribute<JunctionLink> Other_rel := union_data(., JFLineString/last_JL_rel, JFLineString/firstJL_rel);
attribute<JFLS2_Point> JFLS2_Point_Base := cumulate(NrSimplifiedPoints) - NrSimplifiedPoints;
}
unit<uint32> JFLS2_Point := range(uint32, 0, sum(JFLS2/NrSimplifiedPoints))
{
attribute<JFLS2> JFLS2_rel := classify(ID(.), JFLS2/JFLS2_Point_Base);
attribute<uint32> ordinal := ID(.) - JFLS2_rel;
attribute<JFLineString/SimplifiedPoints> SimplifiedPoint_rel := union_data(.
, JFLineString/SimplifiedPointBase[JFLineString/SimplifiedPoints/Sequence_rel] + JFLineString/SimplifiedPoints/ordinal
, (JFLineString/SimplifiedPointBase+JFLineString/SimplifiedPointOrdinalMax)[JFLineString/SimplifiedPoints/Sequence_rel] - JFLineString/SimplifiedPoints/ordinal);
attribute<geography/point_rd> geometry := SimplifiedPoint_rel->Point;
}
unit<uint32> ZoneRing := connected_parts(ID(JunctionLink), JunctionLink/next_rel)
{
attribute<.> per_JunctionLink (JunctionLink) := PartNr;
attribute<JunctionLink> index (JunctionLink) := index(PartNr);
attribute<JunctionLink> order (JunctionLink) := invert(index);
attribute<Zone> FirstZone_rel := first(JunctionLink/rightZone_rel, per_JunctionLink);
attribute<Zone> FirstZone_per_JunctionLink (JunctionLink) := per_JunctionLink->FirstZone_rel;
attribute<bool> CheckUniqueZonePerRing (JunctionLink) := or(JunctionLink/rightZone_rel == FirstZone_per_JunctionLink, not(or(IsDefined(JunctionLink/rightZone_rel), IsDefined(FirstZone_per_JunctionLink))));
attribute<JunctionLink> FirstJunctionLink_rel := min(order, per_JunctionLink);
attribute<uint32> nrJunctions := pcount(per_JunctionLink);
attribute<ZoneRingJunction> ZoneRingJunctionBase := cumulate(nrJunctions) - nrJunctions;
attribute<ZoneRingPoint> ZoneRingPointBase := ZoneRingJunctionBase->ZoneRingPointBase;
attribute<geography/point_rd> geometry(poly) := points2polygon(ClosedZoneRingPoint/geometry, ClosedZoneRingPoint/ZoneRing_rel, ClosedZoneRingPoint/Ring_ordinal);
attribute<units/meter2> area := area(geometry, units/meter2);
}
unit<uint32> ZoneRingJunction := range(uint32, 0, sum(ZoneRing/nrJunctions))
{
attribute<ZoneRing> ZoneRing_rel := classify(ID(.), ZoneRing/ZoneRingJunctionBase);
attribute<uint32> ordinal := id(.) - ZoneRing_rel->ZoneRingJunctionBase;
attribute<JunctionLink> OrderedJunctionLink_rel := rlookup(ZoneRing_rel->FirstJunctionLink_rel + ordinal, ZoneRing/order);
attribute<JunctionLink> NextJunctionLink_rel := OrderedJunctionLink_rel->next_rel;
attribute<.> per_JunctionLink(JunctionLink) := invert(OrderedJunctionLink_rel);
attribute<.> NextOrderedZoneRing_rel := per_JunctionLink[NextJunctionLink_rel];
attribute<uint32> dist_from_first_JL := impedance_table('directed;startPoint(Node_rel);node:TraceBack'
, const(1, ZoneRingJunction), ID(ZoneRingJunction), NextOrderedZoneRing_rel
, ZoneRing/ZoneRingJunctionBase);
attribute<JunctionLink> JunctionLink_rel := OrderedJunctionLink_rel[invert(ZoneRing_rel->ZoneRingJunctionBase + dist_from_first_JL)];
attribute<JFLS2> JFLS2_rel := JunctionLink_rel->outJFLS2;
attribute<uint32> nrPoints := sub_or_null(JFLS2_rel->NrSimplifiedPoints, 1);
attribute<uint32> ZoneRingPointBase := cumulate(nrPoints) - nrPoints;
}
// compose a ring from the JunctionLink geometries, either a JFSL2_Point sequence, the BoundaryPoint of F1, or the BoundaryPoint of F2 for doubleton_links
unit<uint32> ZoneRingPoint := range(uint32, 0, sum(ZoneRingJunction/nrPoints))
{
attribute<ZoneRingJunction> ZoneRingJunction_rel := classify(ID(.), ZoneRingJunction/ZoneRingPointBase);
attribute<ZoneRing> ZoneRing_rel := ZoneRingJunction_rel->ZoneRing_rel;
attribute<uint32> Ring_ordinal := ID(.) - ZoneRing_rel->ZoneRingPointBase;
attribute<uint32> RingJunction_ordinal := ID(.) - ZoneRingJunction_rel->ZoneRingPointBase;
attribute<JunctionLink> JunctionLink_rel := ZoneRingJunction_rel->JunctionLink_rel;
attribute<JFLS2> JFLS2_rel := JunctionLink_rel->outJFLS2;
attribute<JFLS2_Point> JFLS2_Point_rel := JFLS2_rel->JFLS2_Point_Base + RingJunction_ordinal;
attribute<geography/point_rd> geometry := JFLS2_Point_rel->geometry;
}
// add a closing point to each ring equal to the first point of that ring.
unit<uint32> ClosedZoneRingPoint := range(uint32, 0, sum(ZoneRingJunction/nrPoints)+#ZoneRing)
{
attribute<geography/point_rd> geometry := union_data(., ZoneRingPoint/geometry, ZoneRing/ZoneRingPointBase->geometry);
attribute<ZoneRing> ZoneRing_rel := union_data(., ZoneRingPoint/ZoneRing_rel, id(ZoneRing));
attribute<uint32> Ring_ordinal := union_data(., ZoneRingPoint/Ring_ordinal, sum(ZoneRingJunction/nrPoints, ZoneRingJunction/ZoneRing_rel));
attribute<Zone> Zone_rel := ZoneRing_rel->FirstZone_rel;
}
unit<uint32> RingResults := union_unit(ZoneRing, NonEmptyJFRing, NonEmptyJFRing)
{
attribute<geography/point_rd> geometry(poly) := union_data(., ZoneRing/geometry, NonEmptyJFRing/SimplifiedRing, NonEmptyJFRing/SimplifiedReversedRing);
attribute<units/meter2> area := area(geometry, units/meter2);
attribute<Zone> Zone_rel := union_data(., ZoneRing/FirstZone_rel, NonEmptyJFRing/RightZone_rel, NonEmptyJFRing/Left_Zone_rel), IntegrityCheck = "ZoneRing/CheckUniqueZonePerRing";
}
unit<uint32> Results := Zone
{
// Polygonen van dezelfde regio samenvoegen tot multi-polygons, in een afzonderlijke template
// TODO: met afzonderlijke punten doen o.b.v. de te checked aanname dat ze niet overlappen.
attribute<geography/point_rd> geometry(poly) := geos_union_polygon(RingResults/geometry, RingResults/Zone_rel);
// check dat result geen overlappende polygonen bevat
unit<uint32> overlap_check := geos_polygon_connectivity(geometry)
{
attribute<geography/point_rd> overlap (poly) := geometry[first_rel] * geometry[second_rel];
attribute<units/meter2> area := area(overlap, units/meter2);
}
}
}
GeoDMS ©Object Vision BV. Source code distributed under GNU GPL-3. Documentation distributed under CC BY-SA 4.0.