retrained models, support for l2a, refactoring

CHANGELOG.rst

@@ -1,6 +1,17 @@
 Changelog
 =========
 
+[1.0.0] (2024-11-XX)
+--------------------
+Added
+*******
+-support for L2A product level
+
+Changed
+*******
+- retrained all models with new architecture and training data
+- removed backwards compatibility with old band naming scheme
+
 [0.2.2] (2024-10-21)
 --------------------
 Added

MANIFEST.in

@@ -1,5 +1,11 @@
 include requirements.txt
 include README.md
 include LICENSE
-include ukis_csmask/model_4b.onnx
-include ukis_csmask/model_6b.onnx
+include ukis_csmask/model_4b_l1c.onnx
+include ukis_csmask/model_4b_l1c.json
+include ukis_csmask/model_4b_l2a.onnx
+include ukis_csmask/model_4b_l2a.json
+include ukis_csmask/model_6b_l1c.onnx
+include ukis_csmask/model_6b_l1c.json
+include ukis_csmask/model_6b_l2a.onnx
+include ukis_csmask/model_6b_l2a.json

README.md

@@ -8,7 +8,7 @@
 [![Code Style](https://img.shields.io/badge/code%20style-black-000000.svg)](https://black.readthedocs.io/en/stable/)
 [![DOI](https://zenodo.org/badge/328616234.svg)](https://zenodo.org/badge/latestdoi/328616234)
 
-UKIS Cloud Shadow MASK (ukis-csmask) package masks clouds and cloud shadows in Sentinel-2, Landsat-9, Landsat-8, Landsat-7 and Landsat-5 images. Masking is performed with a pre-trained convolution neural network. It is fast and works directly on Level-1C data (no atmospheric correction required). Images just need to be in Top Of Atmosphere (TOA) reflectance and include at least the "blue", "green", "red" and "nir" spectral bands. Best performance (in terms of accuracy and speed) is achieved when images also include "swir16" and "swir22" spectral bands and are resampled to approximately 30 m spatial resolution.
+UKIS Cloud Shadow MASK (ukis-csmask) package masks clouds and cloud shadows in Sentinel-2, Landsat-9, Landsat-8, Landsat-7 and Landsat-5 images. Masking is performed with a pre-trained convolution neural network. It is fast and works with both Level-1C (no atmospheric correction) and Level-2A (atmospherically corrected) data. Images just need to be in reflectance and include at least the "blue", "green", "red" and "nir" spectral bands. Best performance (in terms of accuracy and speed) is achieved when images also include "swir16" and "swir22" spectral bands and are resampled to approximately 30 m spatial resolution.
 
 This [publication](https://doi.org/10.1016/j.rse.2019.05.022) provides further insight into the underlying algorithm and compares it to the widely used [Fmask](http://www.pythonfmask.org/en/latest/) algorithm across a heterogeneous test dataset.
 
@@ -23,8 +23,8 @@ If you use ukis-csmask in your work, please consider citing one of the above pub
 
 ![Examples](img/examples.png)
 
-## Example (Sentinel 2)
-Here's an example on how to compute a cloud and cloud shadow mask from an image. Please note that here we use [ukis-pysat](https://github.com/dlr-eoc/ukis-pysat) for convencience image handling, but you can also work directly with [numpy](https://numpy.org/) arrays.
+## Example (Sentinel-2 Level-1C)
+Here's an example on how to compute a cloud and cloud shadow mask from an image. Please note that here we use [ukis-pysat](https://github.com/dlr-eoc/ukis-pysat) for convencience image handling, but you can also work directly with [numpy](https://numpy.org/) arrays. Further examples can be found [here](examples).
 
 ````python
 from ukis_csmask.mask import CSmask
@@ -64,49 +64,6 @@ csmask_csm.write_to_file("sentinel2_csm.tif", dtype="uint8", compress="PACKBITS"
 csmask_valid.write_to_file("sentinel2_valid.tif", dtype="uint8", compress="PACKBITS", kwargs={"nbits":2})
 ````
 
-## Example (Landsat 8)
-Here's a similar example based on Landsat 8.
-
-````python
-import rasterio
-import numpy as np
-from ukis_csmask.mask import CSmask
-from ukis_pysat.raster import Image, Platform
-
-# set Landsat 8 source path and prefix (example)
-data_path = "/your_data_path/"
-L8_file_prefix = "LC08_L1TP_191015_20210428_20210507_02_T1"
-
-data_path = data_path+L8_file_prefix+"/"
-mtl_file  = data_path+L8_file_prefix+"_MTL.txt"
-
-# stack [B2:'Blue', B3:'Green', B4:'Red', B5:'NIR', B6:'SWIR1', B7:'SWIR2'] as numpy array
-L8_band_files  = [data_path+L8_file_prefix+'_B'+ x + '.TIF' for x in [str(x+2) for x in range(6)]]
-
-# >> adopted from https://gis.stackexchange.com/questions/223910/using-rasterio-or-gdal-to-stack-multiple-bands-without-using-subprocess-commands
-# read metadata of first file
-with rasterio.open(L8_band_files[0]) as src0:
-    meta = src0.meta
-# update meta to reflect the number of layers
-meta.update(count = len(L8_band_files))
-# read each layer and append it to numpy array
-L8_bands = []
-for id, layer in enumerate(L8_band_files, start=1):
-    with rasterio.open(layer) as src1:
-        L8_bands.append(src1.read(1))
-L8_bands = np.stack(L8_bands,axis=2)
-# <<
-
-img = Image(data=L8_bands, crs = meta['crs'], transform = meta['transform'], dimorder="last")
-
-img.dn2toa(
-        platform=Platform.Landsat8,
-        mtl_file=mtl_file,
-        wavelengths = ["blue", "green", "red", "nir", "swir16", "swir22"]
-)
-# >> proceed by analogy with Sentinel 2 example
-````
-
 ## Installation
 The easiest way to install ukis-csmask is through pip. To install ukis-csmask with [default CPU provider](https://onnxruntime.ai/docs/execution-providers/) run the following.
 

examples/landsat_l2a.ipynb

@@ -0,0 +1,137 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Segment clouds and cloud shadows in Landsat images (L2A)\n",
+    "This notebook shows an example on how to use [ukis-csmask]() to segment clouds and cloud shadows in Level-2A images from Landsat-9, Landsat-8, Landsat-7 and Landsat-5 satellites. Images are acquired from [Planetary Computer]() and are prepared according to [ukis-csmask]() requirements.\n",
+    "\n",
+    "> NOTE: to run this notebook, we first need to install some additional dependencies for work with Planetary Computer\n",
+    "```shell\n",
+    "pip install planetary_computer rasterio rioxarray pystac-client odc-stac tqdm\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "703f3744-902d-470b-a80f-9a8d3ea08dfa",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import planetary_computer as pc\n",
+    "\n",
+    "from odc.stac import load\n",
+    "from pystac_client import Client\n",
+    "from tqdm import tqdm\n",
+    "from ukis_csmask.mask import CSmask"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8ca03c78-1e24-479c-9786-a1b43206a08b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stac_api_endpoint = \"https://planetarycomputer.microsoft.com/api/stac/v1\"\n",
+    "collections = [\"landsat-c2-l2\"]\n",
+    "# TODO: add ids for Landsat9, Landsat7, Landsat5\n",
+    "ids = [\"LC08_L2SP_221081_20240508_02_T1\"]\n",
+    "product_level = \"l2a\"\n",
+    "band_order = [\"blue\", \"green\", \"red\", \"nir\", \"swir16\", \"swir22\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b568942-84e6-4baf-b490-a213b3787f80",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# search catalog by scene id\n",
+    "catalog = Client.open(stac_api_endpoint)\n",
+    "search = catalog.search(collections=collections, ids=ids)\n",
+    "items = [item for item in search.items()]\n",
+    "items_cnt = len(items)\n",
+    "print(f\"Search returned {items_cnt} item(s)\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "68eb9c30-06f7-409e-914d-21e00f45de99",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for item in tqdm(items, total=items_cnt, desc=\"Predict images\"):\n",
+    "    # near infrared band has different alias in landsat collections\n",
+    "    bands = [b.replace(\"nir\", \"nir08\") for b in band_order]\n",
+    "\n",
+    "    # load and preprocess image\n",
+    "    arr = (\n",
+    "        load(\n",
+    "            items=[item],\n",
+    "            bands=bands,\n",
+    "            resolution=30,\n",
+    "            dtype=\"float32\",\n",
+    "            patch_url=pc.sign,\n",
+    "        )\n",
+    "        .to_dataarray()\n",
+    "        .squeeze()\n",
+    "        .drop_vars(\"time\")\n",
+    "    )\n",
+    "    arr = arr.rename({\"variable\": \"band\"})\n",
+    "\n",
+    "    # use band-specific rescale factors to convert DN to reflectance\n",
+    "    for idx, band_name in enumerate(bands):\n",
+    "        band_info = item.assets[band_name].extra_fields[\"raster:bands\"][0]\n",
+    "        arr[idx, :, :] = arr.sel(band=str(band_name)).astype(np.float32) * band_info[\"scale\"]\n",
+    "        arr[idx, :, :] += band_info[\"offset\"]\n",
+    "        arr[idx, :, :] = arr[idx, :, :].clip(min=0.0, max=1.0)\n",
+    "\n",
+    "    # compute cloud and cloud shadow mask\n",
+    "    csmask = CSmask(\n",
+    "        img=np.moveaxis(arr, 0, -1),\n",
+    "        band_order=band_order,\n",
+    "        product_level=product_level,\n",
+    "        nodata_value=0,\n",
+    "        invalid_buffer=4,\n",
+    "        intra_op_num_threads=0,\n",
+    "        inter_op_num_threads=0,\n",
+    "        providers=None,\n",
+    "        batch_size=1,\n",
+    "    )\n",
+    "\n",
+    "    # access cloud and cloud shadow mask\n",
+    "    csmask_csm = csmask.csm\n",
+    "\n",
+    "    # access valid mask\n",
+    "    csmask_valid = csmask.valid"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

tests/test_mask.py

@@ -14,11 +14,11 @@
 @pytest.mark.parametrize(
     "img, band_order, nodata_value",
     [
-        (np.empty((256, 256, 6), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"], None),
-        (np.empty((256, 256, 8), np.float32), ["Green", "Red", "Blue", "NIR", "SWIR1", "SWIR2", "NIR2", "ETC"], None),
-        (np.empty((256, 256, 6), np.float32), ["Green", "Red", "Blue", "NIR", "SWIR1", "SWIR2"], -666),
+        (np.empty((256, 256, 6), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"], None),
+        (np.empty((256, 256, 8), np.float32), ["Green", "Red", "Blue", "NIR", "SWIR16", "SWIR22", "NIR2", "ETC"], None),
+        (np.empty((256, 256, 6), np.float32), ["Green", "Red", "Blue", "NIR", "SWIR16", "SWIR22"], -666),
         (np.empty((256, 256, 4), np.float32), ["Red", "Green", "Blue", "NIR"], 0),
-        (np.empty((256, 256, 5), np.float32), ["Red", "Green", "Blue", "NIR", "SWIR2"], 0),
+        (np.empty((256, 256, 5), np.float32), ["Red", "Green", "Blue", "NIR", "SWIR22"], 0),
     ],
 )
 def test_csmask_init(img, band_order, nodata_value):
@@ -28,11 +28,11 @@ def test_csmask_init(img, band_order, nodata_value):
 @pytest.mark.parametrize(
     "img, band_order, nodata_value",
     [
-        (3, ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"], None),
-        (np.empty((256, 256), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"], None),
-        (np.empty((256, 256, 3), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"], None),
-        (np.empty((256, 256, 6), np.uint8), ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"], None),
-        (np.empty((256, 256, 6), np.float32), ["Blue", "Green", "Yellow", "NIR", "SWIR1", "SWIR2"], None),
+        (3, ["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"], None),
+        (np.empty((256, 256), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"], None),
+        (np.empty((256, 256, 3), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"], None),
+        (np.empty((256, 256, 6), np.uint8), ["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"], None),
+        (np.empty((256, 256, 6), np.float32), ["Blue", "Green", "Yellow", "NIR", "SWIR16", "SWIR22"], None),
         (np.empty((256, 256, 6), np.float32), None, None),
     ],
 )
@@ -44,8 +44,8 @@ def test_csmask_init_raises(img, band_order, nodata_value):
 @pytest.mark.parametrize(
     "img, band_order, nodata_value",
     [
-        (np.empty((128, 128, 6), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"], None),
-        (np.empty((64, 64, 6), np.float32), ["Green", "Red", "Blue", "NIR", "SWIR1", "SWIR2"], -666),
+        (np.empty((128, 128, 6), np.float32), ["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"], None),
+        (np.empty((64, 64, 6), np.float32), ["Green", "Red", "Blue", "NIR", "SWIR16", "SWIR22"], -666),
     ],
 )
 def test_csmask_init_warns(img, band_order, nodata_value):
@@ -64,7 +64,7 @@ def test_csmask_init_warns(img, band_order, nodata_value):
     ],
 )
 def test_csmask_csm_6band(data):
-    csmask = CSmask(img=data["img"], band_order=["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"])
+    csmask = CSmask(img=data["img"], band_order=["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"])
     y_pred = csmask.csm
     y_true = reclassify(data["msk"], {"reclass_value_from": [0, 1, 2, 3, 4], "reclass_value_to": [2, 0, 0, 0, 1]})
     y_true = y_true.ravel()
@@ -84,7 +84,7 @@ def test_csmask_csm_6band(data):
     ],
 )
 def test_csmask_valid_6band(data):
-    csmask = CSmask(img=data["img"], band_order=["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"])
+    csmask = CSmask(img=data["img"], band_order=["Blue", "Green", "Red", "NIR", "SWIR16", "SWIR22"])
     y_pred = csmask.valid
     y_true = reclassify(data["msk"], {"reclass_value_from": [0, 1, 2, 3, 4], "reclass_value_to": [0, 1, 1, 1, 0]})
     y_true_inverted = ~y_true.astype(bool)

ukis_csmask/__init__.py

@@ -1 +1 @@
-__version__ = "0.2.2"
+__version__ = "1.0.0"