fix second areal moment calculation, cascade changes down to other stats

.gitignore

@@ -57,7 +57,6 @@ notebooks/splot.ipynb
 notebooks/test.ipynb
 pysal_data/
 spdep-geary.ipynb
-tests/
 tools/changelog.md
 tools/changelog_2.0.0.md
 tools/changelog_2.0.1.md

esda/shape.py

@@ -10,6 +10,14 @@
 from .crand import njit, prange
 
 
+__author__ = (
+    "Martin Fleischmann <martin@fleischmann.net>",
+    "Levi John Wolf <levi.john.wolf@gmail.com>",
+    "Alan Murray <amurray@ucsb.edu>",
+    "Jiwan Baik <jiwon.baik@geog.ucsb.edu>",
+)
+
+
 # -------------------- UTILITIES --------------------#
 def _cast(collection):
     """
@@ -463,43 +471,111 @@
 
 def second_areal_moment(collection):
     """
-    Using equation listed on en.wikipedia.org/Second_Moment_of_area, the second
-    moment of area is actually the cross-moment of area between the X and Y
-    dimensions:
+    Using equation listed on en.wikipedia.org/wiki/Second_moment_of_area#Any_polygon, the second
+    moment of area is the sum of the inertia across the x and y axes:
 
+    The :math:`x` axis is given by:
     .. math::
-        I_xy = (1/24)\\sum^{i=N}^{i=1} (x_iy_{i+1} + 2*x_iy_i + 2*x_{i+1}y_{i+1} +
-        x_{i+1}y_i)(x_iy_i - x_{i+1}y_i)
+        I_x = (1/12)\\sum^{N}_{i=1} (x_i y_{i+1} - x_{i+1}y_i) (x_i^2 + x_ix_{i+1} + x_{i+1}^2)
 
-    where x_i, y_i is the current point and x_{i+1}, y_{i+1} is the next point,
-    and where x_{n+1} = x_1, y_{n+1} = 1.
+    While the :math:`y` axis is in a similar form:
+    .. math::
+        I_y = (1/12)\\sum^{N}_{i=1} (x_i y_{i+1} - x_{i+1}y_i) (y_i^2 + y_iy_{i+1} + y_{i+1}^2)
 
-    This relation is known as the:
-    - second moment of area
-    - moment of inertia of plane area
-    - area moment of inertia
-    - second area moment
+    where :math:`x_i`,:math:`y_i` is the current point and :math:`x_{i+1}`, :math:`y_{i+1}` is the next point,
+    and where :math:`x_{n+1} = x_1, y_{n+1} = y_1`. For multipart polygons with holes, 
+    all parts are treated as separate contributions to the overall centroid, which 
+    provides the same result as if all parts with holes are separately computed, and then
+    merged together using the parallel axis theorem.
+
+    References
+    ----------
+    Hally, D. 1987. The calculations of the moments of polygons. Canadian National
+    Defense Research and Development Technical Memorandum 87/209.
+    https://apps.dtic.mil/dtic/tr/fulltext/u2/a183444.pdf
 
-    and is *not* the mass moment of inertia, a property of the distribution of
-    mass around a shape.
     """
     ga = _cast(collection)
-    result = numpy.zeros(len(ga))
-    n_holes_per_geom = shapely.get_num_interior_rings(ga)
-    for i, geometry in enumerate(ga):
-        n_holes = n_holes_per_geom[i]
-        for hole_ix in range(n_holes):
-            hole = shapely.get_coordinates(shapely.get_interior_ring(ga, hole_ix))
-            result[i] -= _second_moa_ring(hole)
-        for part in shapely.get_parts(geometry):
-            result[i] += _second_moa_ring(shapely.get_coordinates(part))
-    # must divide everything by 24 and flip if polygon is clockwise.
-    signflip = numpy.array([-1, 1])[shapely.is_ccw(ga).astype(int)]
-    return result * (1 / 24) * signflip
+    import geopandas  # function level, to follow module design
+
+    # construct a dataframe of the fundamental parts of all input polygons
+    parts, collection_ix = shapely.get_parts(ga, return_index=True)
+    rings, ring_ix = shapely.get_rings(parts, return_index=True)
+    # get_rings() always returns the exterior first, then the interiors
+    collection_ix = numpy.repeat(
+        collection_ix, shapely.get_num_interior_rings(parts) + 1
+    )
+    # we need to work in polygon-space for the algorithms (centroid, shoelace calculation) to work
+    polygon_rings = shapely.polygons(rings)
+    is_external = numpy.zeros_like(collection_ix).astype(bool)
+    # the first element is always external
+    is_external[0] = True
+    # and each subsequent element is external iff it is different from the preceeding index
+    is_external[1:] = ring_ix[1:] != ring_ix[:-1]
+    # now, our analysis frame contains a bunch of (guaranteed-to-be-simple) polygons
+    # that represent either exterior rings or holes
+    polygon_rings = geopandas.GeoDataFrame(
+        dict(
+            collection_ix=collection_ix,
+            ring_within_geom_ix=ring_ix,
+            is_external_ring=is_external,
+        ),
+        geometry=polygon_rings,
+    )
+    # the polygonal moi can be calculated using the same ring-based strategy,
+    # and this could be parallelized if necessary over the elemental shapes with:
+
+    # from joblib import Parallel, parallel_backend, delayed
+    # with parallel_backend('loky', n_jobs=-1):
+    #     engine = Parallel()
+    #     promise = delayed(_second_moment_of_area_polygon)
+    #     result = engine(promise(geom) for geom in polygon_rings.geometry.values)
+
+    # but we will keep simple for now
+    polygon_rings["moa"] = polygon_rings.geometry.apply(_second_moment_of_area_polygon)
+    # the above algorithm computes an unsigned moa to be insensitive to winding direction.
+    # however, we need to subtract the moa of holes. Hence, the sign of the moa is
+    # -1 when the polygon is an internal ring and 1 otherwise:
+    polygon_rings["sign"] = (1 - polygon_rings.is_external_ring * 2) * -1
+    # shapely already uses the correct formulation for centroids
+    polygon_rings["centroids"] = shapely.centroid(polygon_rings.geometry)
+    # the inertia of parts applies to the overall center of mass:
+    original_centroids = shapely.centroid(ga)
+    polygon_rings["collection_centroid"] = original_centroids[collection_ix]
+    # proportional to the squared distance between the original and part centroids:
+    polygon_rings["radius"] = shapely.distance(
+        polygon_rings.centroid.values, polygon_rings.collection_centroid.values
+    )
+    # now, we take the sum of (I+Ar^2) for each ring, treating the
+    # contribution of holes as negative. Then, we take the sum of all of the contributions
+    return (
+        polygon_rings.groupby(["collection_ix", "ring_within_geom_ix"])
+        .apply(
+            lambda ring_in_part: (
+                (ring_in_part.moa + ring_in_part.radius**2 * ring_in_part.area)
+                * ring_in_part.sign
+            ).sum()
+        )
+        .groupby(level="collection_ix")
+        .sum()
+        .values
+    )
+
+
+def _second_moment_of_area_polygon(polygon):
+    """
+    Compute the absolute value of the moment of area (i.e. ignoring winding direction)
+    for an input polygon.
+    """
+    coordinates = shapely.get_coordinates(polygon)
+    centroid = shapely.centroid(polygon)
+    centroid_coords = shapely.get_coordinates(centroid)
+    moi = _second_moa_ring_xplusy(coordinates - centroid_coords)
+    return abs(moi)
 
 
 @njit
-def _second_moa_ring(points):
+def _second_moa_ring_xplusy(points):
     """
     implementation of the moment of area for a single ring
     """
@@ -511,8 +587,15 @@
         xhyt = x_head * y_tail
         xtyt = x_tail * y_tail
         xhyh = x_head * y_head
-        moi += (xtyh - xhyt) * (xtyh + 2 * xtyt + 2 * xhyh + xhyt)
-    return moi
+        moi += (xtyh - xhyt) * (
+            x_head**2
+            + x_head * x_tail
+            + x_tail**2
+            + y_head**2
+            + y_head * y_tail
+            + y_tail**2
+        )
+    return moi / 12
 
 
 # -------------------- OTHER MEASURES -------------------- #

esda/tests/test_shape.py

@@ -1,15 +1,15 @@
 import pytest
 from numpy import array, testing
-from shapely import geometry
-
+import shapely
 import esda
+import geopandas, numpy
 
 pytest.importorskip("numba")
 
 
 shape = array(
     [
-        geometry.Polygon(
+        shapely.geometry.Polygon(
             [
                 (0, 0),
                 (0.25, 0.25),
@@ -26,6 +26,49 @@
 
 ATOL = 0.001
 
+## for a hole/multipart testbench:
+
+test_geom_translated = shapely.from_wkt(
+    [
+        "POLYGON ((-3.1823503126754247 0.085191513232644, -3.2545854200972997 0.271135116748269, "
+        "-3.2472001661910497 0.296769882373269, -3.2779008497847997 0.3333146821779565, "
+        "-3.286461030448862 0.4668824922365502, -3.312919770683237 0.4887788545412377, "
+        "-3.308738862480112 0.5528352510256127, -3.271751557792612 0.6005037080568627, "
+        "-3.2749711622847997 0.6791780244631127, -3.301475678886362 0.7266938936037377, "
+        "-3.3279496779097997 0.7266252290529565, -3.3439103712691747 0.8217256318849877, "
+        "-3.385307465995737 0.9057634126467065, -3.385490571464487 1.0868776094240502, "
+        "-3.4014665236129247 1.1563432466310815, -3.3219682325972997 1.1584108125490502, "
+        "-3.3168107618941747 1.328538681201394, -3.2181779493941747 1.3360688936037377, "
+        "-3.2150346388472997 1.4000871431154565, -3.1089555372847997 1.4057023774904565, "
+        "-3.105888520683237 1.9232576143068627, -2.702140962089487 1.9063203584474877, "
+        "-2.7046891798629247 0.8271501313967065, -2.749656831230112 0.7665956270021752, "
+        "-2.671684419120737 0.5808656465334252, -2.8658219923629247 0.3313920747560815, "
+        "-2.9108354200972997 0.1093156587404565, -3.1823503126754247 0.085191513232644))"
+    ]
+)[0]
+
+test_simple = shapely.difference(
+    shapely.box(-1, 0, -2, 1), shapely.box(-1, 0, -1.5, 0.5)
+)
+
+test_hole = shapely.difference(
+    shapely.box(0, 0, 1.81, 1.81), shapely.box(0.8, 0.8, 1.6, 1.6)
+)
+
+test_mp = shapely.union(shapely.box(-1, -1, -1.5, -2), shapely.box(0, -1, 1.25, -2))
+
+test_mp_hole = shapely.union(
+    shapely.transform(
+        test_hole, lambda x: numpy.column_stack((-x[:, 0] + 3, x[:, 1] * 0.5 + 3))
+    ),
+    shapely.transform(test_hole, lambda x: numpy.column_stack((x[:, 0] + 4, x[:, 1]))),
+)
+
+testbench = geopandas.GeoDataFrame(
+    geometry=[test_geom_translated, test_simple, test_mp, test_hole, test_mp_hole]
+).reset_index()
+testbench["name"] = ["Hanock County", "Simple", "Multi", "Single Hole", "Multi Hole"]
+
 
 def test_boundary_amplitude():
     observed = esda.shape.boundary_amplitude(shape)
@@ -62,19 +105,28 @@ def test_ipq():
     testing.assert_allclose(observed, 0.387275, atol=ATOL)
 
 
-def test_moa():
-    observed = esda.shape.moa_ratio(shape)
-    testing.assert_allclose(observed, 3.249799, atol=ATOL)
-
-
 def test_moment_of_interia():
     observed = esda.shape.moment_of_inertia(shape)
     testing.assert_allclose(observed, 0.315715, atol=ATOL)
 
 
+def test_second_areal_moment():
+    observed = esda.shape.second_areal_moment(testbench.geometry)
+    testing.assert_allclose(
+        observed,
+        [0.23480628, 0.11458333, 1.57459077, 1.58210246, 14.18946959],
+        atol=ATOL,
+    )
+
+
+def test_moa():
+    observed = esda.shape.moa_ratio(shape)
+    testing.assert_allclose(observed, 5.35261, atol=ATOL)
+
+
 def test_nmi():
     observed = esda.shape.nmi(shape)
-    testing.assert_allclose(observed, 0.487412, atol=ATOL)
+    testing.assert_allclose(observed, 0.802796, atol=ATOL)
 
 
 def test_mbc():