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 point_rd_wms : format = "EPSG:28992"; unit point_rd := range(point_rd_wms, point(300000f,0f), point(625000f,280000f));

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

container grid { unit 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 m := baseunit('m', float32); unit m2 := m * m; unit 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];
  }    } }