Repository: pysal/tobler (language: python)

Issue:
Pick attribute based on largest intersection
I have a function which links variables from one df to the other based on the largest intersection (which of the geometries in `source_df` has the largest area when intersected with geometry in `target_df`). Though it is not as fancy as the rest of Tobler, it may be useful in some cases (like when you have _almost_ the same geometries in both - like tessellation vs parcels). 

Do you think it would make sense to have it in Tobler? Happy to do a PR. 

cc @darribas


```python
def area_max(source_df, target_df, variables):
    """
    Join attributes from source based on the largest intersection. In case of a tie it picks the first one.
    
    Parameters
    ----------
    source_df : GeoDataFrame
        GeoDataFrame containing source values
    target_df : GeoDataFrame
        GeoDataFrame containing source values
    variables : string or list-like
        column(s) in dataframes for variable(s)
    
    Returns
    -------
    GeoDataFrame
    
    Notes
    -----
    TODO: preserve dtype where possible
    
    """    
    target_df = target_df.copy()
    target_ix, source_ix = source_df.sindex.query_bulk(target_df.geometry, predicate='intersects')
    areas = target_df.geometry.values[target_ix].intersection(source_df.geometry.values[source_ix]).area

    main = []
    for i in range(len(target_df)):
        mask = target_ix == i
        if np.any(mask):
            main.append(source_ix[mask][np.argmax(areas[mask])])
        else:
            main.append(np.nan)
    
    main = np.array(main)
    mask = ~np.isnan(main)
    if pd.api.types.is_list_like(variables):
        for v in variables:
            arr = np.empty(len(main), dtype=object)
            arr[:] = np.nan
            arr[mask] = source_df[v].values[main[mask].astype(int)]
            target_df[v] = arr
    else:
        arr = np.empty(len(main), dtype=object)
        arr[:] = np.nan
        arr[mask] = source_df[variables].values[main[mask].astype(int)]
        target_df[variables] = arr
        
    return target_df
```


Discussion / hints:
i can also imagine an analogous case where you have two datasets digitized from different sources, where one version is 'correct' but the other has useful data that you need to transfer to the other polys without needing to reapportion it. (I remember hand-digitizing some school boundary maps and manually entering some data points when i worked at a non-profit a few years ago, and now that data could be transferred to the new canonical boundaries released through other channels). Right now this is often done crudely through, e.g. the `representative_point`, but this max area method would often be preferable in these situations

all that to say i see the value in this and i think it makes sense in `tobler`, though i'm not 100% on where exactly, and what it should be called 
I think this would nicely expand the functionality available to the users. So +1.

Fix patch:
diff --git a/docs/api.rst b/docs/api.rst
index 9931612..4ed2775 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -9,7 +9,7 @@ API reference
 
 Area Weighted
 --------------
-Area weighted approaches use the area of overlap between the source and target geometries to weight the 
+Area weighted approaches use the area of overlap between the source and target geometries to weight the
 variables being assigned to the target
 
 .. currentmodule:: tobler.area_weighted
@@ -18,12 +18,13 @@ variables being assigned to the target
    :toctree: generated/
 
     area_interpolate
+    area_join
 
 
 Dasymetric
 -----------------------
 
-Dasymetric approaches use auxiliary data in addition to use the area of overlap between the source and target geometries to weight the 
+Dasymetric approaches use auxiliary data in addition to use the area of overlap between the source and target geometries to weight the
 variables being assigned to the target
 
 .. currentmodule:: tobler.dasymetric
diff --git a/tobler/area_weighted/__init__.py b/tobler/area_weighted/__init__.py
index dea94ee..93cf55e 100644
--- a/tobler/area_weighted/__init__.py
+++ b/tobler/area_weighted/__init__.py
@@ -1,4 +1,5 @@
 from .area_interpolate import _area_interpolate_binning as area_interpolate
 from .area_interpolate import _area_interpolate as _slow_area_interpolate
 from .area_interpolate import _area_tables, _area_tables_binning, _area_tables_raster
-from .area_interpolate import _check_presence_of_crs
\ No newline at end of file
+from .area_interpolate import _check_presence_of_crs
+from .area_join import area_join
\ No newline at end of file
diff --git a/tobler/area_weighted/area_join.py b/tobler/area_weighted/area_join.py
new file mode 100644
index 0000000..0d56744
--- /dev/null
+++ b/tobler/area_weighted/area_join.py
@@ -0,0 +1,68 @@
+import numpy as np
+import pandas as pd
+import warnings
+
+__author__ = "Martin Fleischmann <martin@martinfleischmann.net>"
+
+
+def area_join(source_df, target_df, variables):
+    """
+    Join variables from source_df based on the largest intersection. In case of a tie it picks the first one.
+
+    Parameters
+    ----------
+    source_df : geopandas.GeoDataFrame
+        GeoDataFrame containing source values
+    target_df : geopandas.GeoDataFrame
+        GeoDataFrame containing source values
+    variables : string or list-like
+        column(s) in source_df dataframe for variable(s) to be joined
+
+    Returns
+    -------
+    joined : geopandas.GeoDataFrame
+         target_df GeoDataFrame with joined variables as additional columns
+
+
+
+    """
+    if not pd.api.types.is_list_like(variables):
+        variables = [variables]
+
+    for v in variables:
+        if v in target_df.columns:
+            raise ValueError(f"Column '{v}' already present in target_df.")
+
+    target_df = target_df.copy()
+    target_ix, source_ix = source_df.sindex.query_bulk(
+        target_df.geometry, predicate="intersects"
+    )
+    areas = (
+        target_df.geometry.values[target_ix]
+        .intersection(source_df.geometry.values[source_ix])
+        .area
+    )
+
+    main = []
+    for i in range(len(target_df)):  # vectorise this loop?
+        mask = target_ix == i
+        if np.any(mask):
+            main.append(source_ix[mask][np.argmax(areas[mask])])
+        else:
+            main.append(np.nan)
+
+    main = np.array(main, dtype=float)
+    mask = ~np.isnan(main)
+
+    for v in variables:
+        arr = np.empty(len(main), dtype=object)
+        arr[mask] = source_df[v].values[main[mask].astype(int)]
+        try:
+            arr = arr.astype(source_df[v].dtype)
+        except TypeError:
+            warnings.warn(
+                f"Cannot preserve dtype of '{v}'. Falling back to `dtype=object`.",
+            )
+        target_df[v] = arr
+
+    return target_df


Test patch:
diff --git a/tobler/tests/test_area_join.py b/tobler/tests/test_area_join.py
new file mode 100644
index 0000000..d3d2574
--- /dev/null
+++ b/tobler/tests/test_area_join.py
@@ -0,0 +1,70 @@
+import geopandas as gpd
+import numpy as np
+from shapely.geometry import Point
+
+import pytest
+
+from tobler.area_weighted import area_join
+
+
+class TestAreaJoin:
+    def setup_method(self):
+        self.grid = gpd.points_from_xy(
+            np.repeat(np.linspace(1, 10, 10), 10), np.tile(np.linspace(1, 10, 10), 10)
+        ).buffer(0.5, cap_style=3)
+        self.source = gpd.GeoDataFrame(
+            {
+                "floats": np.linspace(1, 10, 100),
+                "ints": np.linspace(1, 100, 100, dtype="int"),
+                "strings": np.array(["darribas", "is", "the", "king"] * 25),
+            },
+            geometry=self.grid,
+        )
+
+        self.target = gpd.GeoDataFrame(geometry=self.grid.translate(xoff=2.2, yoff=0.2))
+
+    def test_area_join_float(self):
+        result = area_join(self.source, self.target, "floats")
+        assert (result.columns == ["geometry", "floats"]).all()
+        np.testing.assert_almost_equal(result.floats.mean(), 6.409, 3)
+        assert result.floats.dtype == float
+        assert result.floats.isna().sum() == 20
+
+    def test_area_join_ints(self):
+        with pytest.warns(UserWarning, match="Cannot preserve dtype of"):
+            result = area_join(self.source, self.target, "ints")
+
+        assert (result.columns == ["geometry", "ints"]).all()
+        np.testing.assert_almost_equal(result.ints.mean(), 60.5, 3)
+        assert result.ints.dtype == object
+        assert type(result.ints.iloc[0]) == int
+        assert result.ints.isna().sum() == 20
+
+    def test_area_join_strings(self):
+        result = area_join(self.source, self.target, "strings")
+        assert (result.columns == ["geometry", "strings"]).all()
+        assert result.strings.dtype == object
+        assert type(result.strings.iloc[0]) == str
+        assert result.strings.isna().sum() == 20
+
+    def test_area_join_array(self):
+        with pytest.warns(UserWarning, match="Cannot preserve dtype of"):
+            result = area_join(self.source, self.target, ["floats", "ints", "strings"])
+
+        assert (result.columns == ["geometry", "floats", "ints", "strings"]).all()
+        np.testing.assert_almost_equal(result.floats.mean(), 6.409, 3)
+        assert result.floats.dtype == float
+        assert result.floats.isna().sum() == 20
+        np.testing.assert_almost_equal(result.ints.mean(), 60.5, 3)
+        assert result.ints.dtype == object
+        assert type(result.ints.iloc[0]) == int
+        assert result.ints.isna().sum() == 20
+        assert result.strings.dtype == object
+        assert type(result.strings.iloc[0]) == str
+        assert result.strings.isna().sum() == 20
+
+    def test_area_join_error(self):
+        target = self.target
+        target["floats"] = 0
+        with pytest.raises(ValueError, match="Column 'floats'"):
+            area_join(self.source, target, "floats")
\ No newline at end of file
diff --git a/tobler/tests/test_dasymetric.py b/tobler/tests/test_dasymetric.py
index d7561b6..57fdb17 100644
--- a/tobler/tests/test_dasymetric.py
+++ b/tobler/tests/test_dasymetric.py
@@ -2,8 +2,10 @@
 import sys
 import pandas as pd
 import geopandas
+
 try:
     import quilt3
+
     QUILTMISSING = False
 except ImportError:
     QUILTMISSING = True
@@ -23,27 +25,30 @@ def datasets():
         if not os.path.exists("nlcd_2011.tif"):
             p = quilt3.Package.browse("rasters/nlcd", "s3://spatial-ucr")
             p["nlcd_2011.tif"].fetch()
-        sac1 = load_example('Sacramento1')
-        sac2 = load_example('Sacramento2')
+        sac1 = load_example("Sacramento1")
+        sac2 = load_example("Sacramento2")
         sac1 = geopandas.read_file(sac1.get_path("sacramentot2.shp"))
         sac2 = geopandas.read_file(sac2.get_path("SacramentoMSA2.shp"))
-        sac1['pct_poverty'] = sac1.POV_POP/sac1.POV_TOT
+        sac1["pct_poverty"] = sac1.POV_POP / sac1.POV_TOT
 
         return sac1, sac2
     else:
         pass
 
+
 @pytest.mark.skipif(QUILTMISSING, reason="quilt3 not available.")
 def test_area_interpolate():
     sac1, sac2 = datasets()
     area = area_interpolate(
-        source_df=sac1, target_df=sac2, extensive_variables=["TOT_POP"], intensive_variables=["pct_poverty"]
+        source_df=sac1,
+        target_df=sac2,
+        extensive_variables=["TOT_POP"],
+        intensive_variables=["pct_poverty"],
     )
     assert_almost_equal(area.TOT_POP.sum(), 1796856, decimal=0)
     assert_almost_equal(area.pct_poverty.sum(), 2140, decimal=0)
 
 
-
 @pytest.mark.skipif(QUILTMISSING, reason="quilt3 not available.")
 def test_masked_area_interpolate():
     sac1, sac2 = datasets()
@@ -51,7 +56,7 @@ def test_masked_area_interpolate():
         source_df=sac1,
         target_df=sac2,
         extensive_variables=["TOT_POP"],
-        intensive_variables=['pct_poverty'],
+        intensive_variables=["pct_poverty"],
         raster="nlcd_2011.tif",
     )
     assert masked.TOT_POP.sum() > 1500000