Add reclassify raster (#197)

Author: lehtonenp

docs/raster_processing/reclassify.md

@@ -0,0 +1,3 @@
+# Reclassify raster
+
+::: eis_toolkit.raster_processing.reclassify

eis_toolkit/raster_processing/reclassify.py

@@ -0,0 +1,413 @@
+from numbers import Number
+
+import mapclassify as mc
+import numpy as np
+import rasterio
+from beartype import beartype
+from beartype.typing import Optional, Sequence, Tuple
+
+from eis_toolkit.utilities.checks.raster import check_raster_bands
+
+
+def _reclassify_with_manual_breaks(  # type: ignore[no-any-unimported]
+    band: np.ndarray,
+    breaks: Sequence[int],
+) -> np.ndarray:
+
+    data = np.digitize(band, breaks)
+
+    return data
+
+
+@beartype
+def reclassify_with_manual_breaks(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader,
+    breaks: Sequence[int],
+    bands: Optional[Sequence[int]] = None,
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with manual breaks.
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        breaks: List of break values for the classification.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with manual breaks and metadata.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_manual_breaks(band_data, breaks)
+
+    return out_image, out_meta
+
+
+def _reclassify_with_defined_intervals(  # type: ignore[no-any-unimported]
+    band: np.ndarray,
+    interval_size: int,
+) -> np.ndarray:
+
+    _, edges = np.histogram(band, bins=interval_size)
+
+    data = np.digitize(band, edges)
+
+    return data
+
+
+@beartype
+def reclassify_with_defined_intervals(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader,
+    interval_size: int,
+    bands: Optional[Sequence[int]] = None,
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with defined intervals.
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        interval_size: The number of units in each interval.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with defined intervals and metadata.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_defined_intervals(band_data, interval_size)
+
+    return out_image, out_meta
+
+
+def _reclassify_with_equal_intervals(  # type: ignore[no-any-unimported]
+    band: np.ndarray,
+    number_of_intervals: int,
+) -> np.ndarray:
+
+    percentiles = np.linspace(0, 100, number_of_intervals)
+
+    intervals = np.percentile(band, percentiles)
+
+    data = np.digitize(band, intervals)
+
+    return data
+
+
+@beartype
+def reclassify_with_equal_intervals(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader,
+    number_of_intervals: int,
+    bands: Optional[Sequence[int]] = None,
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with equal intervals.
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        number_of_intervals: The number of intervals.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with equal intervals.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_equal_intervals(band_data, number_of_intervals)
+
+    return out_image, out_meta
+
+
+def _reclassify_with_quantiles(  # type: ignore[no-any-unimported]
+    band: np.ndarray,
+    number_of_quantiles: int,
+) -> np.ndarray:
+
+    intervals = [np.percentile(band, i * 100 / number_of_quantiles) for i in range(number_of_quantiles)]
+    data = np.digitize(band, intervals)
+
+    return data
+
+
+@beartype
+def reclassify_with_quantiles(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader,
+    number_of_quantiles: int,
+    bands: Optional[Sequence[int]] = None,
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with quantiles.
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        number_of_quantiles: The number of quantiles.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with quantiles and metadata.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_quantiles(band_data, number_of_quantiles)
+
+    return out_image, out_meta
+
+
+def _reclassify_with_natural_breaks(  # type: ignore[no-any-unimported]
+    band: np.ndarray,
+    number_of_classes: int,
+) -> np.ndarray:
+
+    breaks = mc.JenksCaspall(band, number_of_classes)
+    data = np.digitize(band, np.sort(breaks.bins))
+
+    return data
+
+
+@beartype
+def reclassify_with_natural_breaks(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader,
+    number_of_classes: int,
+    bands: Optional[Sequence[int]] = None,
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with natural breaks (Jenks Caspall).
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        number_of_classes: The number of classes.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with natural breaks (Jenks Caspall) and metadata.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_natural_breaks(band_data, number_of_classes)
+
+    return out_image, out_meta
+
+
+def _reclassify_with_geometrical_intervals(
+    band: np.ndarray, number_of_classes: int, nodata_value: Number
+) -> np.ndarray:
+
+    # nan_value is either a set integer (e.g. -9999) or np.nan
+    mask = band == nodata_value
+    masked_array = np.ma.masked_array(data=band, mask=mask)
+
+    median_value = np.ma.median(masked_array)
+    max_value = masked_array.max()
+    min_value = masked_array.min()
+
+    values_out = np.ma.zeros_like(masked_array)
+
+    # Determine the tail with larger length
+    if (median_value - min_value) < (max_value - median_value):  # Large end tail longer
+        tail_values = masked_array[np.ma.where((masked_array > median_value))]
+        range_tail = max_value - median_value
+        tail_values = tail_values - median_value + range_tail / 1000.0
+    else:  # Small end tail longer
+        tail_values = masked_array[np.ma.where((masked_array < median_value))]
+        range_tail = median_value - min_value
+        tail_values = tail_values - min_value + range_tail / 1000.0
+
+    min_tail = np.ma.min(tail_values)
+    max_tail = np.ma.max(tail_values)
+
+    # number of classes
+    factor = (max_tail / min_tail) ** (1 / number_of_classes)
+
+    interval_index = 1
+    break_points_tail = [min_tail]
+    break_points = [min_tail]
+    width = [0]
+
+    while break_points[-1] < max_tail:
+        interval_index += 1
+        break_points.append(min_tail * factor ** (interval_index - 1))
+        break_points_tail.append(break_points[-1])
+        width.append(break_points_tail[-1] - break_points_tail[0])
+    k = 0
+
+    for j in range(1, len(width) - 2):
+        values_out[
+            np.ma.where(((median_value + width[j]) < masked_array) & ((masked_array <= (median_value + width[j + 1]))))
+        ] = (j + 1)
+        values_out[
+            np.ma.where(((median_value - width[j]) > masked_array) & ((masked_array >= (median_value - width[j + 1]))))
+        ] = (-j - 1)
+        k = j
+
+    values_out[np.ma.where(((median_value + width[k + 1]) < masked_array))] = k + 1
+    values_out[np.ma.where(((median_value - width[k + 1]) > masked_array))] = -k - 1
+    values_out[np.ma.where(median_value == masked_array)] = 0
+
+    output = np.array(values_out)
+
+    return output
+
+
+@beartype
+def reclassify_with_geometrical_intervals(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader, number_of_classes: int, bands: Optional[Sequence[int]] = None
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with geometrical intervals.
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        number_of_classes: The number of classes. The true number of classes is at most double the amount,
+        depending how symmetrical the input data is.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with geometrical intervals and metadata.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+    nodata_value = raster.nodata
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_geometrical_intervals(band_data, number_of_classes, nodata_value)
+
+    return out_image, out_meta
+
+
+def _reclassify_with_standard_deviation(  # type: ignore[no-any-unimported]
+    band: np.ndarray,
+    number_of_intervals: int,
+) -> np.ndarray:
+
+    band_statistics = []
+
+    stddev = np.nanstd(band)
+    mean = np.nanmean(band)
+    band_statistics.append((mean, stddev))
+    interval_size = 2 * stddev / number_of_intervals
+
+    classified = np.empty_like(band)
+
+    below_mean = band < (mean - stddev)
+    above_mean = band > (mean + stddev)
+
+    classified[below_mean] = -number_of_intervals
+    classified[above_mean] = number_of_intervals
+
+    in_between = ~below_mean & ~above_mean
+    interval = ((band - (mean - stddev)) / interval_size).astype(int)
+    classified[in_between] = interval[in_between] - number_of_intervals // 2
+
+    return classified
+
+
+@beartype
+def reclassify_with_standard_deviation(  # type: ignore[no-any-unimported]
+    raster: rasterio.io.DatasetReader,
+    number_of_intervals: int,
+    bands: Optional[Sequence[int]] = None,
+) -> Tuple[np.ndarray, dict]:
+    """Classify raster with standard deviation.
+
+    If bands are not given, all bands are used for classification.
+
+    Args:
+        raster: Raster to be classified.
+        number_of_intervals: The number of intervals.
+        bands: Selected bands from multiband raster. Indexing begins from one. Defaults to None.
+
+    Returns:
+        Raster classified with standard deviation and metadata.
+
+    Raises:
+        InvalidParameterValueException: Bands contain negative values.
+    """
+
+    if bands is None or len(bands) == 0:
+        bands = range(1, raster.count + 1)
+    else:
+        check_raster_bands(raster, bands)
+
+    out_image = np.empty((len(bands), raster.height, raster.width))
+    out_meta = raster.meta.copy()
+
+    for i, band in enumerate(bands):
+        band_data = raster.read(band)
+        out_image[i] = _reclassify_with_standard_deviation(band_data, number_of_intervals)
+
+    return out_image, out_meta

environment.yml

@@ -21,6 +21,7 @@ dependencies:
   - rtree >= 1.0.1
   - typer >=0.9.0
   - imbalanced-learn >= 0.11.0
+  - mapclassify >= 2.6.1
   - esda >= 2.5.1
   # Dependencies for testing
   - pytest >=7.2.1

tests/raster_processing/reclassify_raster_test.py

@@ -0,0 +1,196 @@
+import numpy as np
+import rasterio
+from beartype.typing import Tuple
+
+from eis_toolkit.raster_processing import reclassify
+from tests.raster_processing.clip_test import raster_path as SMALL_RASTER_PATH
+
+TEST_ARRAY = np.array([[0, 10, 20, 30], [40, 50, 50, 60], [80, 80, 90, 90], [100, 100, 100, 100]])
+
+
+def test_reclassify_with_defined_intervals():
+    """Test raster with defined intervals."""
+    interval_size = 3
+
+    result = reclassify._reclassify_with_defined_intervals(TEST_ARRAY, interval_size)
+
+    expected_output = np.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_defined_intervals_main():
+    """Test raster with defined intervals parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_defined_intervals(
+            raster=raster,
+            interval_size=3,
+            bands=[1],
+        )
+
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)
+
+
+def test_reclassify_with_equal_intervals():
+    """Test raster with equal intervals."""
+    number_of_intervals = 10
+
+    result = reclassify._reclassify_with_equal_intervals(TEST_ARRAY, number_of_intervals)
+
+    expected_output = np.array([[1, 1, 2, 2], [3, 4, 4, 5], [6, 6, 7, 7], [10, 10, 10, 10]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_equal_intervals_main():
+    """Test raster with equal intervals parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_defined_intervals(
+            raster=raster,
+            interval_size=10,
+            bands=[1],
+        )
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)
+
+
+def test_reclassify_with_geometrical_intervals():
+    """Test raster with geometrical intervals."""
+    number_of_classes = 10
+    nodata_value = -9999
+
+    array_with_nan_value = np.array(
+        [[nodata_value, 10, 20, 30], [40, 50, 50, 60], [80, 80, 90, 90], [100, 100, 100, 100]]
+    )
+
+    result = reclassify._reclassify_with_geometrical_intervals(array_with_nan_value, number_of_classes, nodata_value)
+
+    expected_output = np.array([[0, -9, -9, -9], [-9, -9, -9, -9], [0, 0, 8, 8], [9, 9, 9, 9]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_geometrical_intervals_main():
+    """Test raster with geometrical intervals parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_geometrical_intervals(
+            raster=raster,
+            number_of_classes=10,
+            bands=[1],
+        )
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)
+
+
+def test_reclassify_with_manual_breaks():
+    """Test raster with manual break intervals."""
+    breaks = [20, 40, 60, 80]
+
+    result = reclassify._reclassify_with_manual_breaks(TEST_ARRAY, breaks)
+
+    expected_output = np.array([[0, 0, 1, 1], [2, 2, 2, 3], [4, 4, 4, 4], [4, 4, 4, 4]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_manual_breaks_main():
+    """Test raster with manual break intervals parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_manual_breaks(
+            raster=raster,
+            breaks=[2, 5, 9],
+            bands=[1],
+        )
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)
+
+
+def test_reclassify_with_natural_breaks():
+    """Test raster with natural breaks."""
+    number_of_classes = 10
+
+    result = reclassify._reclassify_with_natural_breaks(TEST_ARRAY, number_of_classes)
+
+    expected_output = np.array([[0, 1, 1, 2], [3, 4, 4, 5], [6, 6, 7, 7], [8, 8, 8, 8]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_natural_breaks_main():
+    """Test raster with natural break intervals parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_natural_breaks(
+            raster=raster,
+            number_of_classes=10,
+            bands=[1],
+        )
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)
+
+
+def test_reclassify_with_standard_deviation():
+    """Test raster with standard deviation intervals."""
+    number_of_intervals = 75
+
+    result = reclassify._reclassify_with_standard_deviation(TEST_ARRAY, number_of_intervals)
+
+    expected_output = np.array([[-75, -75, -75, -36], [-25, -14, -14, -3], [20, 20, 31, 31], [75, 75, 75, 75]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_standard_deviation_main():
+    """Test raster with standard_deviation intervals parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_standard_deviation(
+            raster=raster,
+            number_of_intervals=75,
+            bands=[1],
+        )
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)
+
+
+def test_reclassify_with_quantiles():
+    """Test raster with quantile intervals by."""
+    number_of_quantiles = 4
+
+    result = reclassify._reclassify_with_quantiles(TEST_ARRAY, number_of_quantiles)
+
+    expected_output = np.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
+
+    assert isinstance(result, np.ndarray)
+
+    np.testing.assert_allclose(result, expected_output)
+
+
+def test_reclassify_with_quantiles_main():
+    """Test raster with quantiles parameters."""
+    with rasterio.open(SMALL_RASTER_PATH) as raster:
+        result = reclassify.reclassify_with_quantiles(
+            raster=raster,
+            number_of_quantiles=4,
+            bands=[1],
+        )
+    assert isinstance(result, Tuple)
+    assert isinstance(result[0], np.ndarray)
+    assert isinstance(result[1], dict)