Skip to content

Description

Step purpose

Purpose: to apply the Otsu threshold to the normalized difference.

This step is highlighted below:

graph TB style E stroke:#f66,stroke-width:3px subgraph Process STAC item A[STAC Item] -.-> B A[STAC Item] -.-> C A[STAC Item] -.-> F subgraph scatter on bands B["crop(green)"]; C["crop(nir)"]; end B["crop(green)"] -.-> D[Normalized difference]; C["crop(nir)"] -.-> D[Normalized difference]; D == norm_diff.tif ==> E[Otsu threshold] end E == otsu.tif ==> F[Create STAC Catalog] F -.-> G[(storage)]

Code

The otsu Python script is a command-line tool for applying the Otsu threshold to a single input raster image.

It uses the click, rasterio, numpy, skimage.filters, and loguru libraries.

Here's an overview of what the script does:

  • It defines a command-line interface using the click library, with a single argument for providing the file path of the input raster image on which you want to apply the Otsu threshold.

  • The otsu function is the main entry point. It opens the input raster file specified as the argument.

  • It reads the data from the input raster using rasterio and also copies the metadata (e.g., projection, geotransform) of this raster to be used for the output.

  • It applies the Otsu threshold to the input array by calling the threshold function. The threshold_otsu function from skimage.filters is used to calculate the Otsu threshold. The thresholding process marks pixels as True or False based on whether they are greater than the calculated threshold.

  • It creates an output raster named "otsu.tif" using rasterio. This output raster will have the same metadata as the input raster.

  • It writes the binary image to the output raster using dst_dataset.write().

The result, a binary image where pixel values are either True or False based on the thresholding, will be saved as "otsu.tif" in the same directory where the script is executed.

The script is executable as a command-line tool as its usage is:

Usage: app.py [OPTIONS] RASTER

  Applies the Otsu threshold

Options:
  --help  Show this message and exit.

The Python code is provided here:

water-bodies/command-line-tools/otsu/app.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
"""Command line tool to apply the Otsu threshold to a raster"""
import click
import rasterio
import numpy as np
from skimage.filters import threshold_otsu
from loguru import logger


def threshold(data):
    """Returns the Otsu threshold of a numpy array"""
    return data > threshold_otsu(data[np.isfinite(data)])


@click.command(
    short_help="Otsu threshoold",
    help="Applies the Otsu threshold",
)
@click.argument("raster", nargs=1)
def otsu(raster):
    """Applies the Otsu threshold"""

    with rasterio.open(raster) as ds:
        array = ds.read(1)
        out_meta = ds.meta.copy()

    out_meta.update(
        {
            "dtype": "uint8",
            "driver": "COG",
            "tiled": True,
            "compress": "lzw",
            "blockxsize": 256,
            "blockysize": 256,
        }
    )

    logger.info(f"Applying the Otsu threshold to {raster}")

    with rasterio.open("otsu.tif", "w", **out_meta) as dst_dataset:
        logger.info(f"Write otsu.tif")
        dst_dataset.write(threshold(array), indexes=1)

    logger.info("Done!")


if __name__ == "__main__":
    otsu()