RGeo Combining two MultyPolygons returns zero

Question

I want to create a new MultiPolygon from a 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 values ​​of MultiPolygons for reproducing this error can be found in the following form:

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

Answer 1

An interesting problem.

You have the opportunity to have many polygons (34 in total), this makes debugging easier:

 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 

It may be a problematic polygon or a pair of polygons, so try to find them. This code repeats through every combination of indices and checks if the 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 

He returns

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

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

Saving 1 and deleting 11 and 12 is not enough: the union of all polygons is still zero.

There may be lines or very flat polygons:

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

Now that 9 polygons have disappeared, you can calculate the union:

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

The difference in the area is small and can arise from intersecting polygons or (very) small polygons that have been deleted.