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
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.