def create_intersection_matrix(gdf_intersection,
                               gdf_intersection_col="gdf_intermediate_intersection_fraction",
                               gdf_intermediate_col="FIPS",
                               gdf_source_col="ID",
                               normalization=True,
                               normalize_axis=1):
    """
    Pivots stacked dataframe to an intersection matrix.

    Inputs the intersection matrix with the area intersection fraction a[i, j]
    for every intermediate ID i and source ID j. Normalization option available
    (To normalize double-counting)

    Args:
        gdf_intersection (pandas.DataFrame): stacked dataframe with every
            one-to-one mapping of intermediate ids.
        gdf_intersection_col (str): name of the column of the area intersection
            fraction. (Usually the fraction along the column you want to
            normalize. Generally makes sense to do it along the distinct and
            disjoint column like census tracts)
        gdf_intermediate_col (str): ID name of intermediate gdf
        gdf_source_col (str): ID name of source gdf
        normalize_axis (int): 1 for normalizing along intermediate axis
            (generally); 0 for normalizing along source axis

    Returns:
        pandas.DataFrame: matrix which consists the intersection value for
        every intermediate geometry i and source geometry j at index i, j

    """
    intersection_matrix = gdf_intersection.pivot_table(
        values=gdf_intersection_col,
        index=gdf_intermediate_col,
        columns=gdf_source_col,
        fill_value=0
    )

    if normalization is True:
        intersection_matrix = intersection_matrix.divide(intersection_matrix.sum(axis=normalize_axis),
                                                         axis=int(not normalize_axis))

    return intersection_matrix