Skip to content
Snippets Groups Projects
FireClustering.py 8.09 KiB
Newer Older
screbec's avatar
screbec committed
""" FireClustering
This module include all functions used for doing fire clustering
"""

Yang's avatar
Yang committed

screbec's avatar
screbec committed
def remove_self(inds, dist):
Yang's avatar
Yang committed
    """ Remove self from the index and distance arrays
screbec's avatar
screbec committed

    Parameters
    ----------
    inds : np array of np array
        indices of neighbors for each point (self included)
    dist : np array of np array
        distance of neighbors for each point (self included)

    Returns
    -------
    new_inds : np array of np array
        indices of neighbors (self excluded)
    new_dist : np array of np array
        distance of neighbors (self excluded)
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    import numpy as np

    new_inds = []
    new_dist = []
    for i in range(len(inds)):
Yang's avatar
Yang committed
        pos = np.where(inds[i] == i)
        new_inds.append(np.delete(inds[i], pos))
        new_dist.append(np.delete(dist[i], pos))

    return np.array(new_inds, dtype=object), np.array(new_dist, dtype=object)
screbec's avatar
screbec committed


def build_rtree(geoms, fids=False):
Yang's avatar
Yang committed
    """Builds Rtree from a shapely multipolygon shape
    and optionally uses list of fids as identifier"""
screbec's avatar
screbec committed
    import rtree
Yang's avatar
Yang committed

Yang's avatar
Yang committed
    idx = rtree.index.Index()  # create new index
screbec's avatar
screbec committed
    for ind, geom in enumerate(geoms):
        if fids:
            idx.insert(ind, geom.bounds, fids[ind])
        else:
            idx.insert(ind, geom.bounds, ind)
Yang's avatar
Yang committed

screbec's avatar
screbec committed
    return idx

def idx_intersection(idx, bbox):
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    Finds all objects in an index whcih bounding boxes intersect with a geometry's bounding box
Yang's avatar
Yang committed
    """
    intersections = list(idx.intersection(bbox, objects=True))
    if len(intersections) > 0:
screbec's avatar
screbec committed
        fids, bbox = zip(*[(item.object, item.bbox) for item in intersections])
    else:
        fids = []
    return fids


def compute_all_spatial_distances(data, max_thresh_km):
Yang's avatar
Yang committed
    """ Derive neighbors (and distances) for each point (with lat/lon)
screbec's avatar
screbec committed

    Parameters
    ----------
    data : pd DataFrame
Yang's avatar
Yang committed
        point location with 'x' and 'y' columns
screbec's avatar
screbec committed
    max_thresh_km : float
        maximum distance threshold (km) used for classifying neighbors

    Returns
    -------
    inds : np array of np array
        indices of neighbors (self exluded)
    dist : np array of np array
        distance (km) of neighbors (self excluded)
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    import numpy as np
    from sklearn.neighbors import BallTree

Yang's avatar
Yang committed
    x = data.x.values
    y = data.y.values
screbec's avatar
screbec committed

Yang's avatar
Yang committed
    X = np.stack([x, y], axis=-1)
screbec's avatar
screbec committed

Yang's avatar
Yang committed
    bt = BallTree(X, leaf_size=20)
screbec's avatar
screbec committed

Yang's avatar
Yang committed
    inds, dist = bt.query_radius(X, r=max_thresh_km * 1000, return_distance=True)
screbec's avatar
screbec committed
    inds, dist = remove_self(inds, dist)
Yang's avatar
Yang committed
    return inds, dist * 1000
screbec's avatar
screbec committed

screbec's avatar
screbec committed
def sort_neighbors(inds, dists):
Yang's avatar
Yang committed
    """ Do neighbor sorting (based on distance) for all points
screbec's avatar
screbec committed

    Parameters
    ----------
    inds : np array of np array
        indices of neighbors (self exluded)
    dist : np array of np array
        distance (km) of neighbors (self excluded)

    Returns
    -------
    sorted_inds : np array of np array
        indices of neighbors (self exluded)
    sorted_dists : np array of np array
        distance (km) of neighbors (self excluded)
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    import numpy as np

    sorted_inds = []
    sorted_dists = []

    for i in range(len(inds)):
        new_order = np.argsort(dists[i])

        sorted_inds.append(inds[i][new_order])
        sorted_dists.append(dists[i][new_order])

    return sorted_inds, sorted_dists

screbec's avatar
screbec committed
def do_clustering(data, max_thresh_km):
Yang's avatar
Yang committed
    """ Do initial clustering for fire pixels
screbec's avatar
screbec committed

    Parameters
    ----------
    data : list (nx2)
        latitude and longitude values of all fire pixels
    max_thresh_km : float
        maximum distance threshold (km) used for classifying neighbors

    Returns
    -------
    point_to_cluster_id : list
        cluster id for each fire point
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    import numpy as np
    import pandas as pd

    # value to fill in pixels without clustering
    NO_CLUSTER_VAL = -1

    # if number of points is 1 or 2, each point is one cluster
    num_points = len(data)
    if num_points < 3:
        cluster_id = list(range(num_points))
        return cluster_id

    # convert list to pd DataFrame
Yang's avatar
Yang committed
    dfdata = pd.DataFrame(data, columns=["x", "y"])
screbec's avatar
screbec committed

    # initialization
    cluster_id_counter = 0
    point_to_cluster_id = np.full(num_points, fill_value=NO_CLUSTER_VAL, dtype=np.int64)

    # compute and sort neighbor pixels for each pixel
Yang's avatar
Yang committed
    neighbor_inds, neighbor_spatial_dists = compute_all_spatial_distances(
        dfdata, max_thresh_km
    )
    # neighbor_inds, neighbor_spatial_dists = sort_neighbors(neighbor_inds, neighbor_spatial_dists)

    # include all possible pixels in cluster
    to_check = np.full(num_points, fill_value=1, dtype=np.int8)
    while np.sum(to_check) > 0:
Yang's avatar
Yang committed
        start_ind = np.argmax(to_check == 1)  # catch first index to check
Yang's avatar
Yang committed

        neighbors_to_search = list(neighbor_inds[start_ind])
        all_neighbors = neighbors_to_search
Yang's avatar
Yang committed

Yang's avatar
Yang committed
        if (
            len(all_neighbors) == 0
        ):  # if no neighbor, record the current pixel as a separate cluster
            point_to_cluster_id[start_ind] = cluster_id_counter
screbec's avatar
screbec committed
            cluster_id_counter += 1
            to_check[start_ind] = 0
Yang's avatar
Yang committed

screbec's avatar
screbec committed
        else:  # if with neighbor, record all neighbors
            # find all neighbours of neighbours:
            searched_neighbours = [start_ind]
            while len(neighbors_to_search) > 0:
                # take the first of these
                px = neighbors_to_search[0]
                searched_neighbours.append(px)
                px_neighbors = list(neighbor_inds[px])
                all_neighbors = list(set(all_neighbors + px_neighbors))
Yang's avatar
Yang committed
                neighbors_to_search = list(
                    set(all_neighbors).difference(searched_neighbours)
                )
            # now we have all pixels in this cluster in all_neighbors
            point_to_cluster_id[all_neighbors] = cluster_id_counter
            cluster_id_counter += 1
            to_check[all_neighbors] = 0
Yang's avatar
Yang committed

screbec's avatar
screbec committed
    return point_to_cluster_id.tolist()

Yang's avatar
Yang committed

def cal_distance(loc1, loc2):
    """ Calculate the distance between two points
screbec's avatar
screbec committed

    Parameters
    ----------
    loc1 : list[lat,lon]
        position of first point
    loc2 : list[lat,lon]
        position of second point

    Returns
    -------
    distance : float
        distance (km) between two points
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    from FireConsts import EARTH_RADIUS_KM
    import math

    lat1 = math.radians(loc1[0])
    lon1 = math.radians(loc1[1])
    lat2 = math.radians(loc2[0])
    lon2 = math.radians(loc2[1])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

Yang's avatar
Yang committed
    a = (
        math.sin(dlat / 2) ** 2
        + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    )
screbec's avatar
screbec committed
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))

    distance = EARTH_RADIUS_KM * c
    return distance

Yang's avatar
Yang committed

def cal_mindist(c1, c2):
    """ Calculate the minimum distance beween two clusters (may modify the algorithm to speed up this calculation)
screbec's avatar
screbec committed

    Parameters
    ----------
    c1 : list of [lat,lon]
        first cluster
    c2 : list of [lat,lon]
        second cluster

    Returns
    -------
    mindist : float
        the minimum distance (km) between c1 and c2
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed

    import itertools

Yang's avatar
Yang committed
    mindist = min([cal_distance(l1, l2) for l1, l2 in itertools.product(c1, c2)])
screbec's avatar
screbec committed

    return mindist

Yang's avatar
Yang committed
# def filter_centroid(loc_tgt,locs,MAX_THRESH_CENT_KM=50):
#     ''' Get cluster ids if centroid is close to the target cluster centroid
#
#     Parameters
#     ----------
#     loc_tgt : list, [lat,lon]
#         target centroid
#     locs : list of [lat,lon]
#         all existing fire centroid
#     MAX_THRESH_CENT_KM : float
#         maximum threshold to determine close clusters
#
#     Returns
#     -------
#     closeindices : list
#         ids of existing clusters that are close to the target cluster
#     '''
#     closeornot = [cal_distance(loc_tgt,loc) < MAX_THRESH_CENT_KM for loc in locs]
#     closeindices = [i for i, x in enumerate(closeornot) if x == True]
#     return closeindices
screbec's avatar
screbec committed

screbec's avatar
screbec committed
def cal_centroid(data):
Yang's avatar
Yang committed
    """ Calculate the centroid of a list of points
screbec's avatar
screbec committed
    Parameters
    ----------
    data : list of [lat,lon]
        point location

    Returns
    -------
    xct : float
        lat of centroid
    yct : float
        lon of centroid
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    x, y = zip(*(data))
    l = len(x)
Yang's avatar
Yang committed
    xct, yct = sum(x) / l, sum(y) / l
screbec's avatar
screbec committed

Yang's avatar
Yang committed
    return xct, yct