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

Question

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)

Answer 1

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

Answer 2

I recently had two geometries whose intersection produced a GeometryCollection of Polygon and tiny LineString. My "fix" was to erode the source geometry slightly before intersection:

shape.buffer(1e-5).intersection(p.proj_shape_2264)

Now the intersection yields just one Polygon.

Update: I ran into this again in a situation where the application requires the geometry not be eroded. I wrote something similar to boulder_ruby's answer, but I only pas along Polygon and MultiPolygons. Here's my implementation:

  # Returns geometry of type Polygon or MultiPolygon and never a
  # GeometryCollection for the given geometry.
  def as_polygonal(geom)
    # non-intersecting geometries will return an empty
    # GeometryCollection... we only want to "fix" things which have
    # linestrings inside the intersection.
    if geom.geometry_type == RGeo::Feature::GeometryCollection and geom.count > 0
      acceptable_types = [RGeo::Feature::Polygon, RGeo::Feature::MultiPolygon]
      polygonal = geom.select {|g| acceptable_types.include?(g.geometry_type) }

      geom = RGeo::Geos.factory(srid: 4326).parse_wkt('POLYGON EMPTY')
      polygonal.each {|p| geom = geom.union(p)}
    end

    geom
  end