Repository: geopandas/geopandas

Issue:
ENH: Faster .explode() implementation
Hi guys!

Currently the `gdf.explode()` is quite slow as it relies on a [Python loop](https://github.com/geopandas/geopandas/blob/master/geopandas/base.py#L1212). 

@brendan-ward has been working on the [`get_parts()` function](https://github.com/pygeos/pygeos/pull/197) for pygeos which should give quite a nice boost for this functionality.
I was wondering that is the idea to integrate that into `gdf.explode()` function? I'd assume so, but did not find an issue for it here, so doing that now. 🙂 


Discussion:
martinfleis: Nice! Do you have results of that benchmark?
brendan-ward: Before: 1.01±0.05s
After: 33.6±0.6ms
brendan-ward: After refactoring calculation of `inner_index` to use better numpy methods:
13.7±0.2ms.  💥 

Thanks for that suggestion @jorisvandenbossche !
martinfleis: Just a note here. After #1720 we have `explode` implementations in geodataframe.py for `GeoDataFrame` and in geoseries.py for GeoSeries to keep an eye on both.

https://github.com/geopandas/geopandas/blob/7a6fcbbb0bbf02cc4ae0e97d582dc30a1cfa05a0/geopandas/geoseries.py#L634

https://github.com/geopandas/geopandas/blob/7a6fcbbb0bbf02cc4ae0e97d582dc30a1cfa05a0/geopandas/geodataframe.py#L1309
jorisvandenbossche: I don't think it will change a lot (the geodataframe version basically just calls the geoseries version), only the implementation itself was moved from base.py to geoseries.py, so will need to update for that.
brendan-ward: Failing test is unrelated.

Coverage-related tests don't appear to be executing against pygeos (or at least the latest version), so those can probably be ignored.

Patch:
diff --git a/benchmarks/geom_methods.py b/benchmarks/geom_methods.py
index ab6ac612..0da8b4c4 100644
--- a/benchmarks/geom_methods.py
+++ b/benchmarks/geom_methods.py
@@ -2,7 +2,7 @@ import random
 
 import numpy as np
 from geopandas import GeoSeries
-from shapely.geometry import Point, LineString, Polygon
+from shapely.geometry import Point, Polygon, MultiPolygon
 
 
 def with_attributes(**attrs):
@@ -25,10 +25,15 @@ class Bench:
         triangles3 = GeoSeries([Polygon([(random.random(), random.random())
                                          for _ in range(3)])
                                 for _ in range(10000)])
+        triangles4 = GeoSeries([
+            MultiPolygon([
+                Polygon([(random.random(), random.random()) for _ in range(3)])
+            ]) for _ in range(10000)])
         triangle = Polygon([(random.random(), random.random())
                             for _ in range(3)])
         self.triangles, self.triangles2 = triangles, triangles2
         self.triangles_big = triangles3
+        self.multi_triangles = triangles4
         self.triangle = triangle
 
     @with_attributes(param_names=['op'],
@@ -98,7 +103,10 @@ class Bench:
     def time_buffer(self, *args):
         self.points.buffer(2)
 
+    def time_explode(self, *args):
+        self.multi_triangles.explode()
+
 
 # TODO
-# project, interpolate, affine_transform, translate, rotate, scale, skew, explode
+# project, interpolate, affine_transform, translate, rotate, scale, skew
 # cx indexer
diff --git a/ci/envs/37-dev.yaml b/ci/envs/38-dev.yaml
similarity index 86%
rename from ci/envs/37-dev.yaml
rename to ci/envs/38-dev.yaml
index 0f5a270c..5b060fc0 100644
--- a/ci/envs/37-dev.yaml
+++ b/ci/envs/38-dev.yaml
@@ -1,8 +1,8 @@
 name: test
 channels:
-  - defaults
+  - conda-forge
 dependencies:
-  - python=3.7.3
+  - python=3.8
   - cython
   # required
   - fiona
@@ -23,6 +23,7 @@ dependencies:
     - geopy
     - mapclassify>=2.2.0
     # dev versions of packages
+    - git+https://github.com/numpy/numpy.git@master
     - git+https://github.com/pydata/pandas.git@master
     - git+https://github.com/matplotlib/matplotlib.git@master
     - git+https://github.com/Toblerity/Shapely.git@master
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index 413a424e..36501f5c 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -26,14 +26,19 @@ SHAPELY_GE_17 = str(shapely.__version__) >= LooseVersion("1.7.0")
 SHAPELY_GE_18 = str(shapely.__version__) >= LooseVersion("1.8")
 SHAPELY_GE_20 = str(shapely.__version__) >= LooseVersion("2.0")
 
+
 HAS_PYGEOS = None
 USE_PYGEOS = None
 PYGEOS_SHAPELY_COMPAT = None
 
+PYGEOS_GE_09 = None
+
 try:
     import pygeos  # noqa
 
     HAS_PYGEOS = True
+    PYGEOS_GE_09 = str(pygeos.__version__) >= LooseVersion("0.9")
+
 except ImportError:
     HAS_PYGEOS = False
 
diff --git a/geopandas/array.py b/geopandas/array.py
index 114141c7..ea2accc0 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -783,7 +783,9 @@ class GeometryArray(ExtensionArray):
                 "Value should be either a BaseGeometry or None, got %s" % str(value)
             )
         # self.data[idx] = value
-        self.data[idx] = np.array([value], dtype=object)
+        value_arr = np.empty(1, dtype=object)
+        value_arr[:] = [value]
+        self.data[idx] = value_arr
         return self
 
     def fillna(self, value=None, method=None, limit=None):
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 59b6d925..e9381a96 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -15,7 +15,7 @@ from geopandas.plotting import plot_series
 from .array import GeometryArray, GeometryDtype, from_shapely
 from .base import is_geometry_type
 from . import _vectorized as vectorized
-from ._compat import ignore_shapely2_warnings
+from . import _compat as compat
 
 
 _SERIES_WARNING_MSG = """\
@@ -197,7 +197,7 @@ class GeoSeries(GeoPandasBase, Series):
             # https://github.com/pandas-dev/pandas/issues/26469
             kwargs.pop("dtype", None)
             # Use Series constructor to handle input data
-            with ignore_shapely2_warnings():
+            with compat.ignore_shapely2_warnings():
                 s = pd.Series(data, index=index, name=name, **kwargs)
             # prevent trying to convert non-geometry objects
             if s.dtype != object:
@@ -669,6 +669,39 @@ class GeoSeries(GeoPandasBase, Series):
         GeoDataFrame.explode
 
         """
+
+        if compat.USE_PYGEOS and compat.PYGEOS_GE_09:
+            import pygeos  # noqa
+
+            geometries, outer_idx = pygeos.get_parts(
+                self.values.data, return_index=True
+            )
+
+            if len(outer_idx):
+                # Generate inner index as a range per value of outer_idx
+                # 1. identify the start of each run of values in outer_idx
+                # 2. count number of values per run
+                # 3. use cumulative sums to create an incremental range
+                #    starting at 0 in each run
+                run_start = np.r_[True, outer_idx[:-1] != outer_idx[1:]]
+                counts = np.diff(np.r_[np.nonzero(run_start)[0], len(outer_idx)])
+                inner_index = (~run_start).cumsum()
+                inner_index -= np.repeat(inner_index[run_start], counts)
+
+            else:
+                inner_index = []
+
+            # extract original index values based on integer index
+            outer_index = self.index.take(outer_idx)
+
+            index = MultiIndex.from_arrays(
+                [outer_index, inner_index], names=self.index.names + [None]
+            )
+
+            return GeoSeries(geometries, index=index, crs=self.crs).__finalize__(self)
+
+        # else PyGEOS is not available or version <= 0.8
+
         index = []
         geometries = []
         for idx, s in self.geometry.iteritems():


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index b93554bc..a1327c8e 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -51,7 +51,7 @@ jobs:
             os: windows-latest
             postgis: false
             dev: false
-          - env: ci/envs/37-dev.yaml
+          - env: ci/envs/38-dev.yaml
             os: ubuntu-latest
             dev: true
 
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index 96b11a17..59511223 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -159,7 +159,6 @@ def test_to_file_types(tmpdir, df_points):
     """ Test various integer type columns (GH#93) """
     tempfilename = os.path.join(str(tmpdir), "int.shp")
     int_types = [
-        np.int,
         np.int8,
         np.int16,
         np.int32,
@@ -169,7 +168,6 @@ def test_to_file_types(tmpdir, df_points):
         np.uint16,
         np.uint32,
         np.uint64,
-        np.long,
     ]
     geometry = df_points.geometry
     data = dict(
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 3f128f93..9b0f0e0d 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -52,8 +52,8 @@ def equal_geometries(result, expected):
 
 
 def test_points():
-    x = np.arange(10).astype(np.float)
-    y = np.arange(10).astype(np.float) ** 2
+    x = np.arange(10).astype(np.float64)
+    y = np.arange(10).astype(np.float64) ** 2
 
     points = points_from_xy(x, y)
     assert isinstance(points, GeometryArray)
@@ -500,7 +500,7 @@ def test_unary_float(attr):
     na_value = np.nan
     result = getattr(T, attr)
     assert isinstance(result, np.ndarray)
-    assert result.dtype == np.float
+    assert result.dtype == np.dtype("float64")
     expected = [getattr(t, attr) if t is not None else na_value for t in triangles]
     np.testing.assert_allclose(result, expected)
 
diff --git a/geopandas/tests/test_extension_array.py b/geopandas/tests/test_extension_array.py
index 5c2a01d3..37d37a3d 100644
--- a/geopandas/tests/test_extension_array.py
+++ b/geopandas/tests/test_extension_array.py
@@ -16,6 +16,7 @@ expected to be available to pytest by the inherited pandas tests).
 import operator
 
 import numpy as np
+from numpy.testing import assert_array_equal
 import pandas as pd
 from pandas.tests.extension import base as extension_tests
 
@@ -46,7 +47,8 @@ def dtype():
 
 
 def make_data():
-    a = np.array([shapely.geometry.Point(i, i) for i in range(100)], dtype=object)
+    a = np.empty(100, dtype=object)
+    a[:] = [shapely.geometry.Point(i, i) for i in range(100)]
     ga = from_shapely(a)
     return ga
 
@@ -287,6 +289,19 @@ class TestDtype(extension_tests.BaseDtypeTests):
 
 
 class TestInterface(extension_tests.BaseInterfaceTests):
+    def test_array_interface(self, data):
+        # we are overriding this base test because the creation of `expected`
+        # potentionally doesn't work for shapely geometries
+        # TODO can be removed with Shapely 2.0
+        result = np.array(data)
+        assert result[0] == data[0]
+
+        result = np.array(data, dtype=object)
+        # expected = np.array(list(data), dtype=object)
+        expected = np.empty(len(data), dtype=object)
+        expected[:] = list(data)
+        assert_array_equal(result, expected)
+
     def test_contains(self, data, data_missing):
         # overrided due to the inconsistency between
         # GeometryDtype.na_value = np.nan