Skip to content
Snippets Groups Projects
FireVector.py 9.05 KiB
Newer Older
screbec's avatar
screbec committed
""" FireVector
This is the module used for vector related calculations
"""

Yang's avatar
Yang committed

screbec's avatar
screbec committed
def alpha_shape(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    source: http://blog.thehumangeo.com/2014/05/12/drawing-boundaries-in-python/
    """
    from shapely.ops import cascaded_union, polygonize
    from scipy.spatial import Delaunay
    import numpy as np
    import math
    import shapely.geometry as geometry

    if len(points) < 4:
        # When you have a triangle, there is no sense in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull
screbec's avatar
screbec committed
    def add_edge(edges, edge_points, coords, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
Yang's avatar
Yang committed
        edges.add((i, j))
        edge_points.append(coords[[i, j]])

screbec's avatar
screbec committed
    coords = np.array([point for point in points])
    tri = Delaunay(coords)
    edges = set()
    edge_points = []

    # loop over triangles:
    # ia, ib, ic = indices of corner points of the
    # triangle
    # for ia, ib, ic in tri.vertices:
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]
        # Lengths of sides of triangle
Yang's avatar
Yang committed
        a = math.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = math.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = math.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
screbec's avatar
screbec committed
        # Semiperimeter of triangle
Yang's avatar
Yang committed
        s = (a + b + c) / 2.0
screbec's avatar
screbec committed

        # Area of triangle by Heron's formula
        try:
Yang's avatar
Yang committed
            area = math.sqrt(s * (s - a) * (s - b) * (s - c))
screbec's avatar
screbec committed
        except:
            area = 0

        if area > 0:  # modified by YC
Yang's avatar
Yang committed
            circum_r = a * b * c / (4.0 * area)
screbec's avatar
screbec committed
        else:
            circum_r = 0
        # Here's the radius filter.
Yang's avatar
Yang committed
        # print circum_r
Yang's avatar
Yang committed
        # if circum_r < 1.0/alpha:
        if circum_r < alpha:
screbec's avatar
screbec committed
            add_edge(edges, edge_points, coords, ia, ib)
            add_edge(edges, edge_points, coords, ib, ic)
            add_edge(edges, edge_points, coords, ic, ia)
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))

    # return triangles, edge_points
    return cascaded_union(triangles), edge_points

Yang's avatar
Yang committed

def addbuffer(geom_ll, vbuf):
    """ add a geometry in geographic projection with a buffer in meter
screbec's avatar
screbec committed
    Parameters
    ----------
    geom_ll : shapely geometry
        the geometry (in lat /lon)
    vbuf : float
        the buffer value (in meter)
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    # to speed up the calculation, use the following approach
    # the buffer in geographic projection is calculated using the centroid lat value
Yang's avatar
Yang committed
    # from FireConsts import EARTH_RADIUS_KM
    # import numpy as np
    # lat = geom_ll.centroid.y
    # ldeg = (EARTH_RADIUS_KM*np.cos(np.deg2rad(lat))*1000*2*np.pi/360)
    # vbufdeg = vbuf/ldeg    # convert vbuf in m to degs
    # geom_ll_buf = geom_ll.buffer(vbufdeg)
screbec's avatar
screbec committed

Yang's avatar
Yang committed
    return geom_ll.buffer(vbuf)
screbec's avatar
screbec committed

screbec's avatar
screbec committed
def doMultP(locs, buf):
Yang's avatar
Yang committed
    """ deirvve a MultipPolygon (bufferred MultiPoint) shape from given fire locations
screbec's avatar
screbec committed
    Parameters
    ----------
    locs : list (nx2)
        latitude and longitude values of all fire pixels
    Returns
    -------
    multP : MultiPolygon object
        calculated shape
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    from shapely.geometry import MultiPoint

    # MultiPoint shape
Yang's avatar
Yang committed
    geom = MultiPoint([(x, y) for x, y in locs])
screbec's avatar
screbec committed

    # Add buffer to form MultiPolygon shape
Yang's avatar
Yang committed
    geom = addbuffer(geom, buf)
screbec's avatar
screbec committed

Yang's avatar
Yang committed
    return geom
screbec's avatar
screbec committed

screbec's avatar
screbec committed
def doConvH(locs, buf):
Yang's avatar
Yang committed
    """ derive the convex hull given fire locations
screbec's avatar
screbec committed
    Parameters
    ----------
    locs : list (nx2)
        latitude and longitude values of all fire pixels
    Returns
    -------
    hull : Polygon object
        calculated hull shape
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed

    from scipy.spatial import ConvexHull
    from shapely.geometry import Polygon

    # calculate the convex hull using scipy.spatial.ConvexHull
    try:
        qhull = ConvexHull(locs)
    except:
        return None

    # derive qhull object vertices
Yang's avatar
Yang committed
    verts = [(locs[i][0], locs[i][1]) for i in qhull.vertices]
screbec's avatar
screbec committed

    # convert vertices to polygon
Yang's avatar
Yang committed
    hull = Polygon(verts)
screbec's avatar
screbec committed

    # Add buffer
    # hull = hull.buffer(VIIRSbuf)
screbec's avatar
screbec committed
    hull = addbuffer(hull, buf)
screbec's avatar
screbec committed
    return hull

Yang's avatar
Yang committed

def doConcH(locs, buf, alpha=100):
    """ derive the concave hull given fire locations
screbec's avatar
screbec committed
    Parameters
    ----------
    locs : list (nx2)
        latitude and longitude values of all fire pixels
    alpha : int
        alpha value (1/deg)
    Returns
    -------
    hull : Polygon or MultiPolygon object
        calculated hull shape
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed

    # calcluate alpha shape
    try:
        # switch lats and lons to use the alpha_shape
Yang's avatar
Yang committed
        # locs_lonlat = [(v[1],v[0]) for v in locs]
        # concave_hull, edge_points = alpha_shape(locs_lonlat,alpha=alpha)
screbec's avatar
screbec committed

Yang's avatar
Yang committed
        concave_hull, edge_points = alpha_shape(locs, alpha=alpha)
        if (
            concave_hull.area == 0
        ):  # sometimes the concave_hull algorithm returns a empty polygon
            concave_hull = doConvH(locs, buf)
screbec's avatar
screbec committed
            return concave_hull
        else:
            # Add buffer
            # hull = concave_hull.buffer(VIIRSbuf)
Yang's avatar
Yang committed
            hull = addbuffer(concave_hull, buf)
screbec's avatar
screbec committed
            return hull

    except:
        return None


Yang's avatar
Yang committed
def cal_hull(fp_locs, sensor="viirs"):
    """ wrapper to calculate the hull given fire locations.
screbec's avatar
screbec committed
        the returned hull type depends on the pixel number
    Parameters
    ----------
    fp_locs : list (nx2)
        latitude and longitude values of all fire pixels
    Returns
    -------
    hull : object
        calculated hull (a buffer of VIIRS half pixel size included)
Yang's avatar
Yang committed
    """
    from FireConsts import valpha, VIIRSbuf  # ,MCD64buf
screbec's avatar
screbec committed
    import numpy as np
Yang's avatar
Yang committed

screbec's avatar
screbec committed
    # set buffer according to sensor
Yang's avatar
Yang committed
    if sensor == "viirs":
screbec's avatar
screbec committed
        buf = VIIRSbuf
Yang's avatar
Yang committed
    elif sensor == "mcd64":
Yang's avatar
Yang committed
        buf = MCD64buf
screbec's avatar
screbec committed
    else:
        print('please set sensor to viirs or mcd64')
Yang's avatar
Yang committed

screbec's avatar
screbec committed
    # number of points
    nfp = len(fp_locs)

    # For cluster with 1-2 pixel, calculate hull using buffered points (MultiPolygon)
    if nfp < 3:
Yang's avatar
Yang committed
        hull = doMultP(fp_locs, buf)
screbec's avatar
screbec committed
    # For cluster with 3 pixels, calculate hull using convex hull
Yang's avatar
Yang committed
    elif nfp == 3:  # call doConvH to get the hull
        hull = doConvH(fp_locs, buf)
screbec's avatar
screbec committed
        if hull == None:  # in case where convex hull can't be determined, call doMultP
Yang's avatar
Yang committed
            hull = doMultP(fp_locs, buf)
        elif (
            hull.area == 0
        ):  # sometimes the fire pixels are in a traight line, in this case, also call doMultP
            hull = doMultP(fp_locs, buf)
screbec's avatar
screbec committed
    # For cluster with more than 3 pixels, calculate hull using alpha shape
Yang's avatar
Yang committed
    else:  # call doConcH to get the hull
screbec's avatar
screbec committed
        # derive the alpha value used in doConcH (in 1/deg)
Yang's avatar
Yang committed
        # x,y = zip(*fp_locs)
        # vdeg = sum(x)/len(x)
        # km1deg = 6371*np.cos(np.deg2rad(vdeg))*2*np.pi/360
        # valphadeg = 1/(valpha/1000/km1deg)   # in 1/deg
Yang's avatar
Yang committed
        hull = doConcH(fp_locs, buf, alpha=valpha)
        if hull == None:  # in case where convex hull can't be determined, call doMultP
            hull = doMultP(fp_locs, buf)
screbec's avatar
screbec committed
        elif hull.area == 0:
Yang's avatar
Yang committed
            hull = doMultP(fp_locs, buf)
screbec's avatar
screbec committed
    return hull

Yang's avatar
Yang committed

def cal_extpixels(fps, hull, alpha=100):
    """ calculate the exterior pixels around a hull
screbec's avatar
screbec committed
    Parameters
    ----------
    fps : list (nx3)
        list of all fire pixels
    hull : geometry, 'Polygon' | 'MultiPoint'
        the existing hull for the fire
    alpha : int
        alpha value (1/deg) to define the inward buffer
    Returns
    -------
    fps_ext : list (nx2)
        the exterior fire pixels for the fire
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    from shapely.geometry import Point
    from FireConsts import extbuffer

    # use alpha to define an inward buffer and an interior part of the hull
    # hts_buf = hull.buffer(-1/alpha)
    hts_buf = addbuffer(hull, -extbuffer)

    # loop over all pixel locations in fps to determine exterior pixels
    fps_ext = []
    for fp in fps:
        # point locations of each pixel (fp.oc[0]:lat, fp.loc[1]:lon)
Yang's avatar
Yang committed
        pt = Point(fp.loc[0], fp.loc[1])
screbec's avatar
screbec committed

        # exclude points in interior part of the hull
        if not hts_buf.contains(pt):
            fps_ext.append(fp)
    return fps_ext

screbec's avatar
screbec committed
def calConcHarea(hull):
Yang's avatar
Yang committed
    """ calculate area given the concave hull (km2)
screbec's avatar
screbec committed
    Parameters
    ----------
    hull : geometry, 'Polygon' | 'MultiPoint'
        the hull for the fire
    Returns
    -------
    farea : float
        the area (km2) of the polygon enclosed by the vertices
Yang's avatar
Yang committed
    """
screbec's avatar
screbec committed
    import math
    from FireConsts import EARTH_RADIUS_KM

Yang's avatar
Yang committed
    farea = hull.area  # area in deg**2
screbec's avatar
screbec committed

    if farea > 0:
        lat = hull.bounds[1]  # latitude in deg

        # convert deg**2 to km2
Yang's avatar
Yang committed
        farea = (
            farea
            * EARTH_RADIUS_KM ** 2
            * math.cos(math.radians(lat))
            * math.radians(1) ** 2
        )
screbec's avatar
screbec committed

    return farea