def _area_interpolate(source_df, target_df, extensive_variables=None, intensive_variables=None, tables=None, allocate_total=True):
    '\n    Area interpolation for extensive and intensive variables.\n\n    Parameters\n    ----------\n    source_df : geopandas.GeoDataFrame (required)\n        geodataframe with polygon geometries\n    target_df : geopandas.GeoDataFrame (required)\n        geodataframe with polygon geometries\n    extensive_variables : list, (optional)\n        columns in dataframes for extensive variables\n    intensive_variables : list, (optional)\n        columns in dataframes for intensive variables\n    tables : tuple (optional)\n        two 2-D numpy arrays\n        SU: area of intersection of source geometry i with union geometry j\n        UT: binary mapping of union geometry j to target geometry t\n    allocate_total : boolean\n        True if total value of source area should be allocated.\n        False if denominator is area of i. Note that the two cases\n        would be identical when the area of the source polygon is\n        exhausted by intersections. See Notes for more details.\n\n    Returns\n    -------\n    estimates : geopandas.GeoDataFrame\n        new geodaraframe with interpolated variables as columns and target_df geometry\n        as output geometry\n\n    Notes\n    -----\n    The assumption is both dataframes have the same coordinate reference system.\n\n\n    For an extensive variable, the estimate at target polygon j (default case) is:\n\n    v_j = \\sum_i v_i w_{i,j}\n\n    w_{i,j} = a_{i,j} / \\sum_k a_{i,k}\n\n\n    If the area of the source polygon is not exhausted by intersections with\n    target polygons and there is reason to not allocate the complete value of\n    an extensive attribute, then setting allocate_total=False will use the\n    following weights:\n\n\n    v_j = \\sum_i v_i w_{i,j}\n\n    w_{i,j} = a_{i,j} / a_i\n\n    where a_i is the total area of source polygon i.\n\n\n    For an intensive variable, the estimate at target polygon j is:\n\n    v_j = \\sum_i v_i w_{i,j}\n\n    w_{i,j} = a_{i,j} / \\sum_k a_{k,j}\n\n    '
    source_df = source_df.copy()
    target_df = target_df.copy()
    if _check_crs(source_df, target_df):
        pass
    else:
        return None
    if (tables is None):
        (SU, UT) = _area_tables(source_df, target_df)
    else:
        (SU, UT) = tables
    den = source_df[source_df.geometry.name].area.values
    if allocate_total:
        den = SU.sum(axis=1)
    den = (den + (den == 0))
    weights = np.dot(np.diag((1 / den)), SU)
    dfs = []
    extensive = []
    if extensive_variables:
        for variable in extensive_variables:
            vals = _nan_check(source_df, variable)
            vals = _inf_check(source_df, variable)
            estimates = np.dot(np.diag(vals), weights)
            estimates = np.dot(estimates, UT)
            estimates = estimates.sum(axis=0)
            extensive.append(estimates)
        extensive = np.array(extensive)
        extensive = pd.DataFrame(extensive.T, columns=extensive_variables)
    ST = np.dot(SU, UT)
    area = ST.sum(axis=0)
    den = np.diag((1.0 / (area + (area == 0))))
    weights = np.dot(ST, den)
    intensive = []
    if intensive_variables:
        for variable in intensive_variables:
            vals = _nan_check(source_df, variable)
            vals = _inf_check(source_df, variable)
            vals.shape = (len(vals), 1)
            est = (vals * weights).sum(axis=0)
            intensive.append(est)
        intensive = np.array(intensive)
        intensive = pd.DataFrame(intensive.T, columns=intensive_variables)
    if extensive_variables:
        dfs.append(extensive)
    if intensive_variables:
        dfs.append(intensive)
    df = pd.concat(dfs, axis=1)
    df['geometry'] = target_df[target_df.geometry.name].reset_index(drop=True)
    df = gpd.GeoDataFrame(df.replace(np.inf, np.nan))
    return df
Area interpolation for extensive and intensive variables.

Parameters
----------
source_df : geopandas.GeoDataFrame (required)
    geodataframe with polygon geometries
target_df : geopandas.GeoDataFrame (required)
    geodataframe with polygon geometries
extensive_variables : list, (optional)
    columns in dataframes for extensive variables
intensive_variables : list, (optional)
    columns in dataframes for intensive variables
tables : tuple (optional)
    two 2-D numpy arrays
    SU: area of intersection of source geometry i with union geometry j
    UT: binary mapping of union geometry j to target geometry t
allocate_total : boolean
    True if total value of source area should be allocated.
    False if denominator is area of i. Note that the two cases
    would be identical when the area of the source polygon is
    exhausted by intersections. See Notes for more details.

Returns
-------
estimates : geopandas.GeoDataFrame
    new geodaraframe with interpolated variables as columns and target_df geometry
    as output geometry

Notes
-----
The assumption is both dataframes have the same coordinate reference system.


For an extensive variable, the estimate at target polygon j (default case) is:

v_j = \sum_i v_i w_{i,j}

w_{i,j} = a_{i,j} / \sum_k a_{i,k}


If the area of the source polygon is not exhausted by intersections with
target polygons and there is reason to not allocate the complete value of
an extensive attribute, then setting allocate_total=False will use the
following weights:


v_j = \sum_i v_i w_{i,j}

w_{i,j} = a_{i,j} / a_i

where a_i is the total area of source polygon i.


For an intensive variable, the estimate at target polygon j is:

v_j = \sum_i v_i w_{i,j}

w_{i,j} = a_{i,j} / \sum_k a_{k,j}
def _area_interpolate(source_df, target_df, extensive_variables=None, intensive_variables=None, tables=None, allocate_total=True):
    '\n    Area interpolation for extensive and intensive variables.\n\n    Parameters\n    ----------\n    source_df : geopandas.GeoDataFrame (required)\n        geodataframe with polygon geometries\n    target_df : geopandas.GeoDataFrame (required)\n        geodataframe with polygon geometries\n    extensive_variables : list, (optional)\n        columns in dataframes for extensive variables\n    intensive_variables : list, (optional)\n        columns in dataframes for intensive variables\n    tables : tuple (optional)\n        two 2-D numpy arrays\n        SU: area of intersection of source geometry i with union geometry j\n        UT: binary mapping of union geometry j to target geometry t\n    allocate_total : boolean\n        True if total value of source area should be allocated.\n        False if denominator is area of i. Note that the two cases\n        would be identical when the area of the source polygon is\n        exhausted by intersections. See Notes for more details.\n\n    Returns\n    -------\n    estimates : geopandas.GeoDataFrame\n        new geodaraframe with interpolated variables as columns and target_df geometry\n        as output geometry\n\n    Notes\n    -----\n    The assumption is both dataframes have the same coordinate reference system.\n\n\n    For an extensive variable, the estimate at target polygon j (default case) is:\n\n    v_j = \\sum_i v_i w_{i,j}\n\n    w_{i,j} = a_{i,j} / \\sum_k a_{i,k}\n\n\n    If the area of the source polygon is not exhausted by intersections with\n    target polygons and there is reason to not allocate the complete value of\n    an extensive attribute, then setting allocate_total=False will use the\n    following weights:\n\n\n    v_j = \\sum_i v_i w_{i,j}\n\n    w_{i,j} = a_{i,j} / a_i\n\n    where a_i is the total area of source polygon i.\n\n\n    For an intensive variable, the estimate at target polygon j is:\n\n    v_j = \\sum_i v_i w_{i,j}\n\n    w_{i,j} = a_{i,j} / \\sum_k a_{k,j}\n\n    '
    source_df = source_df.copy()
    target_df = target_df.copy()
    if _check_crs(source_df, target_df):
        pass
    else:
        return None
    if (tables is None):
        (SU, UT) = _area_tables(source_df, target_df)
    else:
        (SU, UT) = tables
    den = source_df[source_df.geometry.name].area.values
    if allocate_total:
        den = SU.sum(axis=1)
    den = (den + (den == 0))
    weights = np.dot(np.diag((1 / den)), SU)
    dfs = []
    extensive = []
    if extensive_variables:
        for variable in extensive_variables:
            vals = _nan_check(source_df, variable)
            vals = _inf_check(source_df, variable)
            estimates = np.dot(np.diag(vals), weights)
            estimates = np.dot(estimates, UT)
            estimates = estimates.sum(axis=0)
            extensive.append(estimates)
        extensive = np.array(extensive)
        extensive = pd.DataFrame(extensive.T, columns=extensive_variables)
    ST = np.dot(SU, UT)
    area = ST.sum(axis=0)
    den = np.diag((1.0 / (area + (area == 0))))
    weights = np.dot(ST, den)
    intensive = []
    if intensive_variables:
        for variable in intensive_variables:
            vals = _nan_check(source_df, variable)
            vals = _inf_check(source_df, variable)
            vals.shape = (len(vals), 1)
            est = (vals * weights).sum(axis=0)
            intensive.append(est)
        intensive = np.array(intensive)
        intensive = pd.DataFrame(intensive.T, columns=intensive_variables)
    if extensive_variables:
        dfs.append(extensive)
    if intensive_variables:
        dfs.append(intensive)
    df = pd.concat(dfs, axis=1)
    df['geometry'] = target_df[target_df.geometry.name].reset_index(drop=True)
    df = gpd.GeoDataFrame(df.replace(np.inf, np.nan))
    return df