[GIS] intersection between two multipolygons yielding anomalous GeometryCollection object full of LineString’s & Polygon’s (trying to get intersect area)

intersectionpolygonruby

Normally when I call the #intersection method on a multipolygon to find the intersection shape on another multipolygon I get a multipolygon returned back to me. Either that or something else I can then call the #area function on.

Well, 150 spatial joins into a list of zoning shapefiles (ie this is a rare case) I'm scanning through my code throws an error because it can't call the #area method on a CAPIGeometryCollection. I am having problems linking qgis with my postgresql db and don't know how to create a kml for a "GeometryCollection" object

Here's what the WKT string output for the intersection shape looks like:

"GEOMETRYCOLLECTION (LINESTRING (1480424.245856002 596629.5658560097, 1480362.875856012 596566.8158560097), LINESTRING (1480362.875856012 596566.8158560097, 1480293.995856002 596496.4358560145), LINESTRING (1480293.995856002 596496.4358560145, 1480235.625856012 596435.0658560097), LINESTRING (1480235.625856012 596435.0658560097, 1480186.495856002 596385.8158560097), LINESTRING (1480186.495856002 596385.8158560097, 1480146.125856012 596345.3758560121), LINESTRING (1480146.125856012 596345.3758560121, 1480141.125856012 596341.3158560097), LINESTRING (1480141.125856012 596341.3158560097, 1480090.995856002 596301.3158560097), LINESTRING (1480090.995856002 596301.3158560097, 1480044.875856012 596263.8758560121), LINESTRING (1480044.875856012 596263.8758560121, 1480005.745856002 596231.995856002), POLYGON ((1480512.6034600246 596678.2322904326, 1480512.598981008 596678.1839810014, 1480511.7765895594 596677.7314431852, 1480512.6034600246 596678.2322904326)), POLYGON ((1480005.745856002 596231.995856002, 1479874.7139810026 596136.0839810073, 1479816.8752310127 596088.7502310127, 1479758.3752310127 596040.8127310127, 1479707.547731012 595999.19689875, 1479758.375856012 596040.8158560097, 1479816.875856012 596088.8158560097, 1479875.375856012 596136.6258560121, 1480005.745856002 596231.995856002)))"

I could skip this error with a rescue block or a next identifier and move on per case but that's not very professional. Documentation on this issue is extremely slim so I figured I'd contribute with a stackexchange question

Plus (just checked), I can't load GEOMETRYCOLLECTION wkts into qgis (currently passing in the polygons & line string arghhs each by hand)

(note: someone with >300 points or whatever mark this question with these new question tags: 'geometry-collection' & 'multi-polygon', they're pretty basic to gis I think, esp the multipolygons, thanks)

EDIT:
I could pluck out the polygon(s) from these collections, but I'm wondering what these linestrings mean (i.e. maybe they're the exterior rings of polygons that are part of the intersect shape)

Best Answer

Solution for now: ignore the LineStrings, these are probably produced by an "st_touches" type spatial join where only the perimeters of the multipolygons are touching

RGeo allows you to break out the pieces of a Geometry Collection so I wrote some code to catch the intersect shape if its a GeometryCollection type, sending it to a different area to pluck out the polygons, sum up their size, mark that as the intersect size, and move on

    x = shape.intersection(p.proj_shape_2264)
    if x.geometry_type.to_s == "GeometryCollection"
      area = 0
      (0..(x.count - 1)).to_a.each do |k| 
        collection_shape = x.geometry_n(k)
        next if collection_shape.geometry_type.to_s == "LineString"
        area  += collection_shape.area            
      end
      z.intersect_size = area
    else

Related Solutions

[GIS] GDAL/OGR ogr.Layer.Intersection fails when producing mixed geometry results

What is strange if the intersection result is a point ?

The intersect predicate is

Returns True if the boundary and interior of the object intersect in any way with those of the other.

With a common point between the geometries, the intersects predicate returns TRUE because the boundary of the first geometry intersects the boundary of the second geometry at this point. The result is always this point

from shapely.geometry import LineString, Polygon
from shapely.wkt import loads
Line1 = loads("LINESTRING ({0} {0}, {1} {0})".format(0, 1))
Line2 = loads("LINESTRING ({1} {0}, {2} {0})".format(0, 1, 2))
Line1.intersects(Line2)
True
Line1.touches(Line2)
True
Line1.intersection(Line2).wkt
'POINT (1 0)'
# polygons
polygon1 = loads('POLYGON (( 140 360, 140 480, 220 480, 220 360, 140 360 ))')
polygon2 = loads('POLYGON (( 220 260, 220 360, 300 360, 300 260, 220 260 ))')
polygon1.intersection(polygon2).wkt
'POINT (220 360)'

In your case the resulting layer is an ogr.wkbPoint

If you change the value of the point, the result is

polygon3 = loads('POLYGON (( 220 260, 220 370, 300 360, 300 260, 220 260 ))')
polygon1.intersection(polygon3).wkt
'LINESTRING (220 370, 220 360)'

polygon4 = loads('POLYGON ((220 260, 200 370, 300 360, 300 260, 220 260))')
polygon1.intersection(polygon4).wkt
'POLYGON ((220 368, 220 360, 201.8181818181818 360, 200 370, 220 368))'

UPDATE 2: With Shapely, there is no problem to clip a roads Shapefile using a polygon even if the road vertex points and the edge of the polygon coincide

roads = loads('LineString (220 360, 120 468)')
polygon1.intersection(roads).wkt
'LINESTRING (220 360, 140 446.4)'

[GIS] Multipolygon created from Scratch Layer

From a formal point of view, I think that the object you have created is a multipolygon because it's a set of distinct polygons which don't satisfy one of the rules that define valid polygons (I report it from a previous answer):

The interior of every Polygon is a connected point set

So, referring to your sample object as a whole, this is sufficient to state that we are not talking about polygons here (instead, a multipolygon may be not a connected space).

Choosing the Polygon or Multipolygon option in the Create New Temporary Scratch Layer dialog only depends on what you are going to do with the layer and QGIS obviously doesn't know that (so, I don't think this could be an issue).

I remark that the difference between the two geometry types is purely formal (since shapefiles don't make a distinction between these objects) and I'll try to explain it better with an example which is partially inspired to @tallistroan's answer.

Regardless of what you choose in the Create New Temporary Scratch Layer dialog, you may edit your layer until to get this shape:

which should be a multipolygon for the reasons explained above. Despite your initial choice doesn't affect what you can subsequently do with your layer, you made a formal mistake if you selected the "Polygon" option in the dialog; if this is the case, you will manage a multipolygon which is formally a polygon because this simple test from the Python Console:

layer = iface.activeLayer()
if layer.isValid():
    print layer.wkbType()

gives 3 as output (3 is the identifier for WKBPolygon types, while 6 is for WKBMultiPolygon, as defined in the QGIS API Documentation).

Made this clear, you were probably looking for a multipart object, which is a valid option for both polygons and multipolygons and takes in account a common attribute value which is shared by one or more features for creating a unique feature.

In the following figure I will try to explain what I intend for multipart object in both cases (in the multipart geometries, same colors belong to the same feature):

Best Answer

Related Solutions

[GIS] GDAL/OGR ogr.Layer.Intersection fails when producing mixed geometry results

[GIS] Multipolygon created from Scratch Layer

Related Question