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

[ljwolf]

This addresses #260 by:

1. ensuring the center of mass calculations refer to the centroid, rather than taking raw coordinates. This was done correctly in the moment_of_inertia() function, but not the second_areal_moment() function, which is needed for some of the other compactness measures.
2. calculate the sum of moments, rather than the product, for the computation. This is not specified by Li, but is implicit in the representation of the trapezoidal decomposition from Li et al. (2013). Checking that against the general shoelace-theorem decomposition, we find that the correct formulation to match Li is I_x + I_y, not I_{xy}.
3. propagate contributions from multipolygon parts and holes using the parallel axis theorem.
4. add a test bench in the shape module to ensure these cases are computed correctly.

I've compared this to two different implementations from Alan Murray's lab at UCSB (added to the __authors__), submitted over email. Ours reproduces theirs in simple polygons and agrees with Alan's implementation (and gives plausible output) in Multipolygons and we can't compare polygons with holes.

[ljwolf]

For reference, I implemented the "hierarchical" variant, which first computes the MOA for each part of a multipolygon separately (i.e. around the part centroid) and then combines parts onto the multipolygon centroid using the parallel axis theorem. I got the same result in the test bench, so I preferred the simpler implementation that just refers to the elementary parts and the multipolygon centroid.

def second_areal_moment_hierarchical(collection):
    """
    Using equation listed on en.wikipedia.org/wiki/Second_Moment_of_area, the second
    moment of area is the sum of the inertia across the x and y axes

    .. math::
        I_xy = (1/12)\\sum^{i=N}^{i=1} (x_iy_{i+1} + x_i^2 + x_ix_{i+1} + x_{i+1}^2 + y_i^2 + y_iy_{i+1} + y_{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} = y_1. For multipart polygons with holes, 
    all parts are computed first, and then the multipart polygon itself is computed
    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

    """
    ga = _cast(collection)
    original_centroids = shapely.centroid(ga)
    # construct a dataframe of the fundamental parts of all input polygons
    parts, collection_ix = shapely.get_parts(ga, return_index=True)
    part_centroids = shapely.centroid(parts)
    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 *part* center of mass for this implementation
    polygon_rings['part_centroid'] = part_centroids[ring_ix]
    # proportional to the squared distance between the original and part centroids:
    polygon_rings['radius'] = shapely.distance(polygon_rings.centroid, polygon_rings.part_centroid)
    # 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
    part_moas = 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.geometry.area) * ring_in_part.sign).sum()
        ).to_frame("moa")
    # Now, apply the same parallel axis theorem, relating the parts back to the 
    # overall centroid
    part_moas = geopandas.GeoDataFrame(part_moas, geometry=parts)
    part_moas['centroids'] = part_moas.geometry.centroid
    part_moas['total_centroid'] = original_centroids[part_moas.index.get_level_values("collection_ix")]
    part_moas['radius'] = shapely.distance(part_moas.centroid, part_moas.total_centroid)
    return part_moas.groupby(level='collection_ix').apply(lambda part:
        (part.moa + part.radius**2*part.geometry.area).sum()
        ).values

[martinfleis]

That error on file read in envs with numba does not make any sense 🙃. So I just pushed a minor change encoding that polygon @ljwolf wanted to use as WKT to go around it. I couldn't reproduce it locally and have no idea where it comes from.

[codecov]

Codecov Report

Merging #261 (7bded4b) into main (9df2081) will increase coverage by 0.3%.
The diff coverage is 92.3%.

@@           Coverage Diff           @@
##            main    #261     +/-   ##
=======================================
+ Coverage   73.1%   73.4%   +0.3%     
=======================================
  Files         24      24             
  Lines       3272    3284     +12     
  Branches     520     518      -2     
=======================================
+ Hits        2393    2410     +17     
+ Misses       711     708      -3     
+ Partials     168     166      -2     

Files Changed	Coverage Δ
esda/shape.py	72.5% <92.3%> (+4.9%)	⬆️

[martinfleis]

@ljwolf sorry, couldn't help myself. But tests are fixed now :D.

[ljwolf]

No apology necessary @martinfleis, this is super, thank you!

I'll rebase, make @jGaboardi's final edits, merge this, and then cut a release by Friday?

[martinfleis]

Sounds good!