Fix unhandled GeometryCollection in BaseGeometry.intersection

Solves issue #566.

Author: molinav

Diff for: packages/basemap/src/_geoslib.pyx

@@ -318,6 +318,30 @@ cdef class BaseGeometry:
                 b = _get_coords(gout)
                 p = LineString(b)
                 pout.append(p)
+        elif typeid == GEOS_GEOMETRYCOLLECTION:
+            numgeoms = GEOSGetNumGeometries(g3)
+            pout = []
+            for i from 0 <= i < numgeoms:
+                gout = GEOSGetGeometryN(g3, i)
+                typeid = GEOSGeomTypeId(gout)
+                if typeid == GEOS_POLYGON:
+                    b = _get_coords(gout)
+                    p = Polygon(b)
+                    pout.append(p)
+                elif typeid == GEOS_LINESTRING:
+                    b = _get_coords(gout)
+                    p = LineString(b)
+                    pout.append(p)
+                else:
+                    # More cases might need to be handled here:
+                    # - GEOS_MULTILINESTRING
+                    # - GEOS_MULTIPOLYGON
+                    # - GEOS_GEOMETRYCOLLECTION
+                    # The underlying problem is the need of a generic
+                    # converter from GEOSGeom pointers to `_geoslib`
+                    # objects, since postprocessing `GeometryCollections`
+                    # might need recursiveness.
+                    pass
         else:
             #type = PyBytes_FromString(GEOSGeomType(g3))
             #raise NotImplementedError("intersections of type '%s' not yet implemented" % (type))

Diff for: packages/basemap/test/mpl_toolkits/basemap/test_Basemap.py

@@ -8,6 +8,7 @@
 import numpy as np
 import matplotlib as mpl
 import matplotlib.pyplot as plt
+from matplotlib.collections import LineCollection
 from matplotlib.image import AxesImage
 from matplotlib.patches import Polygon
 from mpl_toolkits.basemap import Basemap
@@ -85,6 +86,46 @@ def test_init_with_optional_casting(self):
         bmap2 = Basemap(lat_1=bmap2_lat_1, **kwds)
         self.assertEqual(bmap1.proj4string, bmap2.proj4string)
 
+    def test_drawcoastlines(self, axs=None, axslen0=10):
+        """Test that no lines are missing when drawing coastlines."""
+
+        axs_obj = plt.gca() if axs is None else axs
+        axs_children = axs_obj.get_children()
+        self.assertEqual(len(axs_children), axslen0)
+
+        bmap = Basemap(projection="merc", resolution="i", lat_ts=20,
+                       llcrnrlat=36.0, llcrnrlon=6.0,
+                       urcrnrlat=47.7, urcrnrlon=19.0)
+
+        collection = bmap.drawcoastlines(linewidth=1, color="red")
+        self.assertIsInstance(collection, LineCollection)
+
+        axs_children = axs_obj.get_children()
+        self.assertEqual(len(axs_children), axslen0 + 1)
+
+        lines = collection.get_paths()
+        self.assertEqual(len(lines), 27)
+
+    def test_drawcountries(self, axs=None, axslen0=10):
+        """Test that no lines are missing when drawing country boundaries."""
+
+        axs_obj = plt.gca() if axs is None else axs
+        axs_children = axs_obj.get_children()
+        self.assertEqual(len(axs_children), axslen0)
+
+        bmap = Basemap(projection="merc", resolution="i", lat_ts=20,
+                       llcrnrlat=36.0, llcrnrlon=6.0,
+                       urcrnrlat=47.7, urcrnrlon=19.0)
+
+        collection = bmap.drawcountries(linewidth=1, color="blue")
+        self.assertIsInstance(collection, LineCollection)
+
+        axs_children = axs_obj.get_children()
+        self.assertEqual(len(axs_children), axslen0 + 1)
+
+        lines = collection.get_paths()
+        self.assertEqual(len(lines), 29)
+
     @unittest.skipIf(PIL is None, reason="pillow unavailable")
     def test_arcgisimage_with_cyl(self, axs=None, axslen0=10):
         """Test showing an ArcGIS image as background."""