RGeo The union of two MultyPolygons is returning nil

Question

I want to create a new MultiPolygon from the union of two geometries, but it returns nil.

multipolygon_1 = RGeo::Geos.factory(srid: 4326).parse_wkt("MULTIPOLYGON ...")
multipolygon_2 = RGeo::Geos.factory(srid: 4326).parse_wkt("MULTIPOLYGON ...")

multipolygon_1 + multipolygon_2 # => nil

The MultiPolygons values to reproduce this error can be found at the following gist:

https://gist.github.com/babasbot/926ae326ff3eb4a79601d56288e82a5f

Answer 1

Interesting problem.

You have the chance to have many polygons (34 in total), it makes it easier to debug :

multi_poly_1 = RGeo::Geos.factory(srid: 4326).parse_wkt(wkt_1)
multi_poly_2 = RGeo::Geos.factory(srid: 4326).parse_wkt(wkt_2)

polygons = multi_poly_1.each.to_a + multi_poly_2.each.to_a

There might be a problematic polygon or pair of polygons, so let's try to find them. This code iterates over every combination of indexes and checks if union is nil :

(0...polygons.size).to_a.combination(2).each do |i1, i2|
  if polygons[i1] + polygons[i2] == nil then
    p [i1, i2]
  end
end

It returns

# [1, 11]
# [1, 12]

1 is a non-empty polygon, 11 and 12 look like lines.

Keeping 1 and removing 11 and 12 isn't enough : the union of all the polygons is still nil.

There might still be lines or very flat polygons :

polygons.reject!{|poly| poly.area < 1E-20}
p polygons.size
# 25

Now that 9 polygons are gone, it's possible to calculate the union :

sum = polygons.inject(&:union)
p sum.area - multi_poly_1.area - multi_poly_2.area
#=> -5.800481622797449e-18

The difference in area isn't big, and might come from intersecting polygons or (very) small polygons that have been deleted.