Closed ljwolf closed 10 months ago
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
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.
Merging #261 (7bded4b) into main (9df2081) will increase coverage by
0.3%
. The diff coverage is92.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%) |
:arrow_up: |
@ljwolf sorry, couldn't help myself. But tests are fixed now :D.
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?
Sounds good!
This addresses #260 by:
moment_of_inertia()
function, but not thesecond_areal_moment()
function, which is needed for some of the other compactness measures.I_x + I_y
, notI_{xy}
.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.