|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "metadata": {}, |
| 6 | + "source": [ |
| 7 | + "# Segment clouds and cloud shadows in Landsat-8 images (L1C)\n", |
| 8 | + "This notebook shows an example on how to use [ukis-csmask](https://github.com/dlr-eoc/ukis-csmask) to segment clouds and cloud shadows in Level-1C images from Landsat-8. Images are loaded from local file system. 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.\n", |
| 9 | + "\n", |
| 10 | + "> NOTE: to run this notebook, we first need to install some additional dependencies for image handling\n", |
| 11 | + "```shell\n", |
| 12 | + "pip install ukis-pysat[complete]\n", |
| 13 | + "```" |
| 14 | + ] |
| 15 | + }, |
| 16 | + { |
| 17 | + "cell_type": "code", |
| 18 | + "execution_count": null, |
| 19 | + "id": "703f3744-902d-470b-a80f-9a8d3ea08dfa", |
| 20 | + "metadata": {}, |
| 21 | + "outputs": [], |
| 22 | + "source": [ |
| 23 | + "import rasterio\n", |
| 24 | + "import numpy as np\n", |
| 25 | + "\n", |
| 26 | + "from pathlib import Path\n", |
| 27 | + "from ukis_csmask.mask import CSmask\n", |
| 28 | + "from ukis_pysat.raster import Image, Platform" |
| 29 | + ] |
| 30 | + }, |
| 31 | + { |
| 32 | + "cell_type": "code", |
| 33 | + "execution_count": null, |
| 34 | + "id": "c9cd86e4", |
| 35 | + "metadata": {}, |
| 36 | + "outputs": [], |
| 37 | + "source": [ |
| 38 | + "# user settings\n", |
| 39 | + "data_path = \"/your_data_path/\"\n", |
| 40 | + "L8_file_prefix = \"LC08_L1TP_191015_20210428_20210507_02_T1\"\n", |
| 41 | + "product_level = \"l1c\"\n", |
| 42 | + "band_order = [\"blue\", \"green\", \"red\", \"nir\", \"swir16\", \"swir22\"]\n", |
| 43 | + "providers = [\"CUDAExecutionProvider\"]\n", |
| 44 | + "out_dir = \"ukis-csmask/examples\"" |
| 45 | + ] |
| 46 | + }, |
| 47 | + { |
| 48 | + "cell_type": "code", |
| 49 | + "execution_count": null, |
| 50 | + "id": "8ca03c78-1e24-479c-9786-a1b43206a08b", |
| 51 | + "metadata": {}, |
| 52 | + "outputs": [], |
| 53 | + "source": [ |
| 54 | + "# set Landsat 8 source path and prefix (example)\n", |
| 55 | + "data_path = data_path + L8_file_prefix + \"/\"\n", |
| 56 | + "mtl_file = data_path + L8_file_prefix + \"_MTL.txt\"\n", |
| 57 | + "\n", |
| 58 | + "# stack [B2:'Blue', B3:'Green', B4:'Red', B5:'NIR', B6:'SWIR1', B7:'SWIR2'] as numpy array\n", |
| 59 | + "L8_band_files = [data_path + L8_file_prefix + \"_B\" + x + \".TIF\" for x in [str(x + 2) for x in range(6)]]\n", |
| 60 | + "\n", |
| 61 | + "# >> adopted from https://gis.stackexchange.com/questions/223910/using-rasterio-or-gdal-to-stack-multiple-bands-without-using-subprocess-commands\n", |
| 62 | + "# read metadata of first file\n", |
| 63 | + "with rasterio.open(L8_band_files[0]) as src0:\n", |
| 64 | + " meta = src0.meta\n", |
| 65 | + "# update meta to reflect the number of layers\n", |
| 66 | + "meta.update(count=len(L8_band_files))\n", |
| 67 | + "# read each layer and append it to numpy array\n", |
| 68 | + "L8_bands = []\n", |
| 69 | + "for id, layer in enumerate(L8_band_files, start=1):\n", |
| 70 | + " with rasterio.open(layer) as src1:\n", |
| 71 | + " L8_bands.append(src1.read(1))\n", |
| 72 | + "L8_bands = np.stack(L8_bands, axis=2)\n", |
| 73 | + "# <<\n", |
| 74 | + "\n", |
| 75 | + "img = Image(data=L8_bands, crs=meta[\"crs\"], transform=meta[\"transform\"], dimorder=\"last\")\n", |
| 76 | + "img.dn2toa(platform=Platform.Landsat8, mtl_file=mtl_file, wavelengths=band_order)\n", |
| 77 | + "img.warp(resampling_method=0, resolution=30, dst_crs=img.dataset.crs)" |
| 78 | + ] |
| 79 | + }, |
| 80 | + { |
| 81 | + "cell_type": "code", |
| 82 | + "execution_count": null, |
| 83 | + "id": "7b568942-84e6-4baf-b490-a213b3787f80", |
| 84 | + "metadata": {}, |
| 85 | + "outputs": [], |
| 86 | + "source": [ |
| 87 | + "# compute cloud and cloud shadow mask\n", |
| 88 | + "csmask = CSmask(\n", |
| 89 | + " img=img.arr,\n", |
| 90 | + " band_order=band_order,\n", |
| 91 | + " product_level=product_level,\n", |
| 92 | + " nodata_value=0,\n", |
| 93 | + " invalid_buffer=4,\n", |
| 94 | + " intra_op_num_threads=0,\n", |
| 95 | + " inter_op_num_threads=0,\n", |
| 96 | + " providers=providers,\n", |
| 97 | + " batch_size=1,\n", |
| 98 | + ")\n", |
| 99 | + "\n", |
| 100 | + "# access cloud and cloud shadow mask as numpy array\n", |
| 101 | + "csm = csmask.csm\n", |
| 102 | + "\n", |
| 103 | + "# access valid mask as numpy array\n", |
| 104 | + "valid = csmask.valid" |
| 105 | + ] |
| 106 | + }, |
| 107 | + { |
| 108 | + "cell_type": "code", |
| 109 | + "execution_count": null, |
| 110 | + "id": "68eb9c30-06f7-409e-914d-21e00f45de99", |
| 111 | + "metadata": {}, |
| 112 | + "outputs": [], |
| 113 | + "source": [ |
| 114 | + "# convert results to ukis-pysat Image\n", |
| 115 | + "# this assigns back the georeference\n", |
| 116 | + "csm = Image(csm, transform=img.dataset.transform, crs=img.dataset.crs, dimorder=\"last\")\n", |
| 117 | + "valid = Image(valid, transform=img.dataset.transform, crs=img.dataset.crs, dimorder=\"last\")\n", |
| 118 | + "\n", |
| 119 | + "# write results to file\n", |
| 120 | + "csm.write_to_file(\n", |
| 121 | + " path_to_file=Path(out_dir) / Path(f\"{L8_file_prefix}_csm.tif\"),\n", |
| 122 | + " dtype=csm.dtype,\n", |
| 123 | + " driver=\"COG\",\n", |
| 124 | + " compress=\"LZW\",\n", |
| 125 | + " kwargs={\"BLOCKSIZE\": 512, \"BIGTIFF\": \"IF_SAFER\"},\n", |
| 126 | + ")\n", |
| 127 | + "valid.write_to_file(\n", |
| 128 | + " path_to_file=Path(out_dir) / Path(f\"{L8_file_prefix}_valid.tif\"),\n", |
| 129 | + " dtype=valid.dtype,\n", |
| 130 | + " driver=\"COG\",\n", |
| 131 | + " compress=\"LZW\",\n", |
| 132 | + " kwargs={\"BLOCKSIZE\": 512, \"BIGTIFF\": \"IF_SAFER\"},\n", |
| 133 | + ")" |
| 134 | + ] |
| 135 | + } |
| 136 | + ], |
| 137 | + "metadata": { |
| 138 | + "kernelspec": { |
| 139 | + "display_name": "Python 3 (ipykernel)", |
| 140 | + "language": "python", |
| 141 | + "name": "python3" |
| 142 | + }, |
| 143 | + "language_info": { |
| 144 | + "codemirror_mode": { |
| 145 | + "name": "ipython", |
| 146 | + "version": 3 |
| 147 | + }, |
| 148 | + "file_extension": ".py", |
| 149 | + "mimetype": "text/x-python", |
| 150 | + "name": "python", |
| 151 | + "nbconvert_exporter": "python", |
| 152 | + "pygments_lexer": "ipython3", |
| 153 | + "version": "3.11.10" |
| 154 | + } |
| 155 | + }, |
| 156 | + "nbformat": 4, |
| 157 | + "nbformat_minor": 5 |
| 158 | +} |
0 commit comments