Fix to calculation of normalized mass moment of inertia (in moment_of_inertia_regions function)

[leehach]

Fix for #437 .

Normalized mass moment of inertia is calculated with:

moment_of_inertia_regions(geometry_array, normalize=True, regions=<geom region assignment>, weights=<geom weight>)

Error was discovered in that the moment of inertia calculated for each geometry in the region was not adjusted for its mass (weight) and area. Depending upon the CRS, the area adjustment could be large. The mass moment of inertia for the reference circle was calculated correctly, because the formula is very straightforward ($mA/2\pi$).

The NMMI is the MMI of the reference circle divided by the (nonnormalized) MMI of the region. The following is a general strategy for confirming the expected NMMI with any input geometry.

Assume a geodataframe gdf with region assignments in the column "region" and weights in the column "mass".

import numpy as np
import geopandas as gpd
import shapely
import esda

# Calculate nonnormalized MOI of region
region_moi = esda.shape.moment_of_inertia_regions(gdf, regions="region", weights="mass").values

# Calculate areas and masses of regions for MOI of reference disc
masses = gdf.groupby("region")["mass"].sum().values
gdf["area"] = gdf.geometry.area
areas = gdf.groupby("region")["area"].sum().values
ref_disc = shapely.Point(0,0).buffer(np.sqrt(areas / np.pi))

# Calculate MOI of reference disc
ref_disc_moi = esda.shape.moment_of_inertia(ref_disc) * masses / areas

# Calculate NMMI manually and using function
nmmi_manual = ref_disc_moi / region_moi
nmmi_function = esda.shape.moment_of_inertia_regions(gdf, normalize=True, regions="region", weights="mass").values

# Calculate proportion difference between manual calculation and function
(nmmi_manual - nmmi_function) / nmmi_function

When executed on the virginia counties built-in dataset with counties grouped by Virginia economic regions, the manually calculated NMMI was 0.3% smaller than the NMMI calculated by the function.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 82.2%. Comparing base (74702c9) to head (ebd7077).
⚠️ Report is 9 commits behind head on main.

Additional details and impacted files

@@          Coverage Diff          @@
##            main    #438   +/-   ##
=====================================
  Coverage   82.2%   82.2%           
=====================================
  Files         27      27           
  Lines       3969    3969           
=====================================
  Hits        3264    3264           
  Misses       705     705           

Files with missing lines	Coverage Δ
esda/shape.py	73.6% <ø> (ø)

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[martinfleis]

One thing to note - the PR titles are automatically pulled to the release notes to form a changelog. I'd suggest to make them a bit more suited for the purpose :).

[leehach]

One thing to note - the PR titles are automatically pulled to the release notes to form a changelog. I'd suggest to make them a bit more suited for the purpose :).

Acknowledged and changed.