From b3e1e2fcbf2bd04509ce1fe399de7f48a0cc4842 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Wed, 23 Aug 2023 16:39:02 +0200 Subject: [PATCH] compat fix --- esda/shape.py | 99 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 59 insertions(+), 40 deletions(-) diff --git a/esda/shape.py b/esda/shape.py index ce8ac7a0..5a165048 100644 --- a/esda/shape.py +++ b/esda/shape.py @@ -14,9 +14,10 @@ "Martin Fleischmann ", "Levi John Wolf ", "Alan Murray ", - "Jiwan Baik " + "Jiwan Baik ", ) + # -------------------- UTILITIES --------------------# def _cast(collection): """ @@ -477,44 +478,48 @@ def second_areal_moment(collection): 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 treated as separate contributions to the overall centroid, which + and where 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. + 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. + 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) - import geopandas # function level, to follow module design - + 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 + # 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 + # 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 + # 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) + 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: - + # 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() @@ -522,31 +527,40 @@ def second_areal_moment(collection): # 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 + 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) + 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] + 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, polygon_rings.collection_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 - 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 + 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. + for an input polygon. """ coordinates = shapely.get_coordinates(polygon) centroid = shapely.centroid(polygon) @@ -568,10 +582,15 @@ def _second_moa_ring_xplusy(points): xhyt = x_head * y_tail xtyt = x_tail * y_tail xhyh = x_head * y_head - 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 + 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 -------------------- #