Skip to content

river module

CalculateAngle(clCleaned)

calculate the orthogonal direction of each pixel of the centerline

Source code in geemap/algorithms/river.py
def CalculateAngle(clCleaned):
    """calculate the orthogonal direction of each pixel of the centerline"""

    w3 = ee.Kernel.fixed(
        9,
        9,
        [
            [135.0, 126.9, 116.6, 104.0, 90.0, 76.0, 63.4, 53.1, 45.0],
            [143.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 36.9],
            [153.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 26.6],
            [166.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 14.0],
            [180.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1e-5],
            [194.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 346.0],
            [206.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 333.4],
            [216.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 323.1],
            [225.0, 233.1, 243.4, 256.0, 270.0, 284.0, 296.6, 306.9, 315.0],
        ],
    )

    combinedReducer = ee.Reducer.sum().combine(ee.Reducer.count(), None, True)

    clAngle = (
        clCleaned.mask(clCleaned)
        .rename(["clCleaned"])
        .reduceNeighborhood(
            reducer=combinedReducer, kernel=w3, inputWeight="kernel", skipMasked=True
        )
    )

    ## mask calculating when there are more than two inputs into the angle calculation
    clAngleNorm = (
        clAngle.select("clCleaned_sum")
        .divide(clAngle.select("clCleaned_count"))
        .mask(clAngle.select("clCleaned_count").gt(2).Not())
    )

    ## if only one input into the angle calculation, rotate it by 90 degrees to get the orthogonal
    clAngleNorm = clAngleNorm.where(
        clAngle.select("clCleaned_count").eq(1), clAngleNorm.add(ee.Image(90))
    )

    return clAngleNorm.rename(["orthDegree"])

CleanCenterline(cl1px, maxBranchLengthToRemove, rmCorners)

clean the 1px centerline: 1. remove branches 2. remove corners to insure 1px width (optional)

Source code in geemap/algorithms/river.py
def CleanCenterline(cl1px, maxBranchLengthToRemove, rmCorners):
    """clean the 1px centerline:
    1. remove branches
    2. remove corners to insure 1px width (optional)
    """

    ## find the number of connecting pixels (8-connectivity)
    nearbyPoints = cl1px.mask(cl1px).reduceNeighborhood(
        reducer=ee.Reducer.count(), kernel=ee.Kernel.circle(1.5), skipMasked=True
    )

    ## define ends
    endsByNeighbors = nearbyPoints.lte(2)

    ## define joint points
    joints = nearbyPoints.gte(4)

    costMap = (
        cl1px.mask(cl1px)
        .updateMask(joints.Not())
        .cumulativeCost(
            source=endsByNeighbors.mask(endsByNeighbors),
            maxDistance=maxBranchLengthToRemove,
            geodeticDistance=True,
        )
    )

    branchMask = costMap.gte(0).unmask(0)
    cl1Cleaned = cl1px.updateMask(branchMask.Not())  # mask short branches;
    ends = ExtractEndpoints(cl1Cleaned)
    cl1Cleaned = cl1Cleaned.updateMask(ends.Not())

    if rmCorners:
        corners = ExtractCorners(cl1Cleaned)
        cl1Cleaned = cl1Cleaned.updateMask(corners.Not())

    return cl1Cleaned

ExtractCorners(CL1px)

calculate corners in the one pixel centerline

Source code in geemap/algorithms/river.py
def ExtractCorners(CL1px):
    """calculate corners in the one pixel centerline"""

    se1w = [[2, 2, 0], [2, 1, 1], [0, 1, 0]]

    se11 = ee.Kernel.fixed(3, 3, splitKernel(se1w, 1))
    se12 = ee.Kernel.fixed(3, 3, splitKernel(se1w, 2))

    result = CL1px
    # // the for loop removes the identified corners from the input image

    i = 0
    while i < 4:  # rotate kernels

        result = result.subtract(hitOrMiss(result, se11, se12))

        se11 = se11.rotate(1)
        se12 = se12.rotate(1)

        i = i + 1

    cornerPoints = CL1px.subtract(result)
    return cornerPoints

ExtractEndpoints(CL1px)

calculate end points in the one pixel centerline

Source code in geemap/algorithms/river.py
def ExtractEndpoints(CL1px):
    """calculate end points in the one pixel centerline"""

    se1w = [[0, 0, 0], [2, 1, 2], [2, 2, 2]]

    se11 = ee.Kernel.fixed(3, 3, splitKernel(se1w, 1))
    se12 = ee.Kernel.fixed(3, 3, splitKernel(se1w, 2))

    result = CL1px

    # // the for loop removes the identified endpoints from the input image
    i = 0
    while i < 4:  # rotate kernels
        result = result.subtract(hitOrMiss(CL1px, se11, se12))
        se11 = se11.rotate(1)
        se12 = se12.rotate(1)
        i = i + 1
    endpoints = CL1px.subtract(result)
    return endpoints

GetWidth(clAngleNorm, segmentInfo, endInfo, DM, crs, bound, scale, sceneID, note)

calculate the width of the river at each centerline pixel, measured according to the orthgonal direction of the river

Source code in geemap/algorithms/river.py
def GetWidth(clAngleNorm, segmentInfo, endInfo, DM, crs, bound, scale, sceneID, note):
    """calculate the width of the river at each centerline pixel, measured according to the orthgonal direction of the river"""

    def GetXsectionEnds(f):
        xc = ee.Number(f.get("x"))
        yc = ee.Number(f.get("y"))
        orthRad = ee.Number(f.get("angle")).divide(180).multiply(math.pi)

        width = ee.Number(f.get("toBankDistance")).multiply(1.5)
        cosRad = width.multiply(orthRad.cos())
        sinRad = width.multiply(orthRad.sin())
        p1 = ee.Geometry.Point([xc.add(cosRad), yc.add(sinRad)], crs)
        p2 = ee.Geometry.Point([xc.subtract(cosRad), yc.subtract(sinRad)], crs)

        xlEnds = ee.Feature(
            ee.Geometry.MultiPoint([p1, p2]).buffer(30),
            {
                "xc": xc,
                "yc": yc,
                "longitude": f.get("lon"),
                "latitude": f.get("lat"),
                "orthogonalDirection": orthRad,
                "MLength": width.multiply(2),
                "p1": p1,
                "p2": p2,
                "crs": crs,
                "image_id": sceneID,
                "note": note,
            },
        )

        return xlEnds

    def SwitchGeometry(f):
        return (
            f.setGeometry(
                ee.Geometry.LineString(
                    coords=[f.get("p1"), f.get("p2")], proj=crs, geodesic=False
                )
            )
            .set("p1", None)
            .set("p2", None)
        )  # remove p1 and p2

    ## convert centerline image to a list. prepare for map function
    clPoints = (
        clAngleNorm.rename(["angle"])
        .addBands(ee.Image.pixelCoordinates(crs))
        .addBands(ee.Image.pixelLonLat().rename(["lon", "lat"]))
        .addBands(DM.rename(["toBankDistance"]))
        .sample(region=bound, scale=scale, projection=None, factor=1, dropNulls=True)
    )

    ## calculate the cross-section lines, returning a featureCollection
    xsectionsEnds = clPoints.map(GetXsectionEnds)

    ## calculate the flags at the xsection line end points
    endStat = endInfo.reduceRegions(
        collection=xsectionsEnds,
        reducer=ee.Reducer.anyNonZero().combine(
            ee.Reducer.count(), None, True
        ),  # test endpoints type
        scale=scale,
        crs=crs,
    )

    ## calculate the width of the river and other flags along the xsection lines
    xsections1 = endStat.map(SwitchGeometry)
    combinedReducer = ee.Reducer.mean()
    xsections = segmentInfo.reduceRegions(
        collection=xsections1, reducer=combinedReducer, scale=scale, crs=crs
    )

    return xsections

Skeletonize(image, iterations, method)

perform skeletonization

Source code in geemap/algorithms/river.py
def Skeletonize(image, iterations, method):
    """perform skeletonization"""

    se1w = [[2, 2, 2], [0, 1, 0], [1, 1, 1]]

    if method == 2:
        se1w = [[2, 2, 2], [0, 1, 0], [0, 1, 0]]

    se11 = ee.Kernel.fixed(3, 3, splitKernel(se1w, 1))
    se12 = ee.Kernel.fixed(3, 3, splitKernel(se1w, 2))

    se2w = [[2, 2, 0], [2, 1, 1], [0, 1, 0]]

    if method == 2:
        se2w = [[2, 2, 0], [2, 1, 1], [0, 1, 1]]

    se21 = ee.Kernel.fixed(3, 3, splitKernel(se2w, 1))
    se22 = ee.Kernel.fixed(3, 3, splitKernel(se2w, 2))

    result = image

    i = 0
    while i < iterations:
        j = 0
        while j < 4:
            result = result.subtract(hitOrMiss(result, se11, se12))
            se11 = se11.rotate(1)
            se12 = se12.rotate(1)
            result = result.subtract(hitOrMiss(result, se21, se22))
            se21 = se21.rotate(1)
            se22 = se22.rotate(1)
            j = j + 1
        i = i + 1

    return result.rename(["clRaw"])

hitOrMiss(image, se1, se2)

perform hitOrMiss transform (adapted from [citation])

Source code in geemap/algorithms/river.py
def hitOrMiss(image, se1, se2):
    """perform hitOrMiss transform (adapted from [citation])"""
    e1 = image.reduceNeighborhood(ee.Reducer.min(), se1)
    e2 = image.Not().reduceNeighborhood(ee.Reducer.min(), se2)
    return e1.And(e2)

maximum_no_of_tasks(MaxNActive, waitingPeriod)

maintain a maximum number of active tasks

Source code in geemap/algorithms/river.py
def maximum_no_of_tasks(MaxNActive, waitingPeriod):
    """maintain a maximum number of active tasks"""
    import time
    import ee

    ee.Initialize()

    time.sleep(10)
    ## initialize submitting jobs
    ts = list(ee.batch.Task.list())

    NActive = 0
    for task in ts:
        if "RUNNING" in str(task) or "READY" in str(task):
            NActive += 1
    ## wait if the number of current active tasks reach the maximum number
    ## defined in MaxNActive
    while NActive >= MaxNActive:
        time.sleep(
            waitingPeriod
        )  # if reach or over maximum no. of active tasks, wait for 2min and check again
        ts = list(ee.batch.Task.list())
        NActive = 0
        for task in ts:
            if "RUNNING" in str(task) or "READY" in str(task):
                NActive += 1
    return ()

merge_collections_std_bandnames_collection1tier1_sr()

merge landsat 5, 7, 8 collection 1 tier 1 SR imageCollections and standardize band names

Source code in geemap/algorithms/river.py
def merge_collections_std_bandnames_collection1tier1_sr():
    """merge landsat 5, 7, 8 collection 1 tier 1 SR imageCollections and standardize band names"""
    ## standardize band names
    bn8 = ["B1", "B2", "B3", "B4", "B6", "pixel_qa", "B5", "B7"]
    bn7 = ["B1", "B1", "B2", "B3", "B5", "pixel_qa", "B4", "B7"]
    bn5 = ["B1", "B1", "B2", "B3", "B5", "pixel_qa", "B4", "B7"]
    bns = ["uBlue", "Blue", "Green", "Red", "Swir1", "BQA", "Nir", "Swir2"]

    # create a merged collection from landsat 5, 7, and 8
    ls5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR").select(bn5, bns)

    ls7 = (
        ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")
        .filterDate("1999-04-15", "2003-05-30")
        .select(bn7, bns)
    )

    ls8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").select(bn8, bns)

    merged = ls5.merge(ls7).merge(ls8)

    return merged

rwc(image, description=None, folder='', file_format='shp', aoi=None, water_method='Jones2019', max_dist=4000, fill_size=333, max_dist_branch_removal=500, return_fc=False)

Calculate river centerlines and widths for one Landsat SR image.

Parameters:

Name Type Description Default
image str | ee.Image

LANDSAT_ID for any Landsat 5, 7, and 8 SR scene. For example, LC08_L1TP_022034_20130422_20170310_01_T1.

required
description str

File name of the output file. Defaults to None.

None
folder str

Folder name within Google Drive to save the exported file. Defaults to "", which is the root directory.

''
file_format str

The supported file format include shp, csv, json, kml, kmz, and TFRecord. Defaults to 'shp'. Defaults to "shp".

'shp'
aoi ee.Geometry.Polygon

A polygon (or rectangle) geometry define the area of interest. Only widths and centerline from this area will be calculated. Defaults to None.

None
water_method str

Water classification method ('Jones2019' or 'Zou2018'). Defaults to "Jones2019".

'Jones2019'
max_dist int

Maximum distance (unit: meters) to check water pixel's connectivity to GRWL centerline. Defaults to 4000.

4000
fill_size int

islands or bars smaller than this value (unit: pixels) will be removed before calculating centerline. Defaults to 333.

333
max_dist_branch_removal int

length of pruning. Spurious branch of the initial centerline will be removed by this length (unit: pixels). Defaults to 500.

500
return_fc(bool, optional

whether to return the result as an ee.FeatureColleciton. Defaults to False.

required
Source code in geemap/algorithms/river.py
def rwc(
    image,
    description=None,
    folder="",
    file_format="shp",
    aoi=None,
    water_method="Jones2019",
    max_dist=4000,
    fill_size=333,
    max_dist_branch_removal=500,
    return_fc=False,
):
    """Calculate river centerlines and widths for one Landsat SR image.

    Args:
        image (str | ee.Image): LANDSAT_ID for any Landsat 5, 7, and 8 SR scene. For example, LC08_L1TP_022034_20130422_20170310_01_T1.
        description (str, optional): File name of the output file. Defaults to None.
        folder (str, optional): Folder name within Google Drive to save the exported file. Defaults to "", which is the root directory.
        file_format (str, optional): The supported file format include shp, csv, json, kml, kmz, and TFRecord. Defaults to 'shp'. Defaults to "shp".
        aoi (ee.Geometry.Polygon, optional): A polygon (or rectangle) geometry define the area of interest. Only widths and centerline from this area will be calculated. Defaults to None.
        water_method (str, optional): Water classification method ('Jones2019' or 'Zou2018'). Defaults to "Jones2019".
        max_dist (int, optional): Maximum distance (unit: meters) to check water pixel's connectivity to GRWL centerline. Defaults to 4000.
        fill_size (int, optional): islands or bars smaller than this value (unit: pixels) will be removed before calculating centerline. Defaults to 333.
        max_dist_branch_removal (int, optional): length of pruning. Spurious branch of the initial centerline will be removed by this length (unit: pixels). Defaults to 500.
        return_fc(bool, optional): whether to return the result as an ee.FeatureColleciton. Defaults to False.
    """

    img = str_to_ee(image)
    if description is None:
        if isinstance(image, str):
            description = image
        else:
            description = img.get("LANDSAT_ID").getInfo()

    gen = rwGenSR(
        aoi=aoi,
        WATER_METHOD=water_method,
        MAXDISTANCE=max_dist,
        FILL_SIZE=fill_size,
        MAXDISTANCE_BRANCH_REMOVAL=max_dist_branch_removal,
    )

    width_fc = gen(img)

    if return_fc:
        return width_fc
    else:
        taskWidth = ee.batch.Export.table.toDrive(
            collection=width_fc,
            description=description,
            folder=folder,
            fileFormat=file_format,
        )
        taskWidth.start()
        print(description, "will be exported to", folder, "as", file_format, "file")

rwc_batch(images, folder='', file_format='shp', aoi=None, water_method='Jones2019', max_dist=4000, fill_size=333, max_dist_branch_remove=500)

Calculate river centerlines and widths for multiple Landsat SR images.

Parameters:

Name Type Description Default
images str | list | ee.ImageCollection

An input csv file containing a list of Landsat IDs (e.g., LC08_L1TP_022034_20130422_20170310_01_T1)

required
folder str

Folder name within Google Drive to save the exported file. Defaults to "", which is the root directory.

''
file_format str

The supported file format include shp, csv, json, kml, kmz, and TFRecord. Defaults to 'shp'. Defaults to "shp".

'shp'
aoi ee.Geometry.Polygon

A polygon (or rectangle) geometry define the area of interest. Only widths and centerline from this area will be calculated. Defaults to None.

None
water_method str

Water classification method ('Jones2019' or 'Zou2018'). Defaults to "Jones2019".

'Jones2019'
max_dist int

Maximum distance (unit: meters) to check water pixel's connectivity to GRWL centerline. Defaults to 4000.

4000
fill_size int

islands or bars smaller than this value (unit: pixels) will be removed before calculating centerline. Defaults to 333.

333
max_dist_branch_removal int

length of pruning. Spurious branch of the initial centerline will be removed by this length (unit: pixels). Defaults to 500.

required
Source code in geemap/algorithms/river.py
def rwc_batch(
    images,
    folder="",
    file_format="shp",
    aoi=None,
    water_method="Jones2019",
    max_dist=4000,
    fill_size=333,
    max_dist_branch_remove=500,
):
    """Calculate river centerlines and widths for multiple Landsat SR images.

    Args:
        images (str | list | ee.ImageCollection): An input csv file containing a list of Landsat IDs (e.g., LC08_L1TP_022034_20130422_20170310_01_T1)
        folder (str, optional): Folder name within Google Drive to save the exported file. Defaults to "", which is the root directory.
        file_format (str, optional): The supported file format include shp, csv, json, kml, kmz, and TFRecord. Defaults to 'shp'. Defaults to "shp".
        aoi (ee.Geometry.Polygon, optional): A polygon (or rectangle) geometry define the area of interest. Only widths and centerline from this area will be calculated. Defaults to None.
        water_method (str, optional): Water classification method ('Jones2019' or 'Zou2018'). Defaults to "Jones2019".
        max_dist (int, optional): Maximum distance (unit: meters) to check water pixel's connectivity to GRWL centerline. Defaults to 4000.
        fill_size (int, optional): islands or bars smaller than this value (unit: pixels) will be removed before calculating centerline. Defaults to 333.
        max_dist_branch_removal (int, optional): length of pruning. Spurious branch of the initial centerline will be removed by this length (unit: pixels). Defaults to 500.

    """
    import pandas as pd

    if isinstance(images, str):

        imageInfo = pd.read_csv(
            images, dtype={"Point_ID": np.unicode_, "LANDSAT_ID": np.unicode_}
        )
        sceneIDList = imageInfo["LANDSAT_ID"].values.tolist()
    elif isinstance(images, ee.ImageCollection):
        sceneIDList = images.aggregate_array("LANDSAT_ID").getInfo()
    elif isinstance(images, list):
        sceneIDList = images
    else:
        raise Exception("images must be a list of Landsat IDs or an ee.ImageCollection")

    for scene in sceneIDList:
        rwc(
            scene,
            scene,
            folder,
            file_format,
            aoi,
            water_method,
            max_dist,
            fill_size,
            max_dist_branch_remove,
        )

splitKernel(kernel, value)

recalculate the kernel according to the given foreground value

Source code in geemap/algorithms/river.py
def splitKernel(kernel, value):
    """recalculate the kernel according to the given foreground value"""
    kernel = np.array(kernel)
    result = kernel
    r = 0
    while r < kernel.shape[0]:
        c = 0
        while c < kernel.shape[1]:
            if kernel[r][c] == value:
                result[r][c] = 1
            else:
                result[r][c] = 0
            c = c + 1
        r = r + 1
    return result.tolist()

str_to_ee(id)

Convert an image id to ee.Image.

Parameters:

Name Type Description Default
id str | ee.Image

A string representing an ee.Image, such as LC08_022034_20180303, LANDSAT/LC08/C01/T1_SR/LC08_022034_20180303, LC08_L1TP_022034_20180303_20180319_01_T1

required

Exceptions:

Type Description
Exception

The image id is not recognized. It must be retrieved using either LANDSAT_ID, system:index, or system:id

Exception

The image id is not recognized. It must be either ee.Image or a string retrieved using either LANDSAT_ID, system:index, or system:id

Returns:

Type Description
ee.Image

An ee.Image.

Source code in geemap/algorithms/river.py
def str_to_ee(id):
    """Convert an image id to ee.Image.

    Args:
        id (str | ee.Image): A string representing an ee.Image, such as LC08_022034_20180303, LANDSAT/LC08/C01/T1_SR/LC08_022034_20180303, LC08_L1TP_022034_20180303_20180319_01_T1

    Raises:
        Exception: The image id is not recognized. It must be retrieved using either LANDSAT_ID, system:index, or system:id
        Exception: The image id is not recognized. It must be either ee.Image or a string retrieved using either LANDSAT_ID, system:index, or system:id

    Returns:
        ee.Image: An ee.Image.
    """
    if isinstance(id, str):
        if len(id) == 43 and "/" in id:
            id = ee.Image(id).get("LANDSAT_ID").getInfo()
            return ee.Image(
                merge_collections_std_bandnames_collection1tier1_sr()
                .filterMetadata("LANDSAT_ID", "equals", id)
                .first()
            )

        elif len(id) == 40:
            return ee.Image(
                merge_collections_std_bandnames_collection1tier1_sr()
                .filterMetadata("LANDSAT_ID", "equals", id)
                .first()
            )
        else:
            raise Exception(
                "The image id is not recognized. It must be retrieved using either LANDSAT_ID, system:index"
            )

    elif isinstance(id, ee.Image):
        id = id.get("LANDSAT_ID").getInfo()
        return ee.Image(
            merge_collections_std_bandnames_collection1tier1_sr()
            .filterMetadata("LANDSAT_ID", "equals", id)
            .first()
        )
    else:
        raise Exception(
            "The image id is not recognized. It must be either ee.Image or a string retrieved using either LANDSAT_ID, system:index, or system:id"
        )

Last update: 2021-10-26