Spatial Overlap

configuration-examples Spatial Overlap

In polygon layers, polygons might overlap. This can cause issues for instance with the point_in_polygon function in the GeoDMS (the sequence of polygons then defines the result, which is often undesired).

It can be useful, before working with a polygon layer, to first find out if in the layer polygons do overlap. For that purpose the following script shows two options.

1. grid approach. This gives a fast indication for substantial overlap. By reducing the grid size, the approach becomes more precise. The grid approach uses the poly2grid function for all polygon geometries in the order of the src file and the reverse order. If these approaches result in different outcomes, there must be overlap. The resulting boolean parameter: anyOverlap indicates if overlap does occur.
2. vector approach. This apporach is more time consuming, but give an exact area of overlap for the overlapping polygons. The polygon_connectivity function is used to find out which polygons are connected. For this set the multiplication function result in the overlap between these polygons. With the area function, the surface of the overlap is calculated.

script

container SpatialOverlap
{
   unit<fpoint>  point_rd_wms : format = "EPSG:28992";
   unit<fpoint>  point_rd     := range(point_rd_wms, point(300000f,0f), point(625000f,280000f)); 

 // the polygon set that might contain overlap is configured here <;
   unit<uint32> src_poly:
      StorageName     = "%SourceDataDir%/CBS/2019/cbs_gem_2019_si.shp"
   ,  StorageType     = "gdal.vect"
   ,  StorageReadOnly = "True"
   {
      attribute<point_rd> geometry (poly);
   }

   container grid
   {
      unit<spoint> domain :=  range(
         gridset(
             point_rd
            ,point(-100f, 100f, point_rd)
            ,point(625000f, 10000f, point_rd)
            ,spoint
         )
            ,point(0s, 0s)
            ,point(3250s, 2700s)
         )
      ,   DialogType = "point_rd";

       // a reverse domain is configured with the same number of rows as the source polygon domain, 
             only with the geometries in reverse order 
       unit<uint32> reverse_poly := range(uint32, 0, #src_poly)
      {
         attribute<point_rd> geometry   (poly)  := src_poly/geometry[#reverse_poly - (id(reverse_poly) + 1)];
         attribute<src_poly> src_poly_rel       := #. - (id(.) + 1);
      }

      // for both the source polygon domain and the reverse domain the poly2grid function is applied
         , the results are made comparable
      attribute<src_poly>     src_poly_rel               (domain) := poly2grid(src_poly/geometry    , domain);
      attribute<reverse_poly> reverse_poly_rel           (domain) := poly2grid(reverse_poly/geometry, domain);
      attribute<src_poly>     reverse_poly_src_poly_rel  (domain) := reverse_poly/src_poly_rel[reverse_poly_rel];
      attribute<bool>         hasOverlap                 (domain) := src_poly_rel <> reverse_poly_src_poly_rel;

      parameter<bool> anyOverlapa := any(hasOverlap);
   }

   container vector
   {
      unit<float32> m  := baseunit('m', float32);
      unit<float32> m2 := m * m;
      unit<ipoint>  point_rd_mm :=
         gridset(
             point_rd
            ,point(0.001f, 0.001f, point_rd)
            ,point(0f    , 0f    , point_rd)
            ,ipoint
         );

       // the  polygon_connectivity results in the set of connected/overlapping polygons 
      unit<uint32> overlap_vector := polygon_connectivity(ipolygon(src_poly/geometry[point_rd_mm]))
      {
         attribute<point_rd> geometry_F1 (poly) := src_poly/geometry[F1];
         attribute<point_rd> geometry_F2 (poly) := src_poly/geometry[F2];
   
         attribute <ipoint>  intersect   (poly) := 
            ipolygon(geometry_F1[point_rd_mm]) * ipolygon(geometry_F2[point_rd_mm]);
         attribute<m2>       area               := area(intersect[point_rd], m2);
      }

      unit<uint32> met_overlap_vector := select_with_org_rel(overlap_vector/area > 0[m2])
      {
         attribute<point_rd> geometry (poly)  := value(overlap_vector/intersect[org_rel], point_rd);
         attribute<m2>       area             := overlap_vector/area[org_rel];
      }
   }
}