def _area_tables_raster(source_df, target_df, raster_path, codes=[21, 22, 23, 24], force_crs_match=True):
    "\n    Construct area allocation and source-target correspondence tables according to a raster 'populated' areas\n\n    Parameters\n    ----------\n    source_df : geopandas.GeoDataFrame\n        geeodataframe with geometry column of polygon type\n    target_df : geopandas.GeoDataFrame\n        geodataframe with geometry column of polygon type\n    raster_path : str\n        the path to the associated raster image.\n    codes : list\n        list of integer code values that should be considered as 'populated'.\n        Since this draw inspiration using the National Land Cover Database (NLCD), the default is 21 (Developed, Open Space), 22 (Developed, Low Intensity), 23 (Developed, Medium Intensity) and 24 (Developed, High Intensity).\n        The description of each code can be found here: https://www.mrlc.gov/sites/default/files/metadata/landcover.html\n        Only taken into consideration for harmonization raster based.\n    force_crs_match : bool (default is True)\n        Whether the Coordinate Reference System (CRS) of the polygon will be reprojected to the CRS of the raster file.\n        It is recommended to let this argument as True.\n\n    Returns\n    -------\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\n    Notes\n    -----\n    The assumption is both dataframes have the same coordinate reference system.\n\n    Union geometry is a geometry formed by the intersection of a source geometry and a target geometry\n\n    SU Maps source geometry to union geometry, UT maps union geometry to target geometry\n\n\n\n    "
    if _check_crs(source_df, target_df):
        pass
    else:
        return None
    source_df = source_df.copy()
    target_df = target_df.copy()
    n_s = source_df.shape[0]
    n_t = target_df.shape[0]
    _left = np.arange(n_s)
    _right = np.arange(n_t)
    source_df.loc[:, '_left'] = _left
    target_df.loc[:, '_right'] = _right
    res_union_pre = gpd.overlay(source_df, target_df, how='union')
    warnings.warn('The CRS for the generated union will be set to be the same as source_df.')
    res_union_pre.crs = source_df.crs
    res_union = _fast_append_profile_in_gdf(res_union_pre, raster_path, force_crs_match=force_crs_match)
    str_codes = [str(i) for i in codes]
    str_list = [('Type_' + i) for i in str_codes]
    str_list_ok = [col for col in res_union.columns if (col in str_list)]
    res_union['Populated_Pixels'] = res_union[str_list_ok].sum(axis=1)
    (n_u, _) = res_union.shape
    SU = np.zeros((n_s, n_u))
    UT = np.zeros((n_u, n_t))
    for (index, row) in res_union.iterrows():
        if ((not np.isnan(row['_left'])) and (not np.isnan(row['_right']))):
            s_id = int(row['_left'])
            t_id = int(row['_right'])
            SU[(s_id, index)] = row['Populated_Pixels']
            UT[(index, t_id)] = 1
    source_df.drop(['_left'], axis=1, inplace=True)
    target_df.drop(['_right'], axis=1, inplace=True)
    return (SU, UT)
Construct area allocation and source-target correspondence tables according to a raster 'populated' areas

Parameters
----------
source_df : geopandas.GeoDataFrame
    geeodataframe with geometry column of polygon type
target_df : geopandas.GeoDataFrame
    geodataframe with geometry column of polygon type
raster_path : str
    the path to the associated raster image.
codes : list
    list of integer code values that should be considered as 'populated'.
    Since this draw inspiration using the National Land Cover Database (NLCD), the default is 21 (Developed, Open Space), 22 (Developed, Low Intensity), 23 (Developed, Medium Intensity) and 24 (Developed, High Intensity).
    The description of each code can be found here: https://www.mrlc.gov/sites/default/files/metadata/landcover.html
    Only taken into consideration for harmonization raster based.
force_crs_match : bool (default is True)
    Whether the Coordinate Reference System (CRS) of the polygon will be reprojected to the CRS of the raster file.
    It is recommended to let this argument as True.

Returns
-------
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

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

Union geometry is a geometry formed by the intersection of a source geometry and a target geometry

SU Maps source geometry to union geometry, UT maps union geometry to target geometry
def _area_tables_raster(source_df, target_df, raster_path, codes=[21, 22, 23, 24], force_crs_match=True):
    "\n    Construct area allocation and source-target correspondence tables according to a raster 'populated' areas\n\n    Parameters\n    ----------\n    source_df : geopandas.GeoDataFrame\n        geeodataframe with geometry column of polygon type\n    target_df : geopandas.GeoDataFrame\n        geodataframe with geometry column of polygon type\n    raster_path : str\n        the path to the associated raster image.\n    codes : list\n        list of integer code values that should be considered as 'populated'.\n        Since this draw inspiration using the National Land Cover Database (NLCD), the default is 21 (Developed, Open Space), 22 (Developed, Low Intensity), 23 (Developed, Medium Intensity) and 24 (Developed, High Intensity).\n        The description of each code can be found here: https://www.mrlc.gov/sites/default/files/metadata/landcover.html\n        Only taken into consideration for harmonization raster based.\n    force_crs_match : bool (default is True)\n        Whether the Coordinate Reference System (CRS) of the polygon will be reprojected to the CRS of the raster file.\n        It is recommended to let this argument as True.\n\n    Returns\n    -------\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\n    Notes\n    -----\n    The assumption is both dataframes have the same coordinate reference system.\n\n    Union geometry is a geometry formed by the intersection of a source geometry and a target geometry\n\n    SU Maps source geometry to union geometry, UT maps union geometry to target geometry\n\n\n\n    "
    if _check_crs(source_df, target_df):
        pass
    else:
        return None
    source_df = source_df.copy()
    target_df = target_df.copy()
    n_s = source_df.shape[0]
    n_t = target_df.shape[0]
    _left = np.arange(n_s)
    _right = np.arange(n_t)
    source_df.loc[:, '_left'] = _left
    target_df.loc[:, '_right'] = _right
    res_union_pre = gpd.overlay(source_df, target_df, how='union')
    warnings.warn('The CRS for the generated union will be set to be the same as source_df.')
    res_union_pre.crs = source_df.crs
    res_union = _fast_append_profile_in_gdf(res_union_pre, raster_path, force_crs_match=force_crs_match)
    str_codes = [str(i) for i in codes]
    str_list = [('Type_' + i) for i in str_codes]
    str_list_ok = [col for col in res_union.columns if (col in str_list)]
    res_union['Populated_Pixels'] = res_union[str_list_ok].sum(axis=1)
    (n_u, _) = res_union.shape
    SU = np.zeros((n_s, n_u))
    UT = np.zeros((n_u, n_t))
    for (index, row) in res_union.iterrows():
        if ((not np.isnan(row['_left'])) and (not np.isnan(row['_right']))):
            s_id = int(row['_left'])
            t_id = int(row['_right'])
            SU[(s_id, index)] = row['Populated_Pixels']
            UT[(index, t_id)] = 1
    source_df.drop(['_left'], axis=1, inplace=True)
    target_df.drop(['_right'], axis=1, inplace=True)
    return (SU, UT)