Interpolate to the union of the polygons from two dataframes

tobler/area_weighted/__init__.py

@@ -2,3 +2,5 @@
 from .area_interpolate import _area_tables_binning
 from .area_join import area_join
 from .area_interpolate_dask import area_interpolate_dask
+from .area_faces import area_faces
+from .area_buffer import area_buffer

tobler/area_weighted/area_buffer.py

@@ -0,0 +1,40 @@
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from .area_interpolate import _area_interpolate_binning as area_interpolate
+
+__author__ = "Serge Rey <[email protected]>"
+
+
+def area_buffer(source_df, buffer_df, in_place=False):
+    """
+    Classify spatial relationship of source_df geometries relative to buffer_df geometries
+
+    Parameters
+    ----------
+    source_df : geopandas.GeoDataFrame
+        GeoDataFrame containing source values
+    buffer_df : geopandas.GeoDataFrame
+        GeoDataFrame containing buffer geometries
+    in_place : boolean
+        If True, the source_df is modified in place, otherwise a copy
+    is returned (default)
+
+    Returns
+    -------
+    source_df : geopandas.GeoDataFrame
+        GeoDataFrame with additional column `right_relation` that
+        takes three possible values ['within', 'partial', 'disjoint']
+        specifying the spatial predicate of source to buffer
+        geometries.
+
+    """
+    within = buffer_df.sindex.query(source_df.geometry, predicate="within")[0]
+    intersects = buffer_df.sindex.query(source_df.geometry, predicate="intersects")[0]
+    partial = [i for i in intersects if i not in within]
+    if not in_place:
+        source_df = source_df.copy()
+    source_df['right_relation'] = 'disjoint'
+    source_df.loc[partial, 'right_relation'] = 'partial'
+    source_df.loc[within, 'right_relation'] = 'within'
+    return source_df

tobler/area_weighted/area_faces.py

@@ -0,0 +1,39 @@
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from .area_interpolate import _area_interpolate_binning as area_interpolate
+
+__author__ = "Serge Rey <[email protected]>"
+
+
+def area_faces(source_df, target_df,
+               extensive_variables=[],
+               intensive_variables=[]):
+    """
+    Interpolation of source_df values to faces formed by the union of
+               the source and target dataframes.
+
+
+    Parameters
+    ----------
+    source_df : geopandas.GeoDataFrame
+        GeoDataFrame containing source values
+    target_df : geopandas.GeoDataFrame
+        GeoDataFrame containing target values
+    extensive_variables : string or list-like
+        column(s) in source_df dataframe for extensive variable(s) to
+    be interpolated
+    intensive_variables : string or list-like
+        column(s) in source_df dataframe for intensive variable(s) to
+    be interpolated
+
+    Returns
+    -------
+    results : geopandas.GeoDataFrame
+        GeoDataFrame with interpolated values as additional columns
+
+    """
+
+    union = gpd.overlay(source_df, target_df, how="union")
+    results = area_interpolate(source_df, union, extensive_variables, intensive_variables)
+    return results

tobler/tests/test_area_faces.py

@@ -0,0 +1,54 @@
+import geopandas as gpd
+import numpy as np
+from shapely.geometry import Polygon
+
+import pytest
+
+from tobler.area_weighted import area_faces, area_buffer
+
+
+class TestAreaFaces:
+    def setup_method(self):
+        polys1 = gpd.GeoSeries([Polygon([(0,0), (10,0), (10,5), (0,5)]),
+                            Polygon([(0,5), (0,10),  (5,10), (5,5)]),
+                            Polygon([(5,5), (5,10), (7,10), (7,5)]),
+                            Polygon([(7,5), (7,10), (10,10), (10,5)]) ]
+                            )
+
+
+        buffer = gpd.GeoSeries([Polygon([ (0,0), (0, 10), (6, 10), (6,0)])])
+
+
+        df1 = gpd.GeoDataFrame({'geometry': polys1})
+        df2 = gpd.GeoDataFrame({'geometry': buffer})
+        df1['population'] = [ 500,  200, 100, 50]
+        df1['pci'] = [75, 100, 40, 30]
+        df1['income'] = df1['population'] * df1['pci']
+        df2['population'] = 10000
+        df2['pci'] = 80
+        self.source = df1
+        self.target = df2
+
+
+    def test_area_faces(self):
+        result = area_faces(self.source, self.target,
+                            extensive_variables=['population']) 
+        assert (result.shape == (6,2))
+        pop_values = np.array([299.99998212, 200., 50., 199.99998808, 50.  , 50.  ])
+        np.testing.assert_almost_equal(result.population.values, pop_values, 2)
+
+    def test_area_faces_2_1(self):
+        result = area_faces(self.target, self.source,
+                            extensive_variables=['population'],
+                            intensive_variables=['pci']) 
+        assert (result.shape == (6,3))
+        pop_values = np.array([5000, 4166.67, 833.33, 0, 0, 0 ])
+        np.testing.assert_almost_equal(result.population.values, pop_values, 2)
+        pci_values = np.array([80, 80, 80, 0, 0, 0 ])
+        np.testing.assert_almost_equal(result.pci.values, pci_values, 2)
+
+    def test_area_buffer(self):
+        result = area_buffer(self.source, self.target)
+        preds = ['partial', 'within', 'partial', 'disjoint']
+        assert (result.right_relation.tolist() == preds)
+