[WIP] shared perimeter-weighted contiguity

libpysal/weights/contiguity.py

@@ -2,16 +2,18 @@
 import warnings
 
 import numpy
+import pandas as pd
 
 from ..cg import voronoi_frames
 from ..io.fileio import FileIO
 from ._contW_lists import ContiguityWeightsLists
+from .raster import da2W, da2WSP
 from .util import get_ids, get_points_array
 from .weights import WSP, W
-from .raster import da2W, da2WSP
 
 try:
     from shapely.geometry import Point as shapely_point
+
     from ..cg.shapes import Point as pysal_point
 
     point_type = (shapely_point, pysal_point)
@@ -141,6 +143,8 @@ def from_dataframe(
         ids=None,
         id_order=None,
         use_index=None,
+        perimeter=False,
+        perim_std=False,
         **kwargs,
     ):
         """
@@ -175,6 +179,12 @@ def from_dataframe(
         use_index   : bool
                       use index of `df` as `ids` to index the spatial weights object.
                       Defaults to False but in future will default to True.
+        perimeter   : bool
+                      if True, use the length of the shared boundary between adjacent units as
+                      the weight value
+        perim_std   : bool
+                      if True, (and `perimeter==True`), then perimeter weights are set to the
+                      ratio of the shared boudary length to the focal unit's perimeter
 
         See Also
         --------
@@ -239,9 +249,13 @@ def from_dataframe(
         if id_order is None:
             id_order = ids
 
-        return cls.from_iterable(
+        w = cls.from_iterable(
             df[geom_col].tolist(), ids=ids, id_order=id_order, **kwargs
         )
+        if perimeter:
+            w = _return_length_weighted_w(w, df, perim_std)
+        return w
+
 
     @classmethod
     def from_xarray(
@@ -427,6 +441,8 @@ def from_dataframe(
         ids=None,
         id_order=None,
         use_index=None,
+        perimeter=False,
+        perim_std=False,
         **kwargs,
     ):
         """
@@ -461,6 +477,12 @@ def from_dataframe(
         use_index   : bool
                       use index of `df` as `ids` to index the spatial weights object.
                       Defaults to False but in future will default to True.
+        perimeter   : bool
+                      if True, use the length of the shared boundary between adjacent units as
+                      the weight value
+        perim_std   : bool
+                      if True, (and `perimeter==True`), then perimeter weights are set to the
+                      ratio of the shared boudary length to the focal unit's perimeter
 
         See Also
         --------
@@ -525,9 +547,13 @@ def from_dataframe(
         if id_order is None:
             id_order = ids
 
-        return cls.from_iterable(
+        w =  cls.from_iterable(
             df[geom_col].tolist(), ids=ids, id_order=id_order, **kwargs
         )
+        if perimeter:
+            w = _return_length_weighted_w(w, df, perim_std)
+        return w
+
 
     @classmethod
     def from_xarray(
@@ -748,3 +774,68 @@ def buildContiguity(polygons, criterion="rook", ids=None):
         return Queen(polygons, ids=ids)
     else:
         raise Exception('Weights criterion "{}" was not found.'.format(criterion))
+
+
+def _return_length_weighted_w(w, data, perimeter_standardize):
+    """Return a W object whose value is the length of the common boundary of two areal units that share border.
+
+    Parameters
+    ----------
+    w : libpsal.weights.Rook
+    data : pandas.DataFrame
+    perimeter_standardize: bool
+        if True, scale the weight value equal to the shared
+        boundary divided by the total boundary of the focal unit.
+
+    Returns
+    --------
+    w : libpysal.weights.W
+        weights object with values equal to the shared border of ij
+
+    """
+    try:
+        import geopandas as gpd
+    except ImportError as e:
+        raise e('You must have geopandas installed to create perimeter-weighted weights')
+    adjlist = w.to_adjlist()
+    islands = pd.DataFrame.from_records(
+        [{"focal": island, "neighbor": island, "weight": 0} for island in w.islands]
+    )
+    merged = adjlist.merge(
+        data.geometry.to_frame("geometry"),
+        left_on="focal",
+        right_index=True,
+        how="left",
+    ).merge(
+        data.geometry.to_frame("geometry"),
+        left_on="neighbor",
+        right_index=True,
+        how="left",
+        suffixes=("_focal", "_neighbor"),
+    )
+    # Transforming from pandas to geopandas
+    merged = gpd.GeoDataFrame(merged, geometry="geometry_focal")
+    total_boundary_length = merged.boundary.length
+    merged["geometry_neighbor"] = gpd.GeoSeries(merged.geometry_neighbor)
+
+    # Getting the shared boundaries
+    merged["shared_boundary"] = merged.geometry_focal.intersection(
+        merged.set_geometry("geometry_neighbor")
+    )
+
+    # Putting it back to a matrix
+    if perimeter_standardize:
+        merged['weight'] = merged.set_geometry("shared_boundary").length / total_boundary_length
+    else:
+        merged["weight"] = merged.set_geometry("shared_boundary").length
+    merged_with_islands = pd.concat((merged, islands))
+    length_weighted_w = W.from_adjlist(
+        merged_with_islands[["focal", "neighbor", "weight"]]
+    )
+    for island in w.islands:
+        length_weighted_w.neighbors[island] = []
+        del length_weighted_w.weights[island]
+
+    length_weighted_w._reset()
+
+    return length_weighted_w

libpysal/weights/tests/test_perimeter.py

@@ -0,0 +1,56 @@
+import os
+import tempfile
+
+import unittest
+import pytest
+import numpy as np
+import geopandas as gpd
+from ..weights import W, WSP
+from ..user import build_lattice_shapefile
+from .. import util
+from ..contiguity import Rook, _return_length_weighted_w
+from ... import examples
+
+NPTA3E = np.testing.assert_array_almost_equal
+
+
+class TestPerimeter(unittest.TestCase):
+    def setUp(self):
+        shp = build_lattice_shapefile(3, 3, "tmp.shp")
+        gdf = gpd.read_file("tmp.shp")
+        dv = [0] * 3
+        dv.extend(list(range(1, 7)))
+        gdf["dv"] = dv
+        gdf = gdf.dissolve(by="dv")
+        self.w0 = Rook.from_dataframe(gdf, perimeter=True, perim_std=False)
+        self.gdf = gdf
+
+        # us case
+        usgdf = gpd.read_file(examples.get_path("us48.shp"))
+        usgdf.set_crs("epsg:4326", inplace=True)
+        usgdf.to_crs(usgdf.estimate_utm_crs(), inplace=True)
+        self.usgdf = usgdf
+        self.wus = Rook.from_dataframe(usgdf, perimeter=True)
+
+    def test_perimeter(self):
+        NPTA3E(self.w0.pct_nonzero, 40.81632653)
+
+    def test_return_length_weighted(self):
+        w1 = _return_length_weighted_w(self.w0, self.gdf, False)
+        NPTA3E(w1.pct_nonzero, 40.81632653)
+        self.assertEqual(w1.weights[0], [1, 1, 1])
+        self.assertEqual(w1.weights[2], [1, 1, 1, 1])
+
+    def test_return_length_weighted_us(self):
+        w1 = _return_length_weighted_w(self.wus, self.usgdf, False)
+        self.assertAlmostEqual(w1[0][7], 354625.0684908081)
+        self.assertAlmostEqual(w1[0][10],  605834.5010118643)
+        NPTA3E(w1[0][7], w1[7][0])
+        w1.transform = "r"
+        self.assertAlmostEqual(w1[0][7],  0.3692243585791264)
+        self.assertAlmostEqual(w1[7][0],  0.12891667056448083)
+        self.assertNotAlmostEquals(w1[0][7], w1[7][0])
+
+
+if __name__ == "__main__":
+    unittest.main()