Repository: pysal/momepy

Issue:
BUG: FormFactor uses incorrect formula
I just realised that `FormFactor` is using an incorrect formula. 

The formula in the original paper by Bourdic et al says `S/(V^2/3)` where `S` is surface. However, I have implemented it as `A/(V^2/3)`, where `A` is an area of a polygon, not the theoretical surface of the whole building.

I think we should fix this with `S` defined as `perimeter * height + area`.

I'll try to get it done today and cut the 0.6.0 release to have it ready for the PySAL meta release later this month.

Discussion:

Patch:
diff --git a/momepy/shape.py b/momepy/shape.py
index 717465d..9da9928 100644
--- a/momepy/shape.py
+++ b/momepy/shape.py
@@ -32,12 +32,34 @@ __all__ = [
 ]
 
 
+def _form_factor(height, geometry, area=None, perimeter=None, volume=None):
+    """helper for FormFactor"""
+    if area is None:
+        area = geometry.area
+    if perimeter is None:
+        perimeter = geometry.length
+    if volume is None:
+        volume = area * height
+
+    surface = (perimeter * height) + area
+    zeros = volume == 0
+    res = np.empty(len(geometry))
+    res[zeros] = np.nan
+    res[~zeros] = surface[~zeros] / (volume[~zeros] ** (2 / 3))
+    return res
+
+
 class FormFactor:
     """
     Calculates form factor of each object in given GeoDataFrame.
 
     .. math::
-        area \\over {volume^{2 \\over 3}}
+        surface \\over {volume^{2 \\over 3}}
+
+    where
+
+    .. math::
+        surface = (perimeter * height) + area
 
     Adapted from :cite:`bourdic2012`.
 
@@ -78,7 +100,11 @@ class FormFactor:
 
     """
 
-    def __init__(self, gdf, volumes, areas=None):
+    def __init__(self, gdf, volumes, areas=None, heights=None):
+        if heights is None:
+            raise ValueError("`heights` cannot be None.")
+            # TODO: this shouldn't be needed but it would be a breaking change now.
+            # remove during the functional refactor
         self.gdf = gdf
 
         gdf = gdf.copy()
@@ -92,11 +118,19 @@ class FormFactor:
             gdf["mm_a"] = areas
             areas = "mm_a"
         self.areas = gdf[areas]
-        zeros = gdf[volumes] == 0
-        res = np.empty(len(gdf))
-        res[zeros] = 0
-        res[~zeros] = gdf[areas][~zeros] / (gdf[volumes][~zeros] ** (2 / 3))
-        self.series = pd.Series(res, index=gdf.index)
+
+        if isinstance(heights, str):
+            heights = gdf[heights]
+
+        self.series = pd.Series(
+            _form_factor(
+                height=np.asarray(heights),
+                geometry=gdf.geometry,
+                area=self.areas,
+                volume=self.volumes,
+            ),
+            index=gdf.index,
+        )
 
 
 class FractalDimension:


Test patch:
diff --git a/tests/test_shape.py b/tests/test_shape.py
index e431dd1..f116d93 100644
--- a/tests/test_shape.py
+++ b/tests/test_shape.py
@@ -20,21 +20,19 @@ class TestShape:
         self.df_buildings["volume"] = mm.Volume(self.df_buildings, "height").series
 
     def test_FormFactor(self):
-        self.df_buildings["ff"] = mm.FormFactor(self.df_buildings, "volume").series
-        check = (self.df_buildings.geometry[0].area) / (
-            self.df_buildings.volume[0] ** (2 / 3)
-        )
-        assert self.df_buildings["ff"][0] == check
+        self.df_buildings["ff"] = mm.FormFactor(
+            self.df_buildings, "volume", heights="height"
+        ).series
+        check = 5.4486362624193
+        assert self.df_buildings["ff"].mean() == approx(check)
 
         self.df_buildings["ff"] = mm.FormFactor(
             self.df_buildings,
             mm.Volume(self.df_buildings, "height").series,
             areas=self.df_buildings.geometry.area,
+            heights=self.df_buildings["height"],
         ).series
-        check = (self.df_buildings.geometry[0].area) / (
-            self.df_buildings.volume[0] ** (2 / 3)
-        )
-        assert self.df_buildings["ff"][0] == check
+        assert self.df_buildings["ff"].mean() == approx(check)
 
     def test_FractalDimension(self):
         self.df_buildings["fd"] = mm.FractalDimension(self.df_buildings).series