diff --git a/.gitignore b/.gitignore index da270fc..718df2f 100644 --- a/.gitignore +++ b/.gitignore @@ -246,3 +246,17 @@ tags *tif* integration_test/ chips +*.h5 + +*.cpg +*.dbf +*.prj +*.sbn +*.sbx +*.shp +*.shx +*.zip +*.qgz +hwds + +notebooks/data diff --git a/README.md b/README.md index 19d7c7d..3a05206 100644 --- a/README.md +++ b/README.md @@ -1,77 +1,30 @@ # SatChip -A package for satellite image AI data prep. This package "chips" data labels and satellite imagery into 264x264 image arrays following the TerraMind extension of the MajorTom specification. +A package for satellite image AI data prep. ## Usage -`SatChip` relies on a two-step process; chip your label train data inputs, then create corresponding chips for different remote sensing data sources. -### Step 1: Chip labels -The `chiplabel` CLI tool takes a GDAL-compatible image, a collection date, and an optional chip directory as input using the following format: +```python +import satchip -```bash -chiplabel PATH/TO/LABELS.tif DATE(UTC FORMAT) --chipdir CHIP_DIR -``` -For example: -```bash -chiplabel LA_damage_20250113_v0.tif 2024-01-01T01:01:01 --chipdir chips -``` -This will produce an output zipped Zarr store label dataset with the name `{LABEL}_{SAMPLE}.zarr.zip` (see the (Tiling Schema)[#tiling_schema] section for details on the `SAMPLE` name) to the `LABEL` directory in the specified chip directory (`--chipdir`). This file will be the input to the remote sensing data chipping step. - -For more information on usage see `chiplabel --help` - -### Step 2: Chip remote sensing data -The `chipdata` CLI tool takes a path to a directory containing chip labels, a dataset name, a date range and a set of optional parameters using the following format: -```bash -chipdata PATH/TO/LABEL DATASET Ymd-Ymd \ - --maxcloudpct MAX_CLOUD_PCT --strategy STRATEGY \ - --chipdir CHIPPUT_DIR --imagedir IMAGE_DIR -``` -For example: -```bash -chipdata LABEL S2L2A 20250112-20250212 --maxcloudpct 20 --chipdir CHIP_DIR --imagedir IMAGES -``` -Similarly to step 1, this will produce an output zipped Zarr store that contains chipped data for your chosen dataset with the name `{LABELS_{SAMPLE}_{DATASET}.zarr.zip`. The arguments are as follows: -- `PATH/TO/LABEL`: the path to your training labels -- `DATASET`: The satellite imagery dataset you would like to create labels for. See the list below for all current options. -- `Ymd-Ymd`: The date range to select imagery from. For example, `20250112-20250212` selects imagery between January 12 and February 12, 2025. -- `MAX_CLOUD_PCT`: For optical data, this optional parameter lets you set the maximum amount of cloud coverage allowed in a chip. Values between 0 and 100 are allowed. Cloud coverage is calculated on a per-chip basis. The default is 100 i.e., no limit. -- `STRATEGY`: Lets you selected what data inside your date range will be used to create chips. Specifying `BEST` (the default) will create a chip for the image closest to the beginning of your date range that has at least 95% spatial coverage. Specifying `ALL` will create chips for all images within your date range that have at least 95% spatial coverage. -- `CHIP_DIR`: Specifies the directory where the image chips will be saved. If not specified, this defaults to your current directory. -- `IMAGE_DIR`: Specifies the directory where the full-size satellite images will be downloaded to. If this argument is not provided, the images will be stored in the `IMAGES` directory within `CHIP_DIR`. - -Currently supported datasets include: -- `S2L2A`: Sentinel-2 L2A data sourced from the [Sentinel-2 AWS Open Data Archive](https://registry.opendata.aws/sentinel-2/) -- `HLS`: Harmonized Landsat Sentinel-2 data sourced from [LP DAAC's Data Archive](https://www.earthdata.nasa.gov/data/projects/hls) -- `S1RTC`: OPERA Sentinel-1 Radiometric Terrain Corrected (RTC) data from [ASF DAAC's Data Archive](https://www.jpl.nasa.gov/go/opera/products/rtc-product/) -- `HYP3S1RTC`: Sentinel-1 Radiometric Terrain Corrected (RTC) data created using [ASF's HyP3 on-demand platform](https://hyp3-docs.asf.alaska.edu/guides/rtc_product_guide/) - -## Tiling Schema +data_paths = { + MODALITY: hwds_path / MODALITY, + "RAW": modality_path / "RAW", + "WGS84": modality_path / "WGS84", + "MERGE": modality_path / "MERGE", + "CHIPS": modality_path / "CHIPS", + "CHIPS_TM": modality_path / "CHIPS_TM", + "PLOTS": modality_path / "PLOTS", + "SPLITS": modality_path / "SPLITS", +} -This package chips images based on the [TerraMesh grid system](https://huggingface.co/datasets/ibm-esa-geospatial/TerraMesh), which builds on the [MajorTOM grid system](https://github.com/ESA-PhiLab/Major-TOM). +modalities = ['HLS'] -The MajorTOM grid system provides a global set of fixed image grids that are 1068x1068 pixels in size. A MajorTOM grid can be defined for any tile size, but we fix the grid to 10x10 Km tiles. Tiles are named using the format: -``` -ROW[U|D]_COL[L|R] -``` -Where, `ROW` is indexed from the equator, with a suffix `U` (up) for tiles north of the equator and `D` (down) for tiles south of it, and `COL` is indexed from the prime meridian, with a suffix `L` (left) for tiles east of the prime meridian and `R` (right) for tiles west of it. - -To support finer subdivisions, the TerraMesh grid system divides each MajorTOM grid into a 4x4 set of sub-tiles, each 264x264 pixels. The subgrid is centered within the parent tile, leaving a 6-pixel border around each sub-tile. Subgrid names extend the base format with two additional indices: -``` -ROW[U|D]_COL[L|R]_SUBCOL_SUBROW -``` -For instance, the bottom-left subgrid of MajorTOM tile `434U_876L` is named `434U_876L_0_3`. See the figure below for a visual description: - -![TerraMesh tiling schema](assets/satchip_schema.svg) +data = satchip.find_data(area, modality) +reprojected_data = satchip.repoject(raw_data, projection='WGS84', modality=) # stack bands and mosaic +mosaics = satchip.mosaic(data, stack_bands=True) # stack bands and mosaic +masks = satchip.generate_masks(mosaics) -## Viewing Chips -Assessing chips after their creation can be challenging due to the large number of small images created. To address this issue, SatChip includes a `chipview` CLI tool that uses Matplotlib to quickly visualize the data included within the created zipped Zarr stores: -```bash -chipview PATH/TO/CHIP.zarr.zip --band BAND +chips = satchip.chip_data(mosaics, masks) +chips = satchip.filter_chips(chips) ``` -Where `PATH/TO/CHIPS.zarr.zip` is the path to the chip file (labels or image data), and `BAND` is an OPTIONAL name of the band you would like to view. If no band is specified, an OPERA-style RGB decomposition will be shown for RTC data, and an RGB composite will be shown for optical data. - -## License -`SatChip` is licensed under the BSD-3-Clause open source license. See the LICENSE file for more details. - -## Contributing -Contributions to the `SatChip` are welcome! If you would like to contribute, please submit a pull request on the GitHub repository. diff --git a/notebooks/hwds-example.ipynb b/notebooks/hwds-example.ipynb new file mode 100644 index 0000000..51d8b10 --- /dev/null +++ b/notebooks/hwds-example.ipynb @@ -0,0 +1,1970 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4dcd6420", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/wbhorn/miniforge3/envs/satchip/lib/python3.14/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", + " from .autonotebook import tqdm as notebook_tqdm\n" + ] + } + ], + "source": [ + "from pathlib import Path\n", + "import shutil\n", + "\n", + "import geopandas as gpd\n", + "import pandas as pd\n", + "\n", + "from satchip import models, download_data, merge_modality, generate_labels, chip_data, view" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "2b0886a1", + "metadata": {}, + "outputs": [], + "source": [ + "shp_path = 'Pristine_Merged'\n", + "df = gpd.read_file(shp_path)\n", + "\n", + "df['SwathDate'] = pd.to_datetime(df['SwathDate'], format='%Y-%m-%d')\n", + "df['HLSDate'] = pd.to_datetime(df['HLSDate'], format='%Y-%m-%d')\n", + "\n", + "#df = df[df['Visibility'] == '1']" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "55da8795", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
SwathDateHLSIDMODISIDVisibilityHLSDateFieldgeometry
02019-07-09102a10212019-07-17S2POLYGON Z ((-99.49952 40.26693 0, -99.47639 40...
12019-08-13106a10622019-08-28LSPOLYGON Z ((-101.59283 40.78718 0, -101.60341 ...
22018-08-06126a12612018-08-11S2POLYGON Z ((-98.2416 41.16581 0, -98.24082 41....
32018-08-06127a12712018-08-11S2POLYGON Z ((-97.78549 41.00092 0, -97.76367 41...
42018-07-29129a12932018-08-07LSPOLYGON Z ((-103.1513 43.15378 0, -103.14676 4...
52018-07-27130a13022018-08-07LSPOLYGON Z ((-102.98845 43.56544 0, -102.97929 ...
62018-07-27132a13232018-08-11LSMULTIPOLYGON Z (((-97.42377 41.60069 0, -97.42...
72018-07-18133a13322018-07-26LSPOLYGON Z ((-97.5217 42.65295 0, -97.5154 42.6...
82018-06-30134a13412018-07-12S2POLYGON Z ((-99.72435 40.41203 0, -99.70477 40...
92017-07-05598a59822017-07-15S2POLYGON Z ((-98.60738 45.24939 0, -98.60167 45...
102018-06-26614a61132018-07-08LSMULTIPOLYGON Z (((-100.65633 44.66629 0, -100....
112019-08-06623d62422019-08-19LSPOLYGON Z ((-100.00341 45.55568 0, -99.98419 4...
122019-08-06623a62312019-08-19LSPOLYGON Z ((-100.33919 46.03454 0, -100.33567 ...
132018-07-02889a88922018-07-13S2POLYGON Z ((-102.70434 46.61383 0, -102.68865 ...
142018-07-02889b89022018-07-13S2POLYGON Z ((-101.87838 46.45158 0, -101.8724 4...
152018-06-26892a89222018-07-08S2POLYGON Z ((-103.21657 46.25444 0, -103.21434 ...
162017-07-20915a91532017-07-26LSPOLYGON Z ((-102.58639 47.46329 0, -102.55608 ...
172019-07-291079a107922019-08-04S2POLYGON Z ((-102.59566 39.43458 0, -102.60361 ...
182019-08-131069a106932019-08-14S2POLYGON Z ((-102.15607 39.11807 0, -102.14741 ...
192018-07-101378a137832018-07-12LSPOLYGON Z ((-95.05656 43.17247 0, -95.04818 43...
202019-07-03628a62832019-07-13S2POLYGON Z ((-102.31279 43.36328 0, -102.30699 ...
212020-06-04638a63832020-06-12S2POLYGON Z ((-103.23338 45.55551 0, -103.21165 ...
222020-06-07116d63932020-06-17S2POLYGON Z ((-101.18151 43.22672 0, -101.17656 ...
232020-06-07116e64022020-06-17S2POLYGON Z ((-100.67301 43.4943 0, -100.68216 4...
242020-07-06648a64822020-07-16S2POLYGON Z ((-97.35114 43.00508 0, -97.34554 43...
252017-08-101052a105222017-08-13LSPOLYGON Z ((-102.68889 38.15458 0, -102.68788 ...
262018-07-261055a105522018-08-04S2POLYGON Z ((-102.09648 39.57415 0, -102.09107 ...
272018-07-281056a105632018-08-07S2MULTIPOLYGON Z (((-101.94001 40.61807 0, -101....
282018-07-281347a134722018-08-04S2POLYGON Z ((-100.73413 39.35341 0, -100.73135 ...
292018-06-191064a106422018-07-03S2POLYGON Z ((-102.42969 39.8437 0, -102.42208 3...
302019-08-241338a133832019-09-05S2POLYGON Z ((-100.22342 39.23137 0, -100.20831 ...
312018-07-271344a134422018-08-02S2POLYGON Z ((-99.62458 38.06997 0, -99.61494 38...
322020-08-101373a137312020-08-17S2POLYGON Z ((-94.97657 42.20582 0, -94.97333 42...
332020-08-101373d137332020-08-17S2POLYGON Z ((-94.80415 42.01749 0, -94.78768 42...
342020-08-101373e137332020-08-17S2POLYGON Z ((-94.55871 41.9804 0, -94.53908 41....
352020-07-111376a137612020-07-19LSPOLYGON Z ((-92.15179 43.02718 0, -92.1428 43....
362019-08-22110d11032019-09-05S2MULTIPOLYGON Z (((-98.39133 40.49204 0, -98.39...
\n", + "
" + ], + "text/plain": [ + " SwathDate HLSID MODISID Visibility HLSDate Field \\\n", + "0 2019-07-09 102a 102 1 2019-07-17 S2 \n", + "1 2019-08-13 106a 106 2 2019-08-28 LS \n", + "2 2018-08-06 126a 126 1 2018-08-11 S2 \n", + "3 2018-08-06 127a 127 1 2018-08-11 S2 \n", + "4 2018-07-29 129a 129 3 2018-08-07 LS \n", + "5 2018-07-27 130a 130 2 2018-08-07 LS \n", + "6 2018-07-27 132a 132 3 2018-08-11 LS \n", + "7 2018-07-18 133a 133 2 2018-07-26 LS \n", + "8 2018-06-30 134a 134 1 2018-07-12 S2 \n", + "9 2017-07-05 598a 598 2 2017-07-15 S2 \n", + "10 2018-06-26 614a 611 3 2018-07-08 LS \n", + "11 2019-08-06 623d 624 2 2019-08-19 LS \n", + "12 2019-08-06 623a 623 1 2019-08-19 LS \n", + "13 2018-07-02 889a 889 2 2018-07-13 S2 \n", + "14 2018-07-02 889b 890 2 2018-07-13 S2 \n", + "15 2018-06-26 892a 892 2 2018-07-08 S2 \n", + "16 2017-07-20 915a 915 3 2017-07-26 LS \n", + "17 2019-07-29 1079a 1079 2 2019-08-04 S2 \n", + "18 2019-08-13 1069a 1069 3 2019-08-14 S2 \n", + "19 2018-07-10 1378a 1378 3 2018-07-12 LS \n", + "20 2019-07-03 628a 628 3 2019-07-13 S2 \n", + "21 2020-06-04 638a 638 3 2020-06-12 S2 \n", + "22 2020-06-07 116d 639 3 2020-06-17 S2 \n", + "23 2020-06-07 116e 640 2 2020-06-17 S2 \n", + "24 2020-07-06 648a 648 2 2020-07-16 S2 \n", + "25 2017-08-10 1052a 1052 2 2017-08-13 LS \n", + "26 2018-07-26 1055a 1055 2 2018-08-04 S2 \n", + "27 2018-07-28 1056a 1056 3 2018-08-07 S2 \n", + "28 2018-07-28 1347a 1347 2 2018-08-04 S2 \n", + "29 2018-06-19 1064a 1064 2 2018-07-03 S2 \n", + "30 2019-08-24 1338a 1338 3 2019-09-05 S2 \n", + "31 2018-07-27 1344a 1344 2 2018-08-02 S2 \n", + "32 2020-08-10 1373a 1373 1 2020-08-17 S2 \n", + "33 2020-08-10 1373d 1373 3 2020-08-17 S2 \n", + "34 2020-08-10 1373e 1373 3 2020-08-17 S2 \n", + "35 2020-07-11 1376a 1376 1 2020-07-19 LS \n", + "36 2019-08-22 110d 110 3 2019-09-05 S2 \n", + "\n", + " geometry \n", + "0 POLYGON Z ((-99.49952 40.26693 0, -99.47639 40... \n", + "1 POLYGON Z ((-101.59283 40.78718 0, -101.60341 ... \n", + "2 POLYGON Z ((-98.2416 41.16581 0, -98.24082 41.... \n", + "3 POLYGON Z ((-97.78549 41.00092 0, -97.76367 41... \n", + "4 POLYGON Z ((-103.1513 43.15378 0, -103.14676 4... \n", + "5 POLYGON Z ((-102.98845 43.56544 0, -102.97929 ... \n", + "6 MULTIPOLYGON Z (((-97.42377 41.60069 0, -97.42... \n", + "7 POLYGON Z ((-97.5217 42.65295 0, -97.5154 42.6... \n", + "8 POLYGON Z ((-99.72435 40.41203 0, -99.70477 40... \n", + "9 POLYGON Z ((-98.60738 45.24939 0, -98.60167 45... \n", + "10 MULTIPOLYGON Z (((-100.65633 44.66629 0, -100.... \n", + "11 POLYGON Z ((-100.00341 45.55568 0, -99.98419 4... \n", + "12 POLYGON Z ((-100.33919 46.03454 0, -100.33567 ... \n", + "13 POLYGON Z ((-102.70434 46.61383 0, -102.68865 ... \n", + "14 POLYGON Z ((-101.87838 46.45158 0, -101.8724 4... \n", + "15 POLYGON Z ((-103.21657 46.25444 0, -103.21434 ... \n", + "16 POLYGON Z ((-102.58639 47.46329 0, -102.55608 ... \n", + "17 POLYGON Z ((-102.59566 39.43458 0, -102.60361 ... \n", + "18 POLYGON Z ((-102.15607 39.11807 0, -102.14741 ... \n", + "19 POLYGON Z ((-95.05656 43.17247 0, -95.04818 43... \n", + "20 POLYGON Z ((-102.31279 43.36328 0, -102.30699 ... \n", + "21 POLYGON Z ((-103.23338 45.55551 0, -103.21165 ... \n", + "22 POLYGON Z ((-101.18151 43.22672 0, -101.17656 ... \n", + "23 POLYGON Z ((-100.67301 43.4943 0, -100.68216 4... \n", + "24 POLYGON Z ((-97.35114 43.00508 0, -97.34554 43... \n", + "25 POLYGON Z ((-102.68889 38.15458 0, -102.68788 ... \n", + "26 POLYGON Z ((-102.09648 39.57415 0, -102.09107 ... \n", + "27 MULTIPOLYGON Z (((-101.94001 40.61807 0, -101.... \n", + "28 POLYGON Z ((-100.73413 39.35341 0, -100.73135 ... \n", + "29 POLYGON Z ((-102.42969 39.8437 0, -102.42208 3... \n", + "30 POLYGON Z ((-100.22342 39.23137 0, -100.20831 ... \n", + "31 POLYGON Z ((-99.62458 38.06997 0, -99.61494 38... \n", + "32 POLYGON Z ((-94.97657 42.20582 0, -94.97333 42... \n", + "33 POLYGON Z ((-94.80415 42.01749 0, -94.78768 42... \n", + "34 POLYGON Z ((-94.55871 41.9804 0, -94.53908 41.... \n", + "35 POLYGON Z ((-92.15179 43.02718 0, -92.1428 43.... \n", + "36 MULTIPOLYGON Z (((-98.39133 40.49204 0, -98.39... " + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "f1a32464", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'HLS_S30': {'modality': {'id': 'HLS_S30',\n", + " 'collection': 'HLSS30',\n", + " 'bands': (Band(id='B02', name='Blue', shortname='B'),\n", + " Band(id='B03', name='Green', shortname='G'),\n", + " Band(id='B04', name='Red', shortname='R'),\n", + " Band(id='B8A', name='NIR Narrow', shortname='N'),\n", + " Band(id='B11', name='SWIR 1', shortname='SW1'),\n", + " Band(id='B12', name='SWIR 2', shortname='SW2'),\n", + " Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))},\n", + " 'raw': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/raw'),\n", + " 'merged': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/merged'),\n", + " 'wgs84': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84'),\n", + " 'stacked': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/stacked'),\n", + " 'warped': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/warped'),\n", + " 'chips': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/chips'),\n", + " 'plots': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/plots')},\n", + " 'HLS_L30': {'modality': {'id': 'HLS_L30',\n", + " 'collection': 'HLSL30',\n", + " 'bands': (Band(id='B02', name='Blue', shortname='B'),\n", + " Band(id='B03', name='Green', shortname='G'),\n", + " Band(id='B04', name='Red', shortname='R'),\n", + " Band(id='B05', name='NIR Narrow', shortname='N'),\n", + " Band(id='B06', name='SWIR 1', shortname='SW1'),\n", + " Band(id='B07', name='SWIR 2', shortname='SW2'),\n", + " Band(id='Fmask', name='Cloud Mask', shortname='fmask'))},\n", + " 'raw': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/raw'),\n", + " 'merged': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged'),\n", + " 'wgs84': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84'),\n", + " 'stacked': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/stacked'),\n", + " 'warped': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/warped'),\n", + " 'chips': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/chips'),\n", + " 'plots': PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/plots')}}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DATA_PATH = Path.cwd() / 'data'\n", + "\n", + "def make_mod_paths(modality):\n", + " base_path = DATA_PATH / modality['id']\n", + "\n", + " return {\n", + " 'modality': modality,\n", + " 'raw': base_path / 'raw',\n", + " 'merged': base_path / 'merged',\n", + " 'wgs84': base_path / 'wgs84',\n", + " 'stacked': base_path / 'stacked',\n", + " 'warped': base_path / 'warped',\n", + " 'chips': base_path / 'chips',\n", + " 'plots': base_path / 'plots'\n", + " }\n", + "\n", + "mod_paths = {\n", + " modality['id']: make_mod_paths(modality) for modality in (models.HLS_S30, models.HLS_L30)\n", + "}\n", + "mod_paths" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "0d0202e9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))} /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/raw\n", + "Logging in to earthaccess\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 13296.49it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 371908.73it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 725937.23it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 102a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/102a.MASK.tif\n", + "Adding event files for 102a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 11018.31it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 324023.48it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 18/18 [00:00<00:00, 431414.13it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 126a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/126a.MASK.tif\n", + "Adding event files for 126a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 12883.53it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 474826.87it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 848286.20it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 127a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/127a.MASK.tif\n", + "Adding event files for 127a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 14230.04it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 547083.13it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 811800.77it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 134a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/134a.MASK.tif\n", + "Adding event files for 134a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 13996.57it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 453438.27it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 715615.85it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 598a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/598a.MASK.tif\n", + "Adding event files for 598a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 14614.30it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 401582.30it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 695829.24it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 889a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/889a.MASK.tif\n", + "Adding event files for 889a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 13644.94it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 461758.24it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 834226.21it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 889b\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/889b.MASK.tif\n", + "Adding event files for 889b\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 26597.66it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 677107.37it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 72/72 [00:00<00:00, 1263556.02it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 892a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/892a.MASK.tif\n", + "Adding event files for 892a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 162/162 [00:00<00:00, 50027.78it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 162/162 [00:00<00:00, 109610.78it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 162/162 [00:00<00:00, 1085426.91it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1079a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1079a.MASK.tif\n", + "Adding event files for 1079a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 13109.48it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 531672.34it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 904161.34it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1069a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1069a.MASK.tif\n", + "Adding event files for 1069a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 28225.99it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 702302.07it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 72/72 [00:00<00:00, 1090216.20it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 628a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/628a.MASK.tif\n", + "Adding event files for 628a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 90/90 [00:00<00:00, 28299.52it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 90/90 [00:00<00:00, 768813.36it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 90/90 [00:00<00:00, 1362770.25it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 638a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/638a.MASK.tif\n", + "Adding event files for 638a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 25067.64it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 610080.58it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 72/72 [00:00<00:00, 1212810.80it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 116d\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/116d.MASK.tif\n", + "Adding event files for 116d\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 12609.18it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 533551.04it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 461758.24it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 116e\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/116e.MASK.tif\n", + "Adding event files for 116e\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 11193.10it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 357807.92it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 18/18 [00:00<00:00, 426539.39it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 648a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/648a.MASK.tif\n", + "Adding event files for 648a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 108/108 [00:00<00:00, 39645.09it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 108/108 [00:00<00:00, 120731.57it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 108/108 [00:00<00:00, 1438047.09it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1055a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1055a.MASK.tif\n", + "Adding event files for 1055a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 27267.71it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 72/72 [00:00<00:00, 774333.05it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 72/72 [00:00<00:00, 1247892.10it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1056a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1056a.MASK.tif\n", + "Adding event files for 1056a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 13693.20it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 517105.97it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 848286.20it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1347a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1347a.MASK.tif\n", + "Adding event files for 1347a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 90/90 [00:00<00:00, 31324.15it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 90/90 [00:00<00:00, 498662.30it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 90/90 [00:00<00:00, 1310720.00it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1064a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1064a.MASK.tif\n", + "Adding event files for 1064a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 14260.95it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 601573.48it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 1027176.49it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1338a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1338a.MASK.tif\n", + "Adding event files for 1338a\n", + "no HLS_S30 data for 1344a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 646.08it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 235194.62it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 18/18 [00:00<00:00, 351151.03it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1373a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1373a.MASK.tif\n", + "Adding event files for 1373a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 10624.47it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 18/18 [00:00<00:00, 387166.52it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 18/18 [00:00<00:00, 585251.72it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1373d\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1373d.MASK.tif\n", + "Adding event files for 1373d\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 14273.08it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 431414.13it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 834226.21it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 1373e\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/1373e.MASK.tif\n", + "Adding event files for 1373e\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 13669.65it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 36/36 [00:00<00:00, 506694.44it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 36/36 [00:00<00:00, 555128.47it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_S30 data for 110d\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/wgs84/110d.MASK.tif\n", + "Adding event files for 110d\n", + "{'id': 'HLS_L30', 'collection': 'HLSL30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B05', name='NIR Narrow', shortname='N'), Band(id='B06', name='SWIR 1', shortname='SW1'), Band(id='B07', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='fmask'))} /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/raw\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 90/90 [00:00<00:00, 29262.59it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 90/90 [00:00<00:00, 778324.45it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 90/90 [00:00<00:00, 1297207.42it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 106a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/106a.MASK.tif\n", + "Adding event files for 106a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 12653.77it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 469511.64it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 30/30 [00:00<00:00, 781547.33it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 129a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/129a.MASK.tif\n", + "Adding event files for 129a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 22708.74it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 710898.98it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 60/60 [00:00<00:00, 1084733.79it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 130a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/130a.MASK.tif\n", + "Adding event files for 130a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 19331.56it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 700997.88it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 60/60 [00:00<00:00, 1048576.00it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 132a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/132a.MASK.tif\n", + "Adding event files for 132a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 24160.74it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 9597.22it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 30/30 [00:00<00:00, 758006.75it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 133a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/133a.MASK.tif\n", + "Adding event files for 133a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 33367.57it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 37803.55it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 60/60 [00:00<00:00, 1103764.21it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 614a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/614a.MASK.tif\n", + "Adding event files for 614a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 12575.37it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 443060.28it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 30/30 [00:00<00:00, 610820.97it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 623d\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/623d.MASK.tif\n", + "Adding event files for 623d\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 21969.29it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 60/60 [00:00<00:00, 706905.17it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 60/60 [00:00<00:00, 793874.57it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 623a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/623a.MASK.tif\n", + "Adding event files for 623a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 12608.13it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 430921.64it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 30/30 [00:00<00:00, 511500.49it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 915a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/915a.MASK.tif\n", + "Adding event files for 915a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 13312.43it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 400729.68it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 30/30 [00:00<00:00, 744551.01it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 1378a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/1378a.MASK.tif\n", + "Adding event files for 1378a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 13177.20it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 30/30 [00:00<00:00, 418036.94it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 30/30 [00:00<00:00, 748982.86it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 1052a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/1052a.MASK.tif\n", + "Adding event files for 1052a\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "QUEUEING TASKS | : 100%|██████████| 15/15 [00:00<00:00, 11670.29it/s]\n", + "PROCESSING TASKS | : 100%|██████████| 15/15 [00:00<00:00, 322638.77it/s]\n", + "COLLECTING RESULTS | : 100%|██████████| 15/15 [00:00<00:00, 491520.00it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found HLS_L30 data for 1376a\n", + "generated: /home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/wgs84/1376a.MASK.tif\n", + "Adding event files for 1376a\n" + ] + } + ], + "source": [ + "event_files = []\n", + "\n", + "for paths in mod_paths.values():\n", + " print(paths['modality'], paths['raw'])\n", + "\n", + " for idx, row in df.iterrows():\n", + " modality = paths['modality']\n", + " best_mod = models.HLS_L30 if row['Field'] == 'LS' else models.HLS_S30\n", + " if modality != best_mod:\n", + " continue\n", + "\n", + " damage_event = models.Event(\n", + " name=row['HLSID'], \n", + " date=row['HLSDate'], \n", + " wgs84_geometry=row['geometry'], \n", + " buffer_m=10000\n", + " )\n", + "\n", + " local = download_data.download_data(damage_event, modality, paths['raw'])\n", + "\n", + " if not local:\n", + " print(f'no {modality[\"id\"]} data for {damage_event.name}')\n", + " continue\n", + " else:\n", + " print(f'Found {modality[\"id\"]} data for {damage_event.name}')\n", + "\n", + " merged_event = merge_modality.merge_modality(\n", + " local, \n", + " modality, \n", + " event=damage_event, \n", + " output_path=paths['merged']\n", + " )\n", + "\n", + " stacked_filename = paths['stacked'] / f'{damage_event.name}.{modality[\"id\"]}.stacked.tif'\n", + " data_bands, fmask = merged_event[:-1], merged_event[-1]\n", + "\n", + " stacked = merge_modality.stack_bands(data_bands, stacked_filename)\n", + "\n", + " label = generate_labels.binary_mask_from_template(stacked, damage_event, paths['wgs84'])\n", + " mask, bands, fmask = merge_modality.warp_to_reference(\n", + " reference_path=label,\n", + " data_files=[stacked, fmask],\n", + " output_dir=paths['warped'],\n", + " bounding_box_wgs84=damage_event.buffered_geometry().bounds,\n", + " )\n", + " print(f'Adding event files for {damage_event.name}')\n", + " event_files.append((damage_event, modality, (mask, bands, fmask)))\n", + "\n", + " view.view_merged(\n", + " bands, \n", + " damage_event, \n", + " modality, \n", + " rgb_bands=[2, 1, 0],\n", + " quite=True, \n", + " save_to_file=paths['plots'] / f'{damage_event.name}.{modality[\"id\"]}.merged.plot.png'\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "a5c47f55", + "metadata": {}, + "source": [ + "Remove any bad aquisitions from `data/{modality}/plots` by remove the .png " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "ed4d4521", + "metadata": {}, + "outputs": [], + "source": [ + "def backup_plots():\n", + " for paths in mod_paths.values():\n", + " plots = paths['plots']\n", + " plots_backup = plots.parent / 'plots-backup'\n", + " shutil.copytree(plots, plots_backup, dirs_exist_ok=True)\n", + "\n", + "# backup_plots()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "519b634b", + "metadata": {}, + "outputs": [], + "source": [ + "def restore_plots():\n", + " for paths in mod_paths.values():\n", + " plots = paths['plots']\n", + " plots_backup = plots.parent / 'plots-backup'\n", + " shutil.copytree(plots_backup, plots, dirs_exist_ok=True)\n", + "\n", + "# restore_plots()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "44902869", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'HLS_S30': ['628a',\n", + " '1079a',\n", + " '1069a',\n", + " '116d',\n", + " '648a',\n", + " '638a',\n", + " '116e',\n", + " '1338a',\n", + " '598a',\n", + " '1347a',\n", + " '102a',\n", + " '110d',\n", + " '1373e',\n", + " '892a',\n", + " '1373d',\n", + " '1064a',\n", + " '1056a',\n", + " '889b',\n", + " '1055a',\n", + " '889a',\n", + " '126a',\n", + " '134a',\n", + " '127a',\n", + " '1373a'],\n", + " 'HLS_L30': ['132a',\n", + " '1378a',\n", + " '1376a',\n", + " '133a',\n", + " '915a',\n", + " '1052a',\n", + " '623d',\n", + " '130a',\n", + " '614a',\n", + " '106a',\n", + " '623a',\n", + " '129a']}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "keepers = {}\n", + "\n", + "for paths in mod_paths.values():\n", + " modality = paths['modality']\n", + " mod_keepers = [plot.name.split('.')[0] for plot in paths['plots'].glob('*.png')]\n", + " keepers[modality['id']] = mod_keepers\n", + "\n", + "keepers" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "6f59f9a1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(Event(name='102a', date=Timestamp('2019-07-17 00:00:00'), wgs84_geometry=, buffer_m=10000), {'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}, (PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/warped/102a.MASK.tif'), PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/warped/102a.HLS_S30.stacked.tif'), PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_S30/warped/102a.HLS_S30.2019-07-17.Fmask.tif')))\n" + ] + } + ], + "source": [ + "print(event_files[0])\n", + "evt, modality, (m, b, f) = event_files[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "19734cc8", + "metadata": {}, + "outputs": [], + "source": [ + "all_base = DATA_PATH / 'chips'\n", + "output_base = DATA_PATH / 'output'\n", + "\n", + "chip_paths = {\n", + " 'all': {\n", + " 'label': all_base / 'LABEL',\n", + " 'hls': all_base / 'HLS',\n", + " 'other': all_base / 'OTHER',\n", + " 'plots': all_base / 'PLOTS',\n", + " },\n", + " 'output': {\n", + " 'label': output_base / 'LABEL',\n", + " 'hls': output_base / 'HLS',\n", + " 'other': output_base / 'OTHER',\n", + " 'plots': all_base / 'PLOTS',\n", + " }\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "3c4cee14", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Chipping: 102a 102a.MASK.tif 102a.HLS_S30.stacked.tif 102a.HLS_S30.2019-07-17.Fmask.tif\n", + "Chipping: 126a 126a.MASK.tif 126a.HLS_S30.stacked.tif 126a.HLS_S30.2018-08-11.Fmask.tif\n", + "Chipping: 127a 127a.MASK.tif 127a.HLS_S30.stacked.tif 127a.HLS_S30.2018-08-11.Fmask.tif\n", + "Chipping: 134a 134a.MASK.tif 134a.HLS_S30.stacked.tif 134a.HLS_S30.2018-07-12.Fmask.tif\n", + "Chipping: 598a 598a.MASK.tif 598a.HLS_S30.stacked.tif 598a.HLS_S30.2017-07-15.Fmask.tif\n", + "Chipping: 889a 889a.MASK.tif 889a.HLS_S30.stacked.tif 889a.HLS_S30.2018-07-13.Fmask.tif\n", + "Chipping: 889b 889b.MASK.tif 889b.HLS_S30.stacked.tif 889b.HLS_S30.2018-07-13.Fmask.tif\n", + "Chipping: 892a 892a.MASK.tif 892a.HLS_S30.stacked.tif 892a.HLS_S30.2018-07-08.Fmask.tif\n", + "Chipping: 1079a 1079a.MASK.tif 1079a.HLS_S30.stacked.tif 1079a.HLS_S30.2019-08-04.Fmask.tif\n", + "Chipping: 1069a 1069a.MASK.tif 1069a.HLS_S30.stacked.tif 1069a.HLS_S30.2019-08-14.Fmask.tif\n", + "Chipping: 628a 628a.MASK.tif 628a.HLS_S30.stacked.tif 628a.HLS_S30.2019-07-13.Fmask.tif\n", + "Chipping: 638a 638a.MASK.tif 638a.HLS_S30.stacked.tif 638a.HLS_S30.2020-06-12.Fmask.tif\n", + "Chipping: 116d 116d.MASK.tif 116d.HLS_S30.stacked.tif 116d.HLS_S30.2020-06-17.Fmask.tif\n", + "Chipping: 116e 116e.MASK.tif 116e.HLS_S30.stacked.tif 116e.HLS_S30.2020-06-17.Fmask.tif\n", + "Chipping: 648a 648a.MASK.tif 648a.HLS_S30.stacked.tif 648a.HLS_S30.2020-07-16.Fmask.tif\n", + "Chipping: 1055a 1055a.MASK.tif 1055a.HLS_S30.stacked.tif 1055a.HLS_S30.2018-08-04.Fmask.tif\n", + "Chipping: 1056a 1056a.MASK.tif 1056a.HLS_S30.stacked.tif 1056a.HLS_S30.2018-08-07.Fmask.tif\n", + "Chipping: 1347a 1347a.MASK.tif 1347a.HLS_S30.stacked.tif 1347a.HLS_S30.2018-08-04.Fmask.tif\n", + "Chipping: 1064a 1064a.MASK.tif 1064a.HLS_S30.stacked.tif 1064a.HLS_S30.2018-07-03.Fmask.tif\n", + "Chipping: 1338a 1338a.MASK.tif 1338a.HLS_S30.stacked.tif 1338a.HLS_S30.2019-09-05.Fmask.tif\n", + "Chipping: 1373a 1373a.MASK.tif 1373a.HLS_S30.stacked.tif 1373a.HLS_S30.2020-08-17.Fmask.tif\n", + "Chipping: 1373d 1373d.MASK.tif 1373d.HLS_S30.stacked.tif 1373d.HLS_S30.2020-08-17.Fmask.tif\n", + "Chipping: 1373e 1373e.MASK.tif 1373e.HLS_S30.stacked.tif 1373e.HLS_S30.2020-08-17.Fmask.tif\n", + "Chipping: 110d 110d.MASK.tif 110d.HLS_S30.stacked.tif 110d.HLS_S30.2019-09-05.Fmask.tif\n", + "Chipping: 106a 106a.MASK.tif 106a.HLS_L30.stacked.tif 106a.HLS_L30.2019-08-28.fmask.tif\n", + "Chipping: 129a 129a.MASK.tif 129a.HLS_L30.stacked.tif 129a.HLS_L30.2018-08-07.fmask.tif\n", + "Chipping: 130a 130a.MASK.tif 130a.HLS_L30.stacked.tif 130a.HLS_L30.2018-08-07.fmask.tif\n", + "Chipping: 132a 132a.MASK.tif 132a.HLS_L30.stacked.tif 132a.HLS_L30.2018-08-11.fmask.tif\n", + "Chipping: 133a 133a.MASK.tif 133a.HLS_L30.stacked.tif 133a.HLS_L30.2018-07-26.fmask.tif\n", + "Chipping: 614a 614a.MASK.tif 614a.HLS_L30.stacked.tif 614a.HLS_L30.2018-07-08.fmask.tif\n", + "Chipping: 623d 623d.MASK.tif 623d.HLS_L30.stacked.tif 623d.HLS_L30.2019-08-19.fmask.tif\n", + "Chipping: 623a 623a.MASK.tif 623a.HLS_L30.stacked.tif 623a.HLS_L30.2019-08-19.fmask.tif\n", + "Chipping: 915a 915a.MASK.tif 915a.HLS_L30.stacked.tif 915a.HLS_L30.2017-07-26.fmask.tif\n", + "Chipping: 1378a 1378a.MASK.tif 1378a.HLS_L30.stacked.tif 1378a.HLS_L30.2018-07-12.fmask.tif\n", + "Chipping: 1052a 1052a.MASK.tif 1052a.HLS_L30.stacked.tif 1052a.HLS_L30.2017-08-13.fmask.tif\n", + "Chipping: 1376a 1376a.MASK.tif 1376a.HLS_L30.stacked.tif 1376a.HLS_L30.2020-07-19.fmask.tif\n" + ] + } + ], + "source": [ + "chipped_events = set() \n", + "event_chip_stacks = []\n", + "\n", + "for evt, modality, (m, b, f) in event_files:\n", + " if evt.name not in keepers[modality['id']]:\n", + " continue\n", + "\n", + " print('Chipping: ', evt.name, m.name, b.name, f.name)\n", + " if evt.name in chipped_events:\n", + " print(f'Already chipped event {evt.name}: skipping...')\n", + " continue\n", + "\n", + " chipped_events.add(evt.name)\n", + " grid = chip_data.make_grid_from_reference(m)\n", + "\n", + " label_chips = chip_data.chip_data(grid, m, chip_paths['all']['label'])\n", + " data_chips = chip_data.chip_data(grid, b, chip_paths['all']['hls'])\n", + " fmask_chips = chip_data.chip_data(grid, f, chip_paths['all']['other'])\n", + "\n", + " event_chips = chip_data.make_chip_stacks(data_chips, fmask_chips, label_chips, modality)\n", + " event_chip_stacks.append((event_chips, evt))" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "a022ca2c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "([ChipStack(id='000.000', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.000.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.000.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.000.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.001', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.001.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.001.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.001.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.002', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.002.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.002.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.002.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.003', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.003.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.003.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.003.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.004', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.004.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.004.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.004.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.005', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.005.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.005.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.005.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.006', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.006.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.006.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.006.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.007', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.007.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.007.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.007.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.008', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.008.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.008.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.008.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='000.009', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/000.009.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/000.009.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/000.009.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.000', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.000.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.000.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.000.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.001', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.001.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.001.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.001.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.002', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.002.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.002.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.002.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.003', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.003.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.003.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.003.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.004', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.004.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.004.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.004.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.005', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.005.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.005.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.005.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.006', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.006.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.006.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.006.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.007', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.007.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.007.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.007.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.008', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.008.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.008.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.008.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='001.009', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/001.009.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/001.009.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/001.009.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.000', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.000.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.000.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.000.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.001', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.001.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.001.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.001.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.002', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.002.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.002.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.002.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.003', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.003.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.003.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.003.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.004', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.004.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.004.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.004.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.005', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.005.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.005.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.005.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.006', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.006.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.006.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.006.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.007', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.007.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.007.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.007.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.008', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.008.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.008.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.008.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))}), ChipStack(id='002.009', data=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/HLS/002.009.102a.HLS_S30.stacked.tif'), validation_mask=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/OTHER/002.009.102a.HLS_S30.2019-07-17.Fmask.tif'), label=PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/chips/LABEL/002.009.102a.MASK.tif'), modality={'id': 'HLS_S30', 'collection': 'HLSS30', 'bands': (Band(id='B02', name='Blue', shortname='B'), Band(id='B03', name='Green', shortname='G'), Band(id='B04', name='Red', shortname='R'), Band(id='B8A', name='NIR Narrow', shortname='N'), Band(id='B11', name='SWIR 1', shortname='SW1'), Band(id='B12', name='SWIR 2', shortname='SW2'), Band(id='Fmask', name='Cloud Mask', shortname='Fmask'))})], Event(name='102a', date=Timestamp('2019-07-17 00:00:00'), wgs84_geometry=, buffer_m=10000))\n" + ] + } + ], + "source": [ + "print(event_chip_stacks[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "637de3af", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "filtering chips for 102a\n", + " damage: 10, no damage: 15\n", + "filtering chips for 126a\n", + " damage: 3, no damage: 4\n", + "filtering chips for 127a\n", + " damage: 16, no damage: 24\n", + "filtering chips for 134a\n", + " damage: 8, no damage: 12\n", + "filtering chips for 598a\n", + " damage: 7, no damage: 10\n", + "filtering chips for 889a\n", + " damage: 10, no damage: 15\n", + "filtering chips for 889b\n", + " damage: 5, no damage: 7\n", + "filtering chips for 892a\n", + " damage: 15, no damage: 22\n", + "filtering chips for 1079a\n", + " damage: 25, no damage: 37\n", + "filtering chips for 1069a\n", + " damage: 3, no damage: 4\n", + "filtering chips for 628a\n", + " damage: 1, no damage: 1\n", + "filtering chips for 638a\n", + " damage: 47, no damage: 70\n", + "filtering chips for 116d\n", + " damage: 6, no damage: 9\n", + "filtering chips for 116e\n", + " damage: 3, no damage: 4\n", + "filtering chips for 648a\n", + " damage: 4, no damage: 6\n", + "filtering chips for 1055a\n", + " damage: 24, no damage: 36\n", + "filtering chips for 1056a\n", + " damage: 22, no damage: 33\n", + "filtering chips for 1347a\n", + " damage: 18, no damage: 27\n", + "filtering chips for 1064a\n", + " damage: 13, no damage: 19\n", + "filtering chips for 1338a\n", + " damage: 7, no damage: 10\n", + "filtering chips for 1373a\n", + " damage: 3, no damage: 4\n", + "filtering chips for 1373d\n", + " damage: 4, no damage: 6\n", + "filtering chips for 1373e\n", + " damage: 7, no damage: 10\n", + "filtering chips for 110d\n", + " damage: 7, no damage: 10\n", + "filtering chips for 106a\n", + " damage: 23, no damage: 33\n", + "filtering chips for 129a\n", + " damage: 13, no damage: 19\n", + "filtering chips for 130a\n", + " damage: 17, no damage: 25\n", + "filtering chips for 132a\n", + " damage: 17, no damage: 25\n", + "filtering chips for 133a\n", + " damage: 3, no damage: 3\n", + "filtering chips for 614a\n", + " damage: 56, no damage: 84\n", + "filtering chips for 623d\n", + " damage: 8, no damage: 11\n", + "filtering chips for 623a\n", + " damage: 5, no damage: 7\n", + "filtering chips for 915a\n", + " damage: 7, no damage: 10\n", + "filtering chips for 1378a\n", + " damage: 4, no damage: 6\n", + "filtering chips for 1052a\n", + " damage: 4, no damage: 4\n", + "filtering chips for 1376a\n", + " damage: 1, no damage: 1\n" + ] + } + ], + "source": [ + "import random\n", + "\n", + "all_chips = []\n", + "filtered_chips = []\n", + "all_no_damage = []\n", + "all_damage = []\n", + "\n", + "for event_chips, event in event_chip_stacks:\n", + " print(f'filtering chips for {event.name}')\n", + " all_chips += event_chips\n", + " view.view_chips(\n", + " event_chips, \n", + " modality, \n", + " rgb_bands=[2, 1, 0], \n", + " save_to_file=chip_paths['all']['plots'] / f'{event.name}.chips.png', \n", + " quite=True\n", + " )\n", + " \n", + " damage_chips = chip_data.filter_damage_chips(event_chips)\n", + " all_damage += damage_chips\n", + "\n", + " no_damage_chips = chip_data.filter_no_damage_chips(event_chips)\n", + "\n", + " random.shuffle(no_damage_chips)\n", + " no_damage_chips = no_damage_chips[: int(len(damage_chips) * 1.5)]\n", + "\n", + " print(f' damage: {len(damage_chips)}, no damage: {len(no_damage_chips)}')\n", + "\n", + " if len(damage_chips) == 0 and len(no_damage_chips) == 0:\n", + " print(f'Filtered out all chips for {event.name}')\n", + " continue\n", + "\n", + " all_no_damage += no_damage_chips\n", + " filtered = damage_chips + no_damage_chips\n", + " \n", + " view.view_chips(\n", + " filtered, \n", + " modality, \n", + " rgb_bands=[2, 1, 0], \n", + " save_to_file=chip_paths['output']['plots'] / f'{event.name}.filtered.png', \n", + " quite=True\n", + " )\n", + " \n", + " filtered_chips += filtered" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "252dff47", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(426, 623, 1049, 1688)" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(all_damage), len(all_no_damage), len(filtered_chips), len(all_chips)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "7d4aeeb6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/output/HLS')" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "chip_paths['output']['hls']" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "6a7dca8f", + "metadata": {}, + "outputs": [], + "source": [ + "chip_paths['output']['hls'].mkdir(exist_ok=True, parents=True)\n", + "chip_paths['output']['label'].mkdir(exist_ok=True, parents=True)\n", + "\n", + "for chip in filtered_chips:\n", + " output_data_path = chip_paths['output']['hls'] / chip.data.name \n", + " output_label_path = chip_paths['output']['label'] / chip.label.name\n", + "\n", + " shutil.copy2(chip.data, output_data_path)\n", + " shutil.copy2(chip.label, output_label_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "b1c32b82", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1049" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(list(chip_paths['output']['hls'].glob('*.tif')))\n", + "len(list(chip_paths['output']['label'].glob('*.tif')))" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "6f82016d", + "metadata": {}, + "outputs": [], + "source": [ + "#len(filtered_chips)\n", + "#filtered_chips[0]\n", + "#CHIPS_PLOTS = PLOT_PATH / 'CHIPS'\n", + "#print(len(filtered_chips), len(all_chips))\n", + "#for chip in filtered_chips:\n", + "# plot_name = f\"{chip.data.name.split('.stacked')[0]}.chip.png\"\n", + "# view.view_chip(chip, models.HLS_S30, [2, 1, 0], save_to_file=CHIPS_PLOTS / plot_name, quite=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "3049d797", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import rasterio\n", + "\n", + "\n", + "def calculate_stats(chip_stacks, n_bands) -> tuple:\n", + " mean = np.zeros(n_bands, dtype=np.float64)\n", + " M2 = np.zeros(n_bands, dtype=np.float64)\n", + " count = np.zeros(n_bands, dtype=np.float64)\n", + "\n", + " for chip in chip_stacks:\n", + " with rasterio.open(chip.data) as src:\n", + " band_data = src.read()\n", + " count, mean, M2 = 0, 0, 0\n", + "\n", + " _, H, W = band_data.shape\n", + "\n", + " batch_count = H * W\n", + " batch_mean = band_data.mean(axis=(1, 2))\n", + " batch_var = band_data.var(axis=(1, 2))\n", + "\n", + " delta = batch_mean - mean\n", + " total_count = count + batch_count\n", + "\n", + " mean = mean + delta * (batch_count / total_count)\n", + " M2 = (\n", + " M2\n", + " + batch_var * batch_count\n", + " + (delta**2) * count * batch_count / total_count\n", + " )\n", + " count = total_count\n", + "\n", + " variance = M2 / count\n", + " std = np.sqrt(variance)\n", + "\n", + " return mean, std\n" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "00880cff", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged/1376a.HLS_L30.2020-07-19.B.tif'),\n", + " PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged/1376a.HLS_L30.2020-07-19.G.tif'),\n", + " PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged/1376a.HLS_L30.2020-07-19.R.tif'),\n", + " PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged/1376a.HLS_L30.2020-07-19.N.tif'),\n", + " PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged/1376a.HLS_L30.2020-07-19.SW1.tif'),\n", + " PosixPath('/home/wbhorn/Repositories/fm/satchip/notebooks/data/HLS_L30/merged/1376a.HLS_L30.2020-07-19.SW2.tif')]" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data_bands" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "f68fa713", + "metadata": {}, + "outputs": [], + "source": [ + "chip_means, chip_stds = calculate_stats(filtered_chips, 6)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "7a6d8b65", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'HLS': {'means': {'B': 267.3031005859375,\n", + " 'G': 500.4424133300781,\n", + " 'R': 332.49700927734375,\n", + " 'N': 4312.5751953125,\n", + " 'SW1': 1867.71044921875,\n", + " 'SW2': 887.797607421875},\n", + " 'stds': {'B': 112.45736694335938,\n", + " 'G': 160.5001678466797,\n", + " 'R': 199.83193969726562,\n", + " 'N': 594.2197265625,\n", + " 'SW1': 496.0038757324219,\n", + " 'SW2': 382.4466857910156}}}" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import json\n", + "\n", + "\n", + "bands = models.HLS_S30['bands'][:-1]\n", + "bands\n", + "\n", + "stats_file = {\n", + " 'HLS':{\n", + " 'means': {},\n", + " 'stds': {},\n", + " }\n", + "}\n", + "\n", + "for mean, std, band in zip(chip_means, chip_stds, bands):\n", + " stats_file['HLS']['means'][band.shortname] = float(mean)\n", + " stats_file['HLS']['stds'][band.shortname] = float(std)\n", + "\n", + "stats_file\n", + "\n", + "(output_base / 'statistics.json').write_text(json.dumps(stats_file, indent=2))\n", + "stats_file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "312e53f0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "satchip", + "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.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/satchip/chip_data.py b/src/satchip/chip_data.py index 847e82d..3d43657 100644 --- a/src/satchip/chip_data.py +++ b/src/satchip/chip_data.py @@ -1,146 +1,190 @@ -import argparse -from collections import Counter -from datetime import datetime from pathlib import Path import numpy as np -import xarray as xr -from shapely.geometry import box -from tqdm import tqdm - -from satchip import utils -from satchip.chip_hls import get_hls_data -from satchip.chip_hyp3s1rtc import get_rtc_paths_for_chips, get_s1rtc_chip_data -from satchip.chip_operas1rtc import get_operartc_data -from satchip.chip_sentinel2 import get_s2l2a_data -from satchip.terra_mind_grid import TerraMindChip, TerraMindGrid - - -def fill_missing_times(data_chip: xr.DataArray, times: np.ndarray) -> xr.DataArray: - missing_times = np.setdiff1d(times, data_chip.time.data) - missing_shape = (len(missing_times), len(data_chip.band), data_chip.y.size, data_chip.x.size) - missing_data = xr.DataArray( - np.full(missing_shape, 0, dtype=data_chip.dtype), - dims=('time', 'band', 'y', 'x'), - coords={ - 'time': missing_times, - 'band': data_chip.band.data, - 'y': data_chip.y.data, - 'x': data_chip.x.data, - }, - ) - return xr.concat([data_chip, missing_data], dim='time').sortby('time') +import rasterio +from rasterio.windows import Window + +from satchip import models + + +def make_grid_from_reference(reference: Path, chip_size: int = 256) -> list[models.GridCell]: + grid = [] + + with rasterio.open(reference) as ref: + n_cols = ref.width // chip_size + n_rows = ref.height // chip_size + + for row in range(n_rows): + for col in range(n_cols): + window = Window(col * chip_size, row * chip_size, chip_size, chip_size) + bounds = ref.window_bounds(window) + + cell_id = f'{row:03d}.{col:03d}' + grid.append(models.GridCell(cell_id, bounds)) -def get_chips(label_paths: list[Path]) -> list[TerraMindChip]: - label_datasets = [utils.load_chip(label_path) for label_path in label_paths] - bounds = utils.get_overall_bounds([ds.bounds for ds in label_datasets]) + return grid - buffered = box(*bounds).buffer(0.5).bounds - grid = TerraMindGrid(latitude_range=(buffered[1], buffered[3]), longitude_range=(buffered[0], buffered[2])) - grid_chips = {chip.name: chip for chip in grid.terra_mind_chips} + +def chip_data(grid: list[models.GridCell], layer: Path, output_path: Path) -> list[models.Chips]: + output_path.mkdir(exist_ok=True, parents=True) chips = [] - for label_dataset in label_datasets: - label_chip_name = label_dataset.sample.item() - assert label_chip_name in grid_chips, f'No TerraMind chip found for label {label_chip_name}' - chip = grid_chips[label_chip_name] - chips.append(chip) + + with rasterio.open(layer) as src: + for grid_cell in grid: + window = src.window(*grid_cell.bounds) + window = Window( + round(window.col_off), + round(window.row_off), + round(window.width), + round(window.height), + ) + + data = src.read(window=window) + + chip_meta = src.meta.copy() + chip_meta.update( + { + 'width': window.width, + 'height': window.height, + 'transform': src.window_transform(window), + } + ) + + chip_name = f'{grid_cell.id}.{layer.name}' + chip_path = output_path / chip_name + + with rasterio.open(chip_path, 'w', **chip_meta) as dst: + dst.write(data) + + chip = models.Chip(grid_cell.id, chip_path) + chips.append(chip) return chips -def chip_data( - chip: TerraMindChip, - platform: str, - opts: utils.ChipDataOpts, - image_dir: Path, -) -> xr.Dataset: - if platform == 'HYP3S1RTC': - rtc_paths = opts['local_hyp3_paths'][chip.name] - chip_dataset = get_s1rtc_chip_data(chip, rtc_paths) - elif platform == 'S1RTC': - chip_dataset = get_operartc_data(chip, image_dir, opts=opts) - elif platform == 'S2L2A': - chip_dataset = get_s2l2a_data(chip, image_dir, opts=opts) - elif platform == 'HLS': - chip_dataset = get_hls_data(chip, image_dir, opts=opts) - else: - raise Exception(f'Unknown platform {platform}') - - return chip_dataset - - -def create_chips( - label_paths: list[Path], - platform: str, - date_start: datetime, - date_end: datetime, - strategy: str, - max_cloud_pct: int, - chip_dir: Path, - image_dir: Path, -) -> list[Path]: - platform_dir = chip_dir / platform - platform_dir.mkdir(parents=True, exist_ok=True) - - opts: utils.ChipDataOpts = {'strategy': strategy, 'date_start': date_start, 'date_end': date_end} - if platform in ['S2L2A', 'HLS']: - opts['max_cloud_pct'] = max_cloud_pct - - chips = get_chips(label_paths) - chip_names = [c.name for c in chips] - if len(chip_names) != len(set(chip_names)): - duplicates = [name for name, count in Counter(chip_names).items() if count > 1] - msg = f'Duplicate sample locations not supported. Duplicate chips: {", ".join(duplicates)}' - raise NotImplementedError(msg) - chip_paths = [ - platform_dir / (x.with_suffix('').with_suffix('').name + f'_{platform}.zarr.zip') for x in label_paths - ] - if platform == 'HYP3S1RTC': - rtc_paths_for_chips = get_rtc_paths_for_chips(chips, image_dir, opts) - opts['local_hyp3_paths'] = rtc_paths_for_chips - - for chip, chip_path in tqdm(list(zip(chips, chip_paths)), desc='Chipping labels'): - dataset = chip_data(chip, platform, opts, image_dir) - utils.save_chip(dataset, chip_path) - return chip_paths - - -def main() -> None: - parser = argparse.ArgumentParser(description='Chip a label image') - parser.add_argument('labelpath', type=Path, help='Path to the label directory') - parser.add_argument( - 'platform', choices=['S1RTC', 'S2L2A', 'HLS', 'HYP3S1RTC'], type=str, help='Dataset to create chips for' - ) - parser.add_argument('daterange', type=str, help='Inclusive date range to search for data in the format Ymd-Ymd') - parser.add_argument('--maxcloudpct', default=100, type=int, help='Maximum percent cloud cover for a data chip') - parser.add_argument('--chipdir', default='.', type=Path, help='Output directory for the chips') - parser.add_argument( - '--imagedir', default=None, type=Path, help='Output directory for image files. Defaults to chipdir/IMAGES' - ) - parser.add_argument( - '--strategy', - default='BEST', - choices=['BEST', 'ALL'], - type=str, - help='Strategy to use when multiple scenes are found (default: BEST)', - ) - args = parser.parse_args() - args.platform = args.platform.upper() - assert 0 <= args.maxcloudpct <= 100, 'maxcloudpct must be between 0 and 100' - date_start, date_end = [datetime.strptime(d, '%Y%m%d') for d in args.daterange.split('-')] - assert date_start < date_end, 'start date must be before end date' - label_paths = list(args.labelpath.glob('*.zarr.zip')) - assert len(label_paths) > 0, f'No label files found in {args.labelpath}' - - if args.imagedir is None: - args.imagedir = args.chipdir / 'IMAGES' - - create_chips( - label_paths, args.platform, date_start, date_end, args.strategy, args.maxcloudpct, args.chipdir, args.imagedir +def make_chip_stacks( + data_chips: list[models.Chip], + validation_mask_chips: list[models.Chip], + label_chips: list[models.Chip], + modality: models.Modality, +) -> list[models.ChipStack]: + chip_stacks = [] + + for data, mask, label in zip( + sorted(data_chips, key=lambda c: c.id), + sorted(validation_mask_chips, key=lambda c: c.id), + sorted(label_chips, key=lambda c: c.id), + ): + chip_stack = models.ChipStack( + id=data.id, + data=data.path, + validation_mask=mask.path, + label=label.path, + modality=modality, + ) + + chip_stacks.append(chip_stack) + + return chip_stacks + + +def filter_damage_chips(chip_stacks: list[models.ChipStack]) -> list[models.ChipStack]: + good_chips = [] + + for chip_stack in chip_stacks: + if 'HLS' in chip_stack.modality['id']: + is_good_chip = is_good_damage_hls_chip(chip_stack) + elif 'RTC' in chip_stack.modality['id']: + is_good_chip = is_good_damage_rtc_chip(chip_stack) + + if is_good_chip: + good_chips.append(chip_stack) + + return good_chips + + +def filter_no_damage_chips(chip_stacks: list[models.ChipStack]) -> list[models.ChipStack]: + good_chips = [] + + for chip_stack in chip_stacks: + if 'HLS' in chip_stack.modality['id']: + is_good_chip = is_good_no_damage_hls_chip(chip_stack) + elif 'RTC' in chip_stack.modality['id']: + is_good_chip = is_good_no_damage_rtc_chip(chip_stack) + + if is_good_chip: + good_chips.append(chip_stack) + + return good_chips + + +def _hls_chip_stats(chip_stack: models.ChipStack) -> tuple[float, float]: + with rasterio.open(chip_stack.validation_mask) as ds: + qc = clear_px_Fmask(ds.read(1)) + with rasterio.open(chip_stack.label) as ds: + event = ds.read(1) + + ny, nx = qc.shape + n_px = 1.0 * ny * nx + + # cloud-free pixels (0 clear, 1 cloud, 255 nodata) + n_cf = len(np.where(qc == 0)[0]) + pct_cf = 100.0 * (n_cf / n_px) + + # pct of chip in event + n_ev = len(np.where(event > 0)[0]) + pct_ev = 100.0 * (n_ev / n_px) if n_ev > 0 else 0 + + return pct_cf, pct_ev + + +def is_good_damage_hls_chip(chip_stack: models.ChipStack) -> bool: + pct_cf, pct_ev = _hls_chip_stats(chip_stack) + return pct_cf > 95 and pct_ev > 1 + + +def is_good_no_damage_hls_chip(chip_stack: models.ChipStack) -> bool: + pct_cf, pct_ev = _hls_chip_stats(chip_stack) + return pct_cf > 95 and pct_ev < 1 + + +def clear_px_Fmask(Fmask: np.ndarray) -> np.ndarray: + fmask_clear = np.array( + [0, 4, 16, 20, 32, 36, 48, 52, 64, 68, 80, 84, 96, 100, 112, 116, + 128, 132, 144, 148, 160, 164, 176, 180, 192, 196, 208, 212, 224, 228, 240, 244], + dtype=Fmask.dtype, ) + # 0 clear, 1 cloud, 255 nodata + cloudmask = np.ones_like(Fmask, dtype=np.uint8) + cloudmask[np.isin(Fmask, fmask_clear)] = 0 + cloudmask[Fmask == 255] = 255 + return cloudmask + + +def _rtc_chip_stats(chip_stack: models.ChipStack) -> tuple[bool, bool]: + with rasterio.open(chip_stack.data) as ds: + rtc_data = ds.read() + with rasterio.open(chip_stack.label) as ds: + event_mask = ds.read(1) + + has_nan_pixels = np.isnan(rtc_data).sum() > 0 + + # pct of chip in event + num_pixels = event_mask.size + num_event_pixels = np.count_nonzero(event_mask > 0) + pct_pixels_over_event = 100.0 * (num_event_pixels / num_pixels) + data_overlaps_event = pct_pixels_over_event > 1 + + return has_nan_pixels, data_overlaps_event + + +def is_good_damage_rtc_chip(chip_stack: models.ChipStack) -> bool: + has_nan_pixels, data_overlaps_event = _rtc_chip_stats(chip_stack) + return not has_nan_pixels and data_overlaps_event -if __name__ == '__main__': - main() +def is_good_no_damage_rtc_chip(chip_stack: models.ChipStack) -> bool: + has_nan_pixels, data_overlaps_event = _rtc_chip_stats(chip_stack) + return not has_nan_pixels and not data_overlaps_event diff --git a/src/satchip/download_data.py b/src/satchip/download_data.py new file mode 100644 index 0000000..73930d8 --- /dev/null +++ b/src/satchip/download_data.py @@ -0,0 +1,36 @@ +from datetime import timedelta +from pathlib import Path + +import earthaccess + +from satchip import models + + +def download_data(event: models.Event, modality: models.Modality, download_path: Path) -> list[Path]: + if not earthaccess.__auth__.authenticated: + print('Logging in to earthaccess') + earthaccess.login() + + results = _search_data(event, modality) + + if not results: + return [] + + local_files = earthaccess.download(results, local_path=download_path, show_progress=True) + + return local_files + + +def _search_data(event: models.Event, modality: models.Modality) -> list[earthaccess.DataGranule]: + start_date = event.date + final_date = start_date + timedelta(days=1) + + collection_id = modality['collection'] + + results = earthaccess.search_data( + short_name=[collection_id], + temporal=(start_date.strftime('%Y-%m-%d'), final_date.strftime('%Y-%m-%d')), + bounding_box=event.buffered_geometry().bounds, + ) + + return results diff --git a/src/satchip/generate_chips.py b/src/satchip/generate_chips.py new file mode 100644 index 0000000..f432766 --- /dev/null +++ b/src/satchip/generate_chips.py @@ -0,0 +1,433 @@ +import shutil +import zipfile +from pathlib import Path + +import cartopy.crs as ccrs +import earthaccess +import gdown +import geopandas as gpd +import hls +import matplotlib.pyplot as plt +import mosaic +import numpy as np +import opera_rtc +import pandas as pd +import rasterio +from modality import Modality +from rasterio import features +from rasterio.windows import Window +from shapely.geometry import box +from sklearn.model_selection import train_test_split + + +QUITE = True +CHIP_SIZE = 256 +RNG_SEED = 42 + + +MODALITY = 'RTC' +ALL_BANDS = ('VV', 'VH', 'mask') +STACK_BANDS = ('VV', 'VH') +CHIP_BANDS = ('BANDS', 'EVENT', 'MASK') + +MODALITY = 'HLS' +ALL_BANDS = ('B', 'G', 'R', 'N', 'SW1', 'SW2', 'Fmask') +STACK_BANDS = ('B', 'G', 'R', 'N', 'SW1', 'SW2') +CHIP_BANDS = ('BANDS', 'EVENT', 'MASK', 'Fmask') +SHOULD_CLEANUP = False + + +def main(modalities: list[Modality]): + hwds_path = Path('hwds') + + print('Making folders') + data_paths = { + 'CHIPS_ALL': hwds_path / 'CHIPS_ALL', + 'CHIPS': hwds_path / 'CHIPS', + 'MERGED': hwds_path / 'MERGED', + 'PLOTS': hwds_path / 'PLOTS', + } + + if SHOULD_CLEANUP: + for item in data_paths['MERGED'].glob('*.tif'): + if item.is_dir(): + continue + + item.unlink() + + for p in ('CHIPS', 'CHIPS_ALL'): + shutil.rmtree(data_paths[p], ignore_errors=True) + + for p in data_paths.values(): + p.mkdir(parents=True, exist_ok=True) + + gdf = _load_event_database(hwds_path) + gdf_utm = gdf.to_crs(32615) + gdf['buffered_event'] = gdf_utm.buffer(3000).to_crs(4326) + gdf['buffered_event_background'] = gdf_utm.buffer(10000).to_crs(4326) + gdf = gdf.to_crs(4326) + + # keepers = [1442, 622, 1079, 628] + keepers = [1442, 622] + gdf = gdf[gdf['swathID'].isin(keepers)] + + earthaccess.login() + + tm_chips = [] + + for i, (swathID, swath) in enumerate(gdf.iterrows(), start=1): + swathID = f'{int(swath["swathID"]):04d}' + print(f'Processing Swath {swathID} ({i} / {len(gdf)})') + + merged = mosaic.data_over_swath(swath, modalities, output_path=data_paths['MERGED']) + + template_path = merged[modalities[0].id]['BANDS'] + event_tif, mask_tif = _generate_masks(template_path, swathID, swath) + + merged_data = { + **merged, + 'EVENT': event_tif, + 'MASK': mask_tif, + } + + if not all(is_valid_data(merged_data, modality) for modality in modalities): + print('Skipping: not enough valid data') + continue + + for modality in modalities: + print(f'Chipping {modality.id}!') + chips = _chip_data(merged_data, data_paths['CHIPS_ALL'], modality) + + good_chips = filter_chips(chips, modality) + print(f'Found {len(good_chips)} good chips') + + for chip in good_chips: + for band, chip_path in chip.items(): + if band not in ('MASK', 'BANDS'): + continue + + dest = data_paths['CHIPS'] / chip_path.name + shutil.copy(chip_path, dest) + + tm_chips += good_chips + + for _, swath in gdf.iterrows(): + swath_id = _make_swath_id(swath['swathID']) + + for modality in modalities: + merged_file = list(data_paths['MERGED'].glob(f'{swath_id}.{modality.id}.*.BANDS.tif')) + + if len(merged_file) == 0: + print(f'no chips for {swath_id}') + continue + + all_chips = list(data_paths['CHIPS_ALL'].glob(f'*.{swath_id}.{modality.id}.*.tif')) + good_chips = list(data_paths['CHIPS'].glob(f'*.{swath_id}.{modality.id}.*.tif')) + + print(f'plotting {swath_id}') + _plot_chips(merged_file[0], all_chips, good_chips, swath, modality, save_to=data_paths['PLOTS']) + + for modality in modalities: + print(f'Calulating stats for modality: {modality.id}') + band_chips = list(data_paths['CHIPS'].glob(f'*.{modality.id}.*.BANDS.tif')) + + means, stds = calculate_stats(chips=band_chips, n_bands=len(modality.stack_bands)) + + means_str = ', '.join(f'{x:.4f}' for x in means) + stds_str = ', '.join(f'{x:.4f}' for x in stds) + + stats_str = f'Means {modality.stack_bands}: {means_str}\nStds {modality.stack_bands}: {stds_str}\n' + (hwds_path / '{modality.id}-statistics.txt').write_text(stats_str) + print(stats_str) + + +def _generate_masks(template_data_path: Path, swathID: str, swath: pd.Series) -> tuple[Path, Path]: + event_path = template_data_path.parent / f'{swathID}.EVENT.tif' + mask_path = template_data_path.parent / f'{swathID}.MASK.tif' + + with rasterio.open(template_data_path) as ds: + profile = ds.profile + + mask_raster = features.rasterize( + shapes=[[swath['geometry'], 1]], + fill=0, + out_shape=ds.shape, + transform=ds.transform, + ) + + with rasterio.open(mask_path, 'w', **profile) as dst: + dst.write(mask_raster, 1) + print('generated:', mask_path) + + event_mask = features.rasterize( + shapes=[ + [swath['buffered_event_background'], 3], + [swath['buffered_event'], 2], + [swath['geometry'], 1], + ], + fill=0, + out_shape=ds.shape, + transform=ds.transform, + ) + + with rasterio.open(event_path, 'w', **profile) as dst: + dst.write(event_mask, 1) + print('generated:', event_path) + + return event_path, mask_path + + +def _load_event_database(data_dir: Path): + # use 60-swath version + hwds_google_drive_id = '1h_JIEcrrUF3OSTrmwAKNPa0eUEhPA2Xx' + drive_url = f'https://drive.google.com/uc?id={hwds_google_drive_id}' + + shp_dir = data_dir / 'SHP' + shp_dir.mkdir(parents=True, exist_ok=True) + + filename = 'hwds_v3_20250205_subset_60.zip' + + zip_path = shp_dir / filename + + if not zip_path.exists(): + gdown.download(drive_url, str(zip_path), quiet=False) + + with zipfile.ZipFile(zip_path, 'r') as zip_ref: + zip_ref.extractall(path=shp_dir) + + shp_path = shp_dir / 'hwds_v3_20250205_subset_60.shp' + gdf = gpd.read_file(shp_path) + + gdf['swathDate'] = pd.to_datetime(gdf['swathDate'], format='%Y-%m-%d') + gdf['ls5hlsDate'] = pd.to_datetime(gdf['ls5hlsDate'], format='%Y-%m-%d') + gdf['s1Date'] = pd.to_datetime(gdf['s1Date'], format='%Y-%m-%d') + + return gdf + + +def create_split_files(band_chips: list[Path], splits_path: Path) -> None: + chip_ids = [p.name.removesuffix('BANDS.tif') for p in band_chips] + + the_rest, test = train_test_split(chip_ids, test_size=0.15, random_state=RNG_SEED) + train, val = train_test_split(the_rest, test_size=0.15, random_state=RNG_SEED) + + splits = {'train': train, 'val': val, 'test': test} + + for split, chip_ids in splits.items(): + split_path = splits_path / f'{split}.txt' + split_path.write_text('\n'.join(chip_ids)) + + +def calculate_stats(chips: list[Path], n_bands: int = 2) -> tuple: + mean = np.zeros(n_bands, dtype=np.float64) + M2 = np.zeros(n_bands, dtype=np.float64) + count = np.zeros(n_bands, dtype=np.float64) + + for chip in chips: + with rasterio.open(chip) as src: + band_data = src.read() + count, mean, M2 = 0, 0, 0 + + _, H, W = band_data.shape + + batch_count = H * W + batch_mean = band_data.mean(axis=(1, 2)) + batch_var = band_data.var(axis=(1, 2)) + + delta = batch_mean - mean + total_count = count + batch_count + + mean = mean + delta * (batch_count / total_count) + M2 = M2 + batch_var * batch_count + (delta**2) * count * batch_count / total_count + count = total_count + + variance = M2 / count + std = np.sqrt(variance) + + return mean, std + + +def _plot_chips( + merged_band_file, + all_chips, + good_chips, + swath, + modality, + save_to: Path | None = None, + quite=QUITE, +): + crs_pc = ccrs.PlateCarree() + + with rasterio.open(merged_band_file) as ds: + bounds = ds.bounds + full_extent = [bounds.left, bounds.right, bounds.bottom, bounds.top] + band_data = ds.read() + + img = get_img(band_data, modality) + + # plot BANDS and geom + fig, ax = plt.subplots( + 1, + 1, + subplot_kw={'projection': crs_pc}, + figsize=(12, 12), + layout='constrained', + ) + + swath_geom = swath['geometry'] + + ax.imshow(img, extent=full_extent, origin='upper', transform=crs_pc) + ax.add_geometries([swath_geom], edgecolor='red', linewidth=2, facecolor='none', crs=crs_pc) + + def show_chips(chips, color, linewidth, z): + for chip in chips: + with rasterio.open(chip) as ds: + chip_bounds = ds.bounds + chip_geom = box( + chip_bounds.left, + chip_bounds.bottom, + chip_bounds.right, + chip_bounds.top, + ) + + ax.add_geometries( + [chip_geom], + edgecolor=color, + linewidth=linewidth, + alpha=1, + zorder=z, + facecolor='none', + crs=crs_pc, + ) + + show_chips(all_chips, 'yellow', 1, z=1) + show_chips(good_chips, 'blue', 3, z=2) + + ax.set_extent(full_extent, crs=crs_pc) + + if save_to: + plt.savefig( + save_to / f'{merged_band_file.name.removesuffix("BANDS.tif")}.png', + dpi=300, + bbox_inches='tight', + ) + + if not quite: + plt.show() + + plt.close(fig) + + +def _make_swath_id(swathID): + return f'{int(swathID):04d}' + + +def _chip_data(merged, output_path: Path, modality: Modality, chip_size=CHIP_SIZE): + chips = {} + grid = [] + + with rasterio.open(merged[modality.id]['BANDS']) as ref: + n_cols = ref.width // chip_size + n_rows = ref.height // chip_size + + for row in range(n_rows): + for col in range(n_cols): + window = Window(col * chip_size, row * chip_size, chip_size, chip_size) + bounds = ref.window_bounds(window) + + tile_id = f'{row:03d}.{col:03d}' + chips[tile_id] = {} + grid.append((tile_id, bounds)) + + for chip_layer in modality.chip_bands: + if chip_layer in merged[modality.id]: + layer_path = merged[modality.id][chip_layer] + else: + layer_path = merged[chip_layer] + + with rasterio.open(layer_path) as src: + for tile_id, bounds in grid: + window = src.window(*bounds) + window = Window( + round(window.col_off), + round(window.row_off), + round(window.width), + round(window.height), + ) + + data = src.read(window=window) + + if chip_layer == 'BANDS': + data = data_transform(data, modality) + + chip_meta = src.meta.copy() + chip_meta.update( + { + 'width': window.width, + 'height': window.height, + 'transform': src.window_transform(window), + } + ) + + chip_name = f'{tile_id}.{layer_path.name}' + chip_path = output_path / chip_name + + with rasterio.open(chip_path, 'w', **chip_meta) as dst: + dst.write(data) + + chips[tile_id][chip_layer] = chip_path + + return chips + + +def is_valid_data(merged, modality): + merged_data = merged[modality.id] + + if modality.id == 'HLS': + is_valid = hls.is_valid_hls(merged_data['Fmask'], merged['EVENT']) + elif modality.id == 'RTC': + is_valid = opera_rtc.is_valid_rtc(merged_data['mask'], merged['EVENT']) + + return is_valid + + +def filter_chips(chips, modality): + if modality.id == 'HLS': + filtered_chips = hls.filter_hls_chips(chips) + elif modality.id == 'RTC': + filtered_chips = opera_rtc.filter_rtc_chips(chips) + + return filtered_chips + + +def get_img(band_data, modality): + if modality.id == 'HLS': + img = hls.get_hls_img(band_data) + elif modality.id == 'RTC': + img = opera_rtc.get_rtc_img(band_data) + + return img + + +def data_transform(data, modality): + if modality.id == 'RTC': + data = 10 * np.log10(np.clip(data, 1e-10, None)) + + return data + + +if __name__ == '__main__': + opera_rtc_mod = Modality( + id='RTC', all_bands=('VV', 'VH', 'mask'), stack_bands=('VV', 'VH'), chip_bands=('BANDS', 'EVENT', 'MASK') + ) + + hls_mod = Modality( + id='HLS', + all_bands=('B', 'G', 'R', 'N', 'SW1', 'SW2', 'Fmask'), + stack_bands=('B', 'G', 'R', 'N', 'SW1', 'SW2'), + chip_bands=('BANDS', 'EVENT', 'MASK', 'Fmask'), + ) + + modalities = [hls_mod, opera_rtc_mod] + + main(modalities) diff --git a/src/satchip/generate_labels.py b/src/satchip/generate_labels.py new file mode 100644 index 0000000..6c23fe8 --- /dev/null +++ b/src/satchip/generate_labels.py @@ -0,0 +1,30 @@ +from pathlib import Path + +import rasterio +from rasterio import features + +from satchip import models + + +def binary_mask_from_template(template_data_path: Path, event: models.Event, output_dir: Path) -> tuple[Path, Path]: + output_dir.mkdir(exist_ok=True, parents=True) + + mask_path = output_dir / f'{event.name}.MASK.tif' + + with rasterio.open(template_data_path) as ds: + profile = ds.profile + + mask_raster = features.rasterize( + shapes=[[event.wgs84_geometry, 1]], + fill=0, + out_shape=ds.shape, + transform=ds.transform, + ) + + profile.update(dtype='uint8', count=1, nodata=255) + + with rasterio.open(mask_path, 'w', **profile) as dst: + dst.write(mask_raster, 1) + print('generated:', mask_path) + + return mask_path diff --git a/src/satchip/hls.py b/src/satchip/hls.py new file mode 100644 index 0000000..7f8478d --- /dev/null +++ b/src/satchip/hls.py @@ -0,0 +1,201 @@ +from datetime import datetime, timedelta +from pathlib import Path + +import earthaccess +import numpy as np +import rasterio +from earthaccess.results import DataGranule + + +def search_hls_data(start_date: datetime, bounding_box: tuple[float, float, float, float]) -> list[DataGranule]: + final_date = start_date + timedelta(days=1) + + # collection_ids = ["C2021957295-LPCLOUD"] # S2 + collection_ids = ['C2021957657-LPCLOUD', 'C2021957295-LPCLOUD'] # S2, L30 + + results = earthaccess.search_data( + concept_id=collection_ids, + temporal=(start_date.strftime('%Y-%m-%d'), final_date.strftime('%Y-%m-%d')), + bounding_box=bounding_box, + cloud_hosted=True, + ) + + return results + + +def band_from_hls_filename(filename): + # HLS.L30.T15TUG.2017167T165321.v2.0.B06.tif + parts = filename.split('.') + + sensor, band = parts[1], parts[-2] + + bands = { + 'L30': { + 'B02': 'B', + 'B03': 'G', + 'B04': 'R', + 'B05': 'N', + 'B06': 'SW1', + 'B07': 'SW2', + 'Fmask': 'Fmask', + }, + 'S30': { + 'B02': 'B', + 'B03': 'G', + 'B04': 'R', + 'B08': 'N', + 'B11': 'SW1', + 'B12': 'SW2', + 'Fmask': 'Fmask', + }, + } + + try: + return bands[sensor][band] + except KeyError: + return '' + + +def make_merged_hls_name(template_filename: str) -> str: + parts = template_filename.split('.') + parts[4] = parts[4][0:7] + parts.pop(3) + parts[-2] = band_from_hls_filename(template_filename) + f_template_merge = '.'.join(parts) + return f_template_merge + + +def clear_px_Fmask(Fmask: np.ndarray) -> np.ndarray: + fmask_clear = np.array( + [ + 0, + 4, + 16, + 20, + 32, + 36, + 48, + 52, + 64, + 68, + 80, + 84, + 96, + 100, + 112, + 116, + 128, + 132, + 144, + 148, + 160, + 164, + 176, + 180, + 192, + 196, + 208, + 212, + 224, + 228, + 240, + 244, + ], + dtype=Fmask.dtype, + ) + + cloudmask = np.ones_like(Fmask, dtype=np.uint8) + cloudmask[np.isin(Fmask, fmask_clear)] = 0 + cloudmask[Fmask == 255] = 255 + + return cloudmask + + +def is_valid_hls(fmask_path: Path, event_path: Path): + with rasterio.open(fmask_path) as ds: + qc = clear_px_Fmask(ds.read(1)) + qc_profile = ds.profile + + with rasterio.open(event_path) as ds: + event_mask = ds.read(1) + event_profile = ds.profile + + print(qc_profile, event_profile) + + ny, nx = np.shape(qc) + mask = np.zeros((ny, nx), 'uint8') + + ok = np.where((event_mask == 2) & (qc == 0)) + n_cf_event = len(ok[0]) + mask[ok] = 1 + + ok = np.where((event_mask == 1) & (qc != 255)) + n_valid_event = len(ok[0]) + + ok = np.where(event_mask == 1) + n_event = len(ok[0]) + + pct_cf_event = 0 + if n_valid_event == 0: + print('No coverage') + else: + pct_cf_event = 100.0 * (n_cf_event / n_event) + print('Percent CF/valid in Event:', pct_cf_event) + + return pct_cf_event > 50 + + +def filter_hls_chips(chips: dict[str, dict]) -> list[dict]: + good_chips = [] + + for tile_id, chip in chips.items(): + with rasterio.open(chip['Fmask']) as ds: + qc = clear_px_Fmask(ds.read(1)) + + with rasterio.open(chip['EVENT']) as ds: + event = ds.read(1) + + ny, nx = qc.shape + n_px = 1.0 * ny * nx + + # cloud-free pixels (0 clear, 1 cloud, 255 nodata) + n_cf = len(np.where(qc == 0)[0]) + pct_cf = 100.0 * (n_cf / n_px) + + # event pixels + n_ev = len(np.where(event > 0)[0]) + + # pct of chip in event + pct_ev = 100.0 * (n_ev / n_px) if n_ev > 0 else 0 + + if pct_cf > 95 and pct_ev > 1: + good_chips.append(chip) + + return good_chips + + +def bytescale(arr, cmin=0, cmax=1, low=0, high=255): + # clip the data to be in the range of cmin to cmax + arr = np.clip(arr, cmin, cmax) + high = float(high) + low = float(low) + cmax = float(cmax) + cmin = float(cmin) + m = (high - low) / (cmax - cmin) # slope + b = high - (m * cmax) # intercept + arr = np.uint8((m * arr) + b) + return arr + + +def get_hls_img(hls_data: np.ndarray) -> np.ndarray: + # B04 + r = bytescale(np.sqrt(np.clip(hls_data[2] / 10000.0, 0, 2)), 0, 0.5) + + # B03 + g = bytescale(np.sqrt(np.clip(hls_data[1] / 10000.0, 0, 2)), 0, 0.5) + + # B02 + b = bytescale(np.sqrt(np.clip(hls_data[0] / 10000.0, 0, 2)), 0, 0.5) + + rgb = np.dstack((r, g, b)) + return rgb diff --git a/src/satchip/merge_modality.py b/src/satchip/merge_modality.py new file mode 100644 index 0000000..cc3679a --- /dev/null +++ b/src/satchip/merge_modality.py @@ -0,0 +1,240 @@ +import datetime +from collections.abc import Iterable +from pathlib import Path + +import numpy as np +import rasterio +from rasterio.crs import CRS +from rasterio.merge import merge +from rasterio.transform import Affine, from_bounds +from rasterio.warp import Resampling, calculate_default_transform, reproject, transform_bounds + +from satchip import models + + +def merge_modality( + modality_files: list[Path], + modality: models.Modality, + event: models.Event, + output_path: Path, + selected_bands: list[models.Band] | None = None, +) -> list[Path]: + reproj_path = output_path / 'wgs84' + + output_path.mkdir(exist_ok=True, parents=True) + reproj_path.mkdir(exist_ok=True, parents=True) + + if len(modality_files) == 0: + print(f'Warning: no data for {event.name}') + return [] + + if selected_bands is None: + selected_bands = modality['bands'] + + merged = [] + + for band in selected_bands: + band_files = [f for f in modality_files if band.id in models.band_id_from_filename(f.name, modality['id'])] + + merged_name = _make_merge_name(event.name, event.date, band.shortname, modality['id']) + + reprojected_files = reproject_files(band_files, reproj_path) + merged_band_path = merge_files(reprojected_files, output_file=output_path / merged_name) + + merged.append(merged_band_path) + + return merged + + +def _make_merge_name(event_name: str, start_date: datetime.datetime, band: str, modality_id: str): + date_str = start_date.date().isoformat() + + return f'{event_name}.{modality_id}.{date_str}.{band}.tif' + + +def merge_files(files: list[Path], output_file: Path) -> Path: + band_datasets = [rasterio.open(band_file) for band_file in files] + + try: + mosaic, out_trans = merge(band_datasets) + mosaic = np.squeeze(mosaic) + + out_meta = band_datasets[0].meta.copy() + + out_meta.update( + { + 'driver': 'GTiff', + 'height': mosaic.shape[0], + 'width': mosaic.shape[1], + 'transform': out_trans, + 'crs': band_datasets[0].crs, + } + ) + + with rasterio.open(output_file, 'w', **out_meta) as dst: + if len(mosaic.shape) == 2: + dst.write(mosaic, 1) + else: + dst.write(mosaic) + + finally: + for ds in band_datasets: + ds.close() + + return output_file + + +def stack_bands(band_files: Iterable[Path], stacked_filename: Path) -> Path: + stacked_filename.parent.mkdir(exist_ok=True, parents=True) + + with rasterio.open(band_files[0]) as src: + meta = src.meta.copy() + + meta.update(count=len(band_files), dtype=np.float32) + + with rasterio.open(stacked_filename, 'w', **meta) as dst: + for idx, band_file in enumerate(band_files, start=1): + with rasterio.open(band_file) as src: + dst.write(src.read(1), idx) + + return stacked_filename + + +def reproject_files(files: list[Path], output_dir: Path) -> list[Path]: + output_dir.mkdir(exist_ok=True, parents=True) + + reprojected_paths = [output_dir / f'{file.name}' for file in files] + + for file, output in zip(files, reprojected_paths): + if output.exists(): + continue + + print(f'reprojecting to wgs84: {output.name}') + reproject_file(file, output) + + return reprojected_paths + + +def reproject_file(local_file: Path, reprojected_file: Path, epsg=4326) -> None: + # https://rasterio.readthedocs.io/en/stable/topics/reproject.html#reprojecting-a-geotiff-dataset + with rasterio.open(local_file) as src: + dst_crs = CRS.from_epsg(epsg) + transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds) + + dst_kwargs = src.meta.copy() + dst_kwargs.update({'crs': dst_crs, 'transform': transform, 'width': width, 'height': height}) + + with rasterio.open(reprojected_file, 'w', **dst_kwargs) as dst: + for i in range(1, src.count + 1): + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs=dst_crs, + ) + + +def warp_to_reference( + reference_path: Path, + data_files: Iterable[Path], + output_dir: Path, + bounding_box_wgs84: tuple(float, float, float, float), +) -> list[Path]: + output_dir.mkdir(parents=True, exist_ok=True) + + dst_transform, width, height, dst_crs = _build_common_grid(bounding_box_wgs84, reference_path) + + files = (reference_path, *data_files) + + output = [] + for data_file in files: + out_path = output_dir / data_file.name + + _warp_single( + data_file, + out_path, + dst_transform, + width, + height, + dst_crs, + ) + + output.append(out_path) + + return output + + +def _warp_single( + input_path: Path, output_path: Path, dst_transform: Affine, width: int, height: int, dst_crs: CRS +) -> Path: + resampling = _get_resampling_method(input_path) + + with rasterio.open(input_path) as src: + dst_data = np.zeros((src.count, height, width), dtype=src.dtypes[0]) + + reproject( + source=rasterio.band(src, list(range(1, src.count + 1))), + destination=dst_data, + src_transform=src.transform, + src_crs=src.crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=resampling, + dst_nodata=src.nodata, + ) + + out_meta = src.meta.copy() + out_meta.update( + { + 'driver': 'GTiff', + 'height': height, + 'width': width, + 'transform': dst_transform, + 'crs': dst_crs, + } + ) + + with rasterio.open(output_path, 'w', **out_meta) as dest: + dest.write(dst_data) + + return output_path + + +def _get_resampling_method(filepath: Path) -> Resampling: + filename = filepath.name.lower() + + if 'fmask' in filename or 'mask' in filename: + return Resampling.nearest + else: + return Resampling.bilinear + + +def _build_common_grid( + bounding_box_4326: tuple(float, float, float, float), reference_path: Path +) -> tuple[Affine, int, int, CRS]: + dst_crs = CRS.from_epsg(4326) + minx, miny, maxx, maxy = bounding_box_4326 + + with rasterio.open(reference_path) as ref: + bounds_4326 = transform_bounds(ref.crs, dst_crs, *ref.bounds, densify_pts=21) + ref_width = ref.width + ref_height = ref.height + + ref_bbox_width = bounds_4326[2] - bounds_4326[0] + ref_bbox_height = bounds_4326[3] - bounds_4326[1] + res_x = ref_bbox_width / ref_width + res_y = ref_bbox_height / ref_height + + minx = np.floor(minx / res_x) * res_x + miny = np.floor(miny / res_y) * res_y + maxx = np.ceil(maxx / res_x) * res_x + maxy = np.ceil(maxy / res_y) * res_y + + width = int(round((maxx - minx) / res_x)) + height = int(round((maxy - miny) / res_y)) + + dst_transform = from_bounds(minx, miny, maxx, maxy, width, height) + + return dst_transform, width, height, dst_crs diff --git a/src/satchip/modality.py b/src/satchip/modality.py new file mode 100644 index 0000000..4f13ebe --- /dev/null +++ b/src/satchip/modality.py @@ -0,0 +1,9 @@ +from dataclasses import dataclass + + +@dataclass(frozen=True) +class Modality: + id: str + all_bands: tuple[str, ...] + stack_bands: tuple[str, ...] + chip_bands: tuple[str, ...] diff --git a/src/satchip/models.py b/src/satchip/models.py new file mode 100644 index 0000000..a72984d --- /dev/null +++ b/src/satchip/models.py @@ -0,0 +1,166 @@ +from dataclasses import dataclass +from datetime import datetime +from pathlib import Path +from typing import NamedTuple, TypedDict + +import pyproj +import shapely +from rasterio.coords import BoundingBox + + +class Band(NamedTuple): + id: str + name: str + shortname: str + + +class Modality(TypedDict): + id: str + bands: tuple[Band, ...] + collection: str + + +@dataclass(frozen=True) +class Event: + name: str + date: datetime + wgs84_geometry: shapely.geometry.Polygon + buffer_m: int = 0 + + def buffered_geometry(self) -> shapely.geometry.Polygon: + if self.buffer_m == 0: + return self.wgs84_geometry + + to_utm = pyproj.Transformer.from_crs(4326, 32615, always_xy=True) + projected = shapely.ops.transform(to_utm.transform, self.wgs84_geometry) + buffered = projected.buffer(self.buffer_m) + + to_wgs84 = pyproj.Transformer.from_crs(32615, 4326, always_xy=True) + return shapely.ops.transform(to_wgs84.transform, buffered) + + +@dataclass(frozen=True) +class GridCell: + id: str + bounds: BoundingBox + + +@dataclass(frozen=True) +class Layer: + name: str + modality: Modality + path: Path + + +@dataclass(frozen=True) +class Chip: + id: str + path: Path + + +@dataclass(frozen=True) +class ChipStack: + id: str + data: Path + validation_mask: Path + label: Path + modality: Modality + + +class ModalityError(Exception): + pass + + +def bands_by_shortname(bands: tuple[Band, ...]) -> dict[str, Band]: + return {band.shortname: band for band in bands} + + +def bands_by_id(bands: tuple[Band, ...]) -> dict[str, Band]: + return {band.id: band for band in bands} + + +MODALITIES: dict[str, Modality] = { + 'OPERA_RTC': { + 'id': 'OPERA_RTC', + 'collection': 'OPERA_L2_RTC-S1_V1', + 'bands': ( + Band('VV', 'VV', 'VV'), + Band('VH', 'VH', 'VH'), + Band('mask', 'Validitiy Mask', 'mask'), + ), + }, + 'HLS_S30': { + 'id': 'HLS_S30', + 'collection': 'HLSS30', + 'bands': ( + Band('B02', 'Blue', 'B'), + Band('B03', 'Green', 'G'), + Band('B04', 'Red', 'R'), + Band('B8A', 'NIR Narrow', 'N'), + Band('B11', 'SWIR 1', 'SW1'), + Band('B12', 'SWIR 2', 'SW2'), + Band('Fmask', 'Cloud Mask', 'Fmask'), + ), + }, + 'HLS_L30': { + 'id': 'HLS_L30', + 'collection': 'HLSL30', + 'bands': ( + Band('B02', 'Blue', 'B'), + Band('B03', 'Green', 'G'), + Band('B04', 'Red', 'R'), + Band('B05', 'NIR Narrow', 'N'), + Band('B06', 'SWIR 1', 'SW1'), + Band('B07', 'SWIR 2', 'SW2'), + Band('Fmask', 'Cloud Mask', 'fmask'), + ), + }, +} + +MODALITY_IDS = list(MODALITIES.keys()) + +# https://hyp3-docs.asf.alaska.edu/guides/opera_rtc_product_guide/ +OPERA_RTC = MODALITIES['OPERA_RTC'] +OPERA_RTC_BANDS = bands_by_shortname(MODALITIES['OPERA_RTC']['bands']) + +# https://www.earthdata.nasa.gov/data/projects/hls/spectral-bands +HLS_S30 = MODALITIES['HLS_S30'] +HLS_S30_BANDS = bands_by_shortname(MODALITIES['HLS_S30']['bands']) + +HLS_L30 = MODALITIES['HLS_L30'] +HLS_L30_BANDS = bands_by_shortname(MODALITIES['HLS_L30']['bands']) + + +def band_id_from_filename(filename: str, modality_id: str) -> Band | None: + if 'HLS_L30' in modality_id: + sensor, band = band_from_hls_filename(filename) + elif 'HLS_S30' in modality_id: + sensor, band = band_from_hls_filename(filename) + elif 'OPERA_RTC' in modality_id: + sensor, band = band_from_rtc_filename(filename) + else: + raise ModalityError(f'Modality not found {modality_id}, must be ({MODALITY_IDS})') + + if sensor not in modality_id: + return '' + + bands = bands_by_id(MODALITIES[modality_id]['bands']) + + try: + return bands[band] + except KeyError: + return '' + + +def band_from_hls_filename(filename: str) -> tuple[str, str]: + # HLS.L30.T15TUG.2017167T165321.v2.0.B06.tif + parts = filename.split('.') + + sensor_key, band = parts[1], parts[-2] + + return sensor_key, band + + +def band_from_rtc_filename(filename: str) -> tuple[str, str]: + # OPERA_L2_RTC-S1_T063-133415-IW2_20170620T001327Z_20250925T045340Z_S1A_30_v1.0_VV.tif + return 'RTC', filename.split('_')[-1].split('.')[0] diff --git a/src/satchip/mosaic.py b/src/satchip/mosaic.py new file mode 100644 index 0000000..a5d1bc1 --- /dev/null +++ b/src/satchip/mosaic.py @@ -0,0 +1,265 @@ +import datetime +import shutil +from pathlib import Path + +import earthaccess +import hls +import numpy as np +import opera_rtc +import rasterio +from modality import Modality +from rasterio.crs import CRS +from rasterio.merge import merge +from rasterio.transform import from_bounds +from rasterio.warp import Resampling, calculate_default_transform, reproject, transform_bounds + + +def data_over_swath(swath, modalities: list[Modality], output_path: Path): + data_paths = { + 'RAW': output_path / 'RAW', + 'REPROJECTED': output_path / 'REPROJECTED', + 'MOSAIC': output_path / 'MOSAIC', + } + + for p in ('MOSAIC',): + shutil.rmtree(data_paths[p], ignore_errors=True) + + for p in data_paths.values(): + p.mkdir(parents=True, exist_ok=True) + + swathID = f'{int(swath["swathID"]):04d}' + + stacked = {} + for modality in modalities: + print(f'Localizing data for {modality.id}.') + + bounding_box = swath['buffered_event_background'].bounds + + start_date = swath['ls5hlsDate'] if modality.id == 'HLS' else swath['s1Date'] + + results = search_data(bounding_box, start_date, modality) + + local_files = earthaccess.download(results, local_path=data_paths['RAW'], show_progress=True) + + data_tifs = [f for f in local_files if f.name.endswith('.tif')] + reprojected_tifs = _reproject_files(data_tifs, output_path=data_paths['REPROJECTED']) + + if len(data_tifs) == 0: + print(f'Skipping: no data for swath {swathID}') + continue + + mod_merged = {} + + for band in modality.all_bands: + band_files = [f for f in reprojected_tifs if band in band_from_filename(f.name, modality)] + merged_name = make_merge_name(swathID, start_date, band, modality) + + merged_band_path = _merge(band_files, output_file=data_paths['MOSAIC'] / merged_name) + + mod_merged[band] = merged_band_path + + print(f'Stacking bands for {modality.id}') + stacked_data = _stack_bands(mod_merged, data_bands=modality.stack_bands, stacked_name='BANDS') + stacked[modality.id] = stacked_data + + print(f'Warp band {band} data to same area') + warped = _warp_over_swath( + data=stacked, + bounding_box_4326=bounding_box, + output_dir=output_path, + ) + + return warped + + +def search_data(bounding_box: tuple, start_date: datetime.datetime, modality: Modality): + results = {} + + if modality.id == 'HLS': + results = hls.search_hls_data(start_date=start_date, bounding_box=bounding_box) + elif modality.id == 'RTC': + results = opera_rtc.search_rtc_data(start_date=start_date, bounding_box=bounding_box) + + return results + + +def band_from_filename(filename, modality): + if modality.id == 'HLS': + band = hls.band_from_hls_filename(filename) + elif modality.id == 'RTC': + band = opera_rtc.band_from_rtc_filename(filename) + + return band + + +def make_merge_name(swathID: str, start_date: datetime.datetime, band: str, modality: Modality): + date_str = start_date.date().isoformat() + + return f'{swathID}.{modality.id}.{date_str}.{band}.tif' + + +def _reproject_files(files: list[Path], output_path: Path) -> list[Path]: + reprojected_paths = [output_path / f'{granule.name}' for granule in files] + + for granule, output_path in zip(files, reprojected_paths): + if output_path.exists(): + continue + + print(f'reprojecting to wgs84: {output_path.name}') + _reproject_file(granule, output_path) + + return reprojected_paths + + +def _reproject_file(local_file: Path, reprojected_file: Path, epsg=4326) -> None: + # https://rasterio.readthedocs.io/en/stable/topics/reproject.html#reprojecting-a-geotiff-dataset + with rasterio.open(local_file) as src: + dst_crs = CRS.from_epsg(epsg) + transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds) + + dst_kwargs = src.meta.copy() + dst_kwargs.update({'crs': dst_crs, 'transform': transform, 'width': width, 'height': height}) + + with rasterio.open(reprojected_file, 'w', **dst_kwargs) as dst: + for i in range(1, src.count + 1): + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs=dst_crs, + ) + + +def _merge(band_files: list[Path], output_file: Path) -> Path: + band_datasets = [rasterio.open(rtc_tif) for rtc_tif in band_files] + + master_crs = band_datasets[0].crs + for ds in band_datasets[1:]: + if ds.crs != master_crs: + ds.crs = master_crs + + try: + mosaic, out_trans = merge(band_datasets) + mosaic = np.squeeze(mosaic) + + out_meta = band_datasets[0].meta.copy() + + out_meta.update( + { + 'driver': 'GTiff', + 'height': mosaic.shape[0], + 'width': mosaic.shape[1], + 'transform': out_trans, + 'crs': band_datasets[0].crs, + } + ) + + with rasterio.open(output_file, 'w', **out_meta) as dst: + dst.write(mosaic, 1) + finally: + for ds in band_datasets: + ds.close() + + return output_file + + +def _stack_bands(merged: dict[str, Path], data_bands: tuple[str], stacked_name: str) -> None: + with rasterio.open(merged[data_bands[0]]) as src: + meta = src.meta.copy() + + band = data_bands[0] + meta.update(count=len(data_bands), dtype=np.float32) + stacked_file_name = _rename(merged[band], f'{band}.tif', f'{stacked_name}.tif') + + with rasterio.open(stacked_file_name, 'w', **meta) as dst: + for idx, band in enumerate(data_bands, start=1): + with rasterio.open(merged[band]) as src: + dst.write(src.read(1), idx) + + merged[stacked_name] = stacked_file_name + return merged + + +def _rename(path: Path, extension: str, mask_name: str) -> Path: + return path.parent / path.name.replace(extension, mask_name) + + +def _warp_over_swath(data, bounding_box_4326, output_dir): + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + reference_path = next(iter(next(iter(data.values())).values())) + dst_transform, width, height, dst_crs = _build_common_grid(bounding_box_4326, reference_path) + + output = {} + for sensor, bands in data.items(): + output[sensor] = {} + for band_name, input_path in bands.items(): + out_path = output_dir / input_path.name + _warp_single(input_path, out_path, dst_transform, width, height, dst_crs, band_name) + output[sensor][band_name] = out_path + + return output + + +def _warp_single(input_path, output_path, dst_transform, width, height, dst_crs, band_name): + CATEGORICAL_BANDS = {'Fmask', 'mask'} + resampling = Resampling.nearest if band_name in CATEGORICAL_BANDS else Resampling.bilinear + + with rasterio.open(input_path) as src: + dst_data = np.zeros((src.count, height, width), dtype=src.dtypes[0]) + + reproject( + source=rasterio.band(src, list(range(1, src.count + 1))), + destination=dst_data, + src_transform=src.transform, + src_crs=src.crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + resampling=resampling, + dst_nodata=src.nodata, + ) + + out_meta = src.meta.copy() + out_meta.update( + { + 'driver': 'GTiff', + 'height': height, + 'width': width, + 'transform': dst_transform, + 'crs': dst_crs, + } + ) + + with rasterio.open(output_path, 'w', **out_meta) as dest: + dest.write(dst_data) + + return output_path + + +def _build_common_grid(bounding_box_4326, reference_path): + dst_crs = CRS.from_epsg(4326) + minx, miny, maxx, maxy = bounding_box_4326 + + with rasterio.open(reference_path) as ref: + bounds_4326 = transform_bounds(ref.crs, dst_crs, *ref.bounds, densify_pts=21) + ref_width = ref.width + ref_height = ref.height + + ref_bbox_width = bounds_4326[2] - bounds_4326[0] + ref_bbox_height = bounds_4326[3] - bounds_4326[1] + res_x = ref_bbox_width / ref_width + res_y = ref_bbox_height / ref_height + + minx = np.floor(minx / res_x) * res_x + miny = np.floor(miny / res_y) * res_y + maxx = np.ceil(maxx / res_x) * res_x + maxy = np.ceil(maxy / res_y) * res_y + + width = int(round((maxx - minx) / res_x)) + height = int(round((maxy - miny) / res_y)) + dst_transform = from_bounds(minx, miny, maxx, maxy, width, height) + + return dst_transform, width, height, dst_crs diff --git a/src/satchip/old/chip_data.py b/src/satchip/old/chip_data.py new file mode 100644 index 0000000..ec13168 --- /dev/null +++ b/src/satchip/old/chip_data.py @@ -0,0 +1,180 @@ +import argparse +import glob +from collections import Counter +from datetime import datetime +from pathlib import Path + +import earthaccess +import numpy as np +import xarray as xr +from satchip.chip_hls import get_hls_data +from satchip.chip_hyp3s1rtc import get_rtc_paths_for_chips, get_s1rtc_chip_data +from satchip.chip_operas1rtc import get_operartc_data +from satchip.chip_sentinel2 import get_s2l2a_data +from satchip.terra_mind_grid import TerraMindChip, TerraMindGrid +from shapely.geometry import box +from tqdm import tqdm + +from satchip import utils + + +def fill_missing_times(data_chip: xr.DataArray, times: np.ndarray) -> xr.DataArray: + missing_times = np.setdiff1d(times, data_chip.time.data) + missing_shape = (len(missing_times), len(data_chip.band), data_chip.y.size, data_chip.x.size) + missing_data = xr.DataArray( + np.full(missing_shape, 0, dtype=data_chip.dtype), + dims=('time', 'band', 'y', 'x'), + coords={ + 'time': missing_times, + 'band': data_chip.band.data, + 'y': data_chip.y.data, + 'x': data_chip.x.data, + }, + ) + return xr.concat([data_chip, missing_data], dim='time').sortby('time') + + +def get_chips(label_paths: list[Path]) -> list[TerraMindChip]: + label_datasets = [utils.load_chip(label_path) for label_path in label_paths] + bounds = utils.get_overall_bounds([ds.bounds for ds in label_datasets]) + + buffered = box(*bounds).buffer(0.5).bounds + grid = TerraMindGrid(latitude_range=(buffered[1], buffered[3]), longitude_range=(buffered[0], buffered[2])) + grid_chips = {chip.name: chip for chip in grid.terra_mind_chips} + + chips = [] + for label_dataset in label_datasets: + label_chip_name = label_dataset.sample.item() + assert label_chip_name in grid_chips, f'No TerraMind chip found for label {label_chip_name}' + chip = grid_chips[label_chip_name] + chips.append(chip) + + return chips + + +def chip_data( + chip: TerraMindChip, + platform: str, + opts: utils.ChipDataOpts, + image_dir: Path, +) -> xr.Dataset: + if platform == 'HYP3S1RTC': + rtc_paths = opts['local_hyp3_paths'][chip.name] + chip_dataset = get_s1rtc_chip_data(chip, rtc_paths) + elif platform == 'S1RTC': + chip_dataset = get_operartc_data(chip, image_dir, opts=opts) + elif platform == 'S2L2A': + chip_dataset = get_s2l2a_data(chip, image_dir, opts=opts) + elif platform == 'HLS': + chip_dataset = get_hls_data(chip, image_dir, opts=opts) + else: + raise Exception(f'Unknown platform {platform}') + + return chip_dataset + + +def create_chips( + label_paths: list[Path], + platform: str, + dates: utils.DateRange | list[datetime], + strategy: str, + max_cloud_pct: int, + chip_dir: Path, + image_dir: Path, +) -> list[Path]: + platform_dir = chip_dir / platform + platform_dir.mkdir(parents=True, exist_ok=True) + + opts: utils.ChipDataOpts = {'strategy': strategy, 'dates': dates} + if platform in ['S2L2A', 'HLS']: + opts['max_cloud_pct'] = max_cloud_pct + + if platform in ['S1RTC', 'HLS']: + earthaccess.login() + + chips = get_chips(label_paths) + chip_names = [c.name for c in chips] + if len(chip_names) != len(set(chip_names)): + duplicates = [name for name, count in Counter(chip_names).items() if count > 1] + msg = f'Duplicate sample locations not supported. Duplicate chips: {", ".join(duplicates)}' + raise NotImplementedError(msg) + + chip_paths = [ + platform_dir / (x.with_suffix('').with_suffix('').name + f'_{platform}.zarr.zip') for x in label_paths + ] + + if platform == 'HYP3S1RTC': + rtc_paths_for_chips = get_rtc_paths_for_chips(chips, image_dir, opts) + opts['local_hyp3_paths'] = rtc_paths_for_chips + + for chip, chip_path in tqdm(list(zip(chips, chip_paths)), desc='Chipping labels'): + dataset = chip_data(chip, platform, opts, image_dir) + utils.save_chip(dataset, chip_path) + return chip_paths + + +def main() -> None: + parser = argparse.ArgumentParser(description='Chip a label image') + + parser.add_argument('labels', type=str, help='Path or Glob pattern for label chips') + + parser.add_argument( + 'platform', choices=['S1RTC', 'S2L2A', 'HLS', 'HYP3S1RTC'], type=str, help='Dataset to create chips for' + ) + + date_group = parser.add_mutually_exclusive_group(required=True) + + date_group.add_argument( + '--daterange', + nargs=2, + type=str, + metavar=('START', 'END'), + help='Inclusive date range in the format YYYY-mm-dd YYYY-mm-dd', + ) + date_group.add_argument( + '--dates', nargs='+', type=str, help='Space-separated list of specific dates in format YYYY-mm-dd' + ) + + parser.add_argument('--maxcloudpct', default=100, type=int, help='Maximum percent cloud cover for a data chip') + parser.add_argument('--chipdir', default='.', type=Path, help='Output directory for the chips') + parser.add_argument( + '--imagedir', default=None, type=Path, help='Output directory for image files. Defaults to chipdir/IMAGES' + ) + parser.add_argument( + '--strategy', + default='BEST', + choices=['BEST', 'ALL', 'SPECIFIC'], + type=str, + help='Strategy to use when multiple scenes are found (default: BEST)', + ) + args = parser.parse_args() + + args.platform = args.platform.upper() + assert 0 <= args.maxcloudpct <= 100, 'maxcloudpct must be between 0 and 100' + + dates: utils.DateRange | list[datetime] + if args.daterange: + start_date = datetime.strptime(args.daterange[0], '%Y-%m-%d') + end_date = datetime.strptime(args.daterange[1], '%Y-%m-%d') + assert start_date < end_date, 'start date must be before end date' + + dates = utils.DateRange(start_date, end_date) + else: + dates = [datetime.strptime(d, '%Y-%m-%d') for d in args.dates] + args.strategy = 'SPECIFIC' + + if '*' in args.labels or '?' in args.labels or '[' in args.labels: + label_paths = [Path(p) for p in glob.glob(args.labels)] + else: + label_paths = list(Path(args.labels).glob('*.zarr.zip')) + + assert len(label_paths) > 0, f'No label files found in {args.labels}' + + if args.imagedir is None: + args.imagedir = args.chipdir / 'IMAGES' + + create_chips(label_paths, args.platform, dates, args.strategy, args.maxcloudpct, args.chipdir, args.imagedir) + + +if __name__ == '__main__': + main() diff --git a/src/satchip/chip_hls.py b/src/satchip/old/chip_hls.py similarity index 84% rename from src/satchip/chip_hls.py rename to src/satchip/old/chip_hls.py index 4452356..5163f0d 100644 --- a/src/satchip/chip_hls.py +++ b/src/satchip/old/chip_hls.py @@ -8,13 +8,15 @@ import shapely import xarray as xr from earthaccess.results import DataGranule +from satchip.chip_xr_base import create_dataset_chip, create_template_da +from satchip.terra_mind_grid import TerraMindChip from shapely.geometry import Polygon from satchip import utils -from satchip.chip_xr_base import create_dataset_chip, create_template_da -from satchip.terra_mind_grid import TerraMindChip +earthaccess.login() + HLS_L_BANDS = OrderedDict( { 'B01': 'COASTAL', @@ -78,19 +80,25 @@ def get_scenes( Returns: The best HLS items. """ - assert strategy in ['BEST', 'ALL'], 'Strategy must be either BEST or ALL' + assert strategy in ['BEST', 'ALL', 'SPECIFIC'], 'Strategy must be either BEST or ALL' overlapping_items = [x for x in items if get_pct_intersect(x['umm'], roi) > 95] best_first = sorted(overlapping_items, key=lambda x: (-get_pct_intersect(x['umm'], roi), get_date(x['umm']))) + valid_scenes = [] + for item in best_first: product_id = get_product_id(item['umm']) n_products = len(list(image_dir.glob(f'{product_id}*'))) + if n_products < 15: earthaccess.download([item], image_dir, pqdm_kwargs={'disable': True}) + fmask_path = image_dir / f'{product_id}.v2.0.Fmask.tif' assert fmask_path.exists(), f'File not found: {fmask_path}' + qual_da = rioxarray.open_rasterio(fmask_path).rio.clip_box(*roi.bounds, crs='EPSG:4326') # type: ignore bit_masks = np.unpackbits(qual_da.data[0][..., np.newaxis], axis=-1) + # Looks for a 1 in the 4th, 6th and 7th bit of the Fmask (reverse order). See table 9 and appendix A of: # https://lpdaac.usgs.gov/documents/1698/HLS_User_Guide_V2.pdf bad_pixels = (bit_masks[..., 4] == 1) | (bit_masks[..., 6] == 1) | (bit_masks[..., 7] == 1) @@ -104,22 +112,38 @@ def get_scenes( return valid_scenes +def search_for_data(dates: utils.DateRange | list[datetime], bounds: utils.Bounds) -> list: + results = [] + if isinstance(dates, utils.DateRange): + results = earthaccess.search_data( + short_name=['HLSL30', 'HLSS30'], bounding_box=bounds, temporal=(dates.start, dates.end + timedelta(days=1)) + ) + else: + for date in dates: + day_results = earthaccess.search_data( + short_name=['HLSL30', 'HLSS30'], bounding_box=bounds, temporal=(date, date + timedelta(days=1)) + ) + results.extend(day_results) + + return results + + def get_hls_data(chip: TerraMindChip, image_dir: Path, opts: utils.ChipDataOpts) -> xr.Dataset: """Returns XArray DataArray of a Harmonized Landsat Sentinel-2 image for the given bounds and closest collection after date. """ - date_start = opts['date_start'] - date_end = opts['date_end'] + timedelta(days=1) # inclusive end - earthaccess.login() - results = earthaccess.search_data( - short_name=['HLSL30', 'HLSS30'], bounding_box=chip.bounds, temporal=(date_start, date_end) - ) - assert len(results) > 0, f'No HLS scenes found for chip {chip.name} between {date_start} and {date_end}.' + dates = opts['dates'] + + results = search_for_data(dates, utils.Bounds(*chip.bounds)) + assert len(results) > 0, f'No HLS scenes found for chip {chip.name} {utils.dates_error_msg(dates)}.' + roi = shapely.box(*chip.bounds) roi_buffered = roi.buffer(0.01) max_cloud_pct = opts.get('max_cloud_pct', 100) strategy = opts.get('strategy', 'BEST').upper() - timesteps = get_scenes(results, roi, max_cloud_pct, strategy, image_dir) + + timesteps = get_scenes(results, roi_buffered, max_cloud_pct, strategy, image_dir) + template = create_template_da(chip) timestep_arrays = [] for scene in timesteps: diff --git a/src/satchip/chip_hyp3s1rtc.py b/src/satchip/old/chip_hyp3s1rtc.py similarity index 100% rename from src/satchip/chip_hyp3s1rtc.py rename to src/satchip/old/chip_hyp3s1rtc.py index f94625c..c349737 100644 --- a/src/satchip/chip_hyp3s1rtc.py +++ b/src/satchip/old/chip_hyp3s1rtc.py @@ -7,11 +7,11 @@ import rioxarray import shapely import xarray as xr - -from satchip import utils from satchip.chip_xr_base import create_dataset_chip, create_template_da from satchip.terra_mind_grid import TerraMindChip +from satchip import utils + S1RTC_BANDS = ('VV', 'VH') diff --git a/src/satchip/chip_label.py b/src/satchip/old/chip_label.py similarity index 100% rename from src/satchip/chip_label.py rename to src/satchip/old/chip_label.py index 246f8e0..6234a9f 100644 --- a/src/satchip/chip_label.py +++ b/src/satchip/old/chip_label.py @@ -5,11 +5,11 @@ import numpy as np import rasterio as rio import xarray as xr +from satchip.chip_xr_base import create_dataset_chip +from satchip.terra_mind_grid import TerraMindGrid from tqdm import tqdm from satchip import utils -from satchip.chip_xr_base import create_dataset_chip -from satchip.terra_mind_grid import TerraMindGrid def is_valuable(chip: np.ndarray) -> bool: diff --git a/src/satchip/chip_operas1rtc.py b/src/satchip/old/chip_operas1rtc.py similarity index 80% rename from src/satchip/chip_operas1rtc.py rename to src/satchip/old/chip_operas1rtc.py index 5811998..4eb8151 100644 --- a/src/satchip/chip_operas1rtc.py +++ b/src/satchip/old/chip_operas1rtc.py @@ -8,12 +8,12 @@ import xarray as xr from earthaccess.results import DataGranule from osgeo import gdal - -from satchip import utils from satchip.chip_hls import get_geometry, get_product_id from satchip.chip_xr_base import create_dataset_chip, create_template_da from satchip.terra_mind_grid import TerraMindChip +from satchip import utils + gdal.UseExceptions() @@ -87,25 +87,54 @@ def get_scenes(groups: list[RTCGroup], roi: shapely.geometry.Polygon, strategy: return intersecting[:1] elif strategy == 'ALL': return intersecting + elif strategy == 'SPECIFIC': + return intersecting else: raise ValueError(f'Strategy must be either BEST or ALL. Got {strategy}') +@utils.retry_on_connection_error(max_retries=5, backoff_factor=2) +def search_for_data(dates: utils.DateRange | list[datetime], bounds: utils.Bounds) -> list: + results = [] + + if isinstance(dates, utils.DateRange): + results = earthaccess.search_data( + short_name=['OPERA_L2_RTC-S1_V1'], + bounding_box=bounds, + temporal=(dates.start, dates.end + timedelta(days=1)), + ) + else: + for date in dates: + day_results = earthaccess.search_data( + short_name=['OPERA_L2_RTC-S1_V1'], bounding_box=bounds, temporal=(date, date + timedelta(days=1)) + ) + results.extend(day_results) + + return results + + def get_operartc_data(chip: TerraMindChip, image_dir: Path, opts: utils.ChipDataOpts) -> xr.Dataset: """Returns XArray DataArray of a OPERA S1-RTC for the given chip and selection startegy.""" - date_start = opts['date_start'] - date_end = opts['date_end'] + timedelta(days=1) # inclusive end - earthaccess.login() - results = earthaccess.search_data( - short_name=['OPERA_L2_RTC-S1_V1'], bounding_box=chip.bounds, temporal=(date_start, date_end) - ) - results = filter_to_dualpol(results) - rtc_groups = group_rtcs(results) + + dates = opts['dates'] + roi = shapely.box(*chip.bounds) roi_buffered = roi.buffer(0.01) + + results = search_for_data(dates, utils.Bounds(*roi_buffered.bounds)) + dualpol = filter_to_dualpol(results) + + rtc_groups = group_rtcs(dualpol) strategy = opts.get('strategy', 'BEST').upper() timesteps = get_scenes(rtc_groups, roi_buffered, strategy) - assert len(timesteps) > 0, f'No OPERA RTC scenes found for chip {chip.name} between {date_start} and {date_end}.' + + assert len(timesteps) > 0, f'No OPERA RTC scenes found for chip {chip.name} {utils.dates_error_msg(dates)}' + + if isinstance(dates, list): + print(dates, timesteps) + breakpoint() + assert len(timesteps) == len(dates) + vrts = [timestep.download(image_dir) for timestep in timesteps] template = create_template_da(chip) timestep_arrays = [] diff --git a/src/satchip/chip_sentinel2.py b/src/satchip/old/chip_sentinel2.py similarity index 85% rename from src/satchip/chip_sentinel2.py rename to src/satchip/old/chip_sentinel2.py index f4725b4..39370ed 100644 --- a/src/satchip/chip_sentinel2.py +++ b/src/satchip/old/chip_sentinel2.py @@ -10,11 +10,11 @@ import xarray as xr from pystac.item import Item from pystac_client import Client - -from satchip import utils from satchip.chip_xr_base import create_dataset_chip, create_template_da from satchip.terra_mind_grid import TerraMindChip +from satchip import utils + S2_BANDS = OrderedDict( { @@ -91,11 +91,11 @@ def get_scenes( The best Sentinel-2 L2A item. """ strategy = strategy.upper() - assert strategy in ['BEST', 'ALL'], 'Strategy must be either BEST or ALL' assert len(items) > 0, 'No Sentinel-2 L2A scenes found for chip.' items = [item for item in items if get_pct_intersect(item.geometry, roi) > 0.95] best_first = sorted(items, key=lambda x: (-get_pct_intersect(x.geometry, roi), x.datetime)) valid_scenes = [] + for item in best_first: scl_href = item.assets['scl'].href local_path = fetch_s3_file(scl_href, image_dir) @@ -131,6 +131,36 @@ def get_s2_version(item: Item) -> int: return latest_items +def search_for_data(dates: utils.DateRange | list[datetime], roi: shapely.Polygon) -> list: + results = [] + client = Client.open('https://earth-search.aws.element84.com/v1') + + if isinstance(dates, utils.DateRange): + date_end = dates.end + timedelta(days=1) + date_range = f'{datetime.strftime(dates.start, "%Y-%m-%d")}/{datetime.strftime(date_end, "%Y-%m-%d")}' + search = client.search( + collections=['sentinel-2-l2a'], + intersects=roi, + datetime=date_range, + max_items=1000, + ) + + results = list(search.item_collection()) + else: + for date in dates: + search = client.search( + collections=['sentinel-2-l2a'], + intersects=roi, + datetime=datetime.strftime(date, '%Y-%m-%d'), + max_items=1000, + ) + results.extend(search.item_collection()) + + results = get_latest_image_versions(results) + + return results + + def get_s2l2a_data(chip: TerraMindChip, image_dir: Path, opts: utils.ChipDataOpts) -> xr.Dataset: """Get XArray DataArray of Sentinel-2 L2A image for the given bounds and best collection parameters. @@ -146,28 +176,20 @@ def get_s2l2a_data(chip: TerraMindChip, image_dir: Path, opts: utils.ChipDataOpt Returns: XArray Dataset containing the Sentinel-2 L2A image data. """ - date_start = opts['date_start'] - date_end = opts['date_end'] + timedelta(days=1) # inclusive end - date_range = f'{datetime.strftime(date_start, "%Y-%m-%d")}/{datetime.strftime(date_end, "%Y-%m-%d")}' + dates = opts['dates'] + roi = shapely.box(*chip.bounds) roi_buffered = roi.buffer(0.01) - client = Client.open('https://earth-search.aws.element84.com/v1') - search = client.search( - collections=['sentinel-2-l2a'], - intersects=roi, - datetime=date_range, - max_items=1000, - ) - assert len(search.item_collection()) > 0, ( - f'No Sentinel-2 L2A scenes found for chip {chip.name} between {date_start} and {date_end}.' - ) - assert len(search.item_collection()) < 1000, ( - 'Too many Sentinel-2 L2A scenes found for chip. Please narrow the date range.' - ) - items = list(search.item_collection()) - items = get_latest_image_versions(items) + + items = search_for_data(dates, roi) max_cloud_pct = opts.get('max_cloud_pct', 100) strategy = opts.get('strategy', 'BEST') + + assert len(items) > 0, ( + f'No Sentinel-2 L2A scenes found for chip {chip.name} between {utils.dates_error_msg(dates)}.' + ) + assert len(items) < 1000, 'Too many Sentinel-2 L2A scenes found for chip. Please narrow the date range.' + timesteps = get_scenes(items, roi, strategy, max_cloud_pct, image_dir) urls = [item.assets[band.lower()].href for item in timesteps for band in S2_BANDS.values()] diff --git a/src/satchip/chip_view.py b/src/satchip/old/chip_view.py similarity index 100% rename from src/satchip/chip_view.py rename to src/satchip/old/chip_view.py diff --git a/src/satchip/chip_xr_base.py b/src/satchip/old/chip_xr_base.py similarity index 100% rename from src/satchip/chip_xr_base.py rename to src/satchip/old/chip_xr_base.py index b441ef1..e9541e2 100644 --- a/src/satchip/chip_xr_base.py +++ b/src/satchip/old/chip_xr_base.py @@ -3,9 +3,9 @@ import numpy as np import xarray as xr +from satchip.terra_mind_grid import TerraMindChip import satchip -from satchip.terra_mind_grid import TerraMindChip def _check_spec(dataset: xr.Dataset) -> None: diff --git a/src/satchip/major_tom_grid.py b/src/satchip/old/major_tom_grid.py similarity index 100% rename from src/satchip/major_tom_grid.py rename to src/satchip/old/major_tom_grid.py diff --git a/src/satchip/terra_mind_grid.py b/src/satchip/old/terra_mind_grid.py similarity index 100% rename from src/satchip/terra_mind_grid.py rename to src/satchip/old/terra_mind_grid.py index 71596b6..cc4876d 100644 --- a/src/satchip/terra_mind_grid.py +++ b/src/satchip/old/terra_mind_grid.py @@ -3,8 +3,8 @@ import numpy as np import pyproj from rasterio import Affine - from satchip.major_tom_grid import MajorTomGrid + from satchip.utils import get_epsg4326_bbox, get_epsg4326_point diff --git a/src/satchip/old/utils.py b/src/satchip/old/utils.py new file mode 100644 index 0000000..fd2467b --- /dev/null +++ b/src/satchip/old/utils.py @@ -0,0 +1,114 @@ +import datetime +import functools +import time +import warnings +from collections.abc import Callable +from pathlib import Path +from typing import NamedTuple, ParamSpec, TypeVar, TypedDict + +import xarray as xr +import zarr +from pyproj import CRS, Transformer +from requests.exceptions import ConnectionError + + +class RtcImageSet(TypedDict): + VV: Path + VH: Path + + +class Bounds(NamedTuple): + minx: float + miny: float + maxx: float + maxy: float + + +class DateRange(NamedTuple): + start: datetime.datetime + end: datetime.datetime + + +class ChipDataRequiredOpts(TypedDict): + strategy: str + dates: DateRange | list[datetime.datetime] + + +class ChipDataOpts(ChipDataRequiredOpts, total=False): + max_cloud_pct: int + local_hyp3_paths: dict[str, list[RtcImageSet]] + + +def dates_error_msg(dates: DateRange | list[datetime.datetime]) -> str: + return f'between {dates.start} and {dates.end}' if isinstance(dates, DateRange) else f'for dates {dates}' + + +def get_overall_bounds(bounds: list) -> Bounds: + minx = min([b[0] for b in bounds]) + miny = min([b[1] for b in bounds]) + maxx = max([b[2] for b in bounds]) + maxy = max([b[3] for b in bounds]) + return Bounds(minx, miny, maxx, maxy) + + +def get_epsg4326_point(x: float, y: float, in_epsg: int) -> tuple[float, float]: + if in_epsg == 4326: + return x, y + in_crs = CRS.from_epsg(in_epsg) + out_crs = CRS.from_epsg(4326) + transformer = Transformer.from_crs(in_crs, out_crs, always_xy=True) + newx, newy = transformer.transform(x, y) + return round(newx, 5), round(newy, 5) + + +def get_epsg4326_bbox( + bounds: tuple[float, float, float, float], in_epsg: int, buffer: float = 0.1 +) -> tuple[float, float, float, float]: + minx, miny = get_epsg4326_point(bounds[0], bounds[1], in_epsg) + maxx, maxy = get_epsg4326_point(bounds[2], bounds[3], in_epsg) + bbox = minx - buffer, miny - buffer, maxx + buffer, maxy + buffer + return bbox + + +def save_chip(dataset: xr.Dataset, save_path: str | Path) -> None: + """Save a zipped zarr archive""" + store = zarr.storage.ZipStore(save_path, mode='w') + with warnings.catch_warnings(): + warnings.filterwarnings('ignore', message='Duplicate name:', module='zipfile') + dataset.to_zarr(store) # type: ignore[call-overload] + store.close() + + +def load_chip(label_path: str | Path) -> xr.Dataset: + """Load a zipped zarr archive""" + store = zarr.storage.ZipStore(label_path, read_only=True) + dataset = xr.open_zarr(store) + return dataset + + +P = ParamSpec('P') +R = TypeVar('R') + + +def retry_on_connection_error( + max_retries: int = 3, backoff_factor: float = 1 +) -> Callable[[Callable[P, R]], Callable[P, R]]: + def decorator(func: Callable[P, R]) -> Callable[P, R]: + @functools.wraps(func) + def wrapper(*args: P.args, **kwargs: P.kwargs) -> R: + for attempt in range(max_retries): + try: + return func(*args, **kwargs) + except (ConnectionError, OSError) as e: + if attempt == max_retries - 1: + print(f'Failed after {max_retries} attempts: {e}') + raise e + + wait_time = backoff_factor * (2**attempt) + time.sleep(wait_time) + + raise RuntimeError('Unexpected exit from retry loop') + + return wrapper + + return decorator diff --git a/src/satchip/opera_rtc.py b/src/satchip/opera_rtc.py new file mode 100644 index 0000000..f67768d --- /dev/null +++ b/src/satchip/opera_rtc.py @@ -0,0 +1,104 @@ +from datetime import datetime, timedelta +from pathlib import Path + +import earthaccess +import numpy as np +import rasterio +from earthaccess.results import DataGranule + + +def search_rtc_data(start_date: datetime, bounding_box: tuple[float, float, float, float]) -> list[DataGranule]: + final_date = start_date + timedelta(days=1) + + results = earthaccess.search_data( + short_name=['OPERA_L2_RTC-S1_V1'], + temporal=(start_date.strftime('%Y-%m-%d'), final_date.strftime('%Y-%m-%d')), + bounding_box=bounding_box, + ) + + return results + + +def band_from_rtc_filename(filename): + # OPERA_L2_RTC-S1_T063-133415-IW2_20170620T001327Z_20250925T045340Z_S1A_30_v1.0_VV.tif + return filename.split('_')[-1].split('.')[0] + + +def make_merged_rtc_name(template_filename: str) -> str: + """ + https://hyp3-docs.asf.alaska.edu/guides/opera_rtc_product_guide/#naming-convention + swathID.OPERA_L2_RTC-S1_[BurstID]_[StartDateTime]_[ProductGenerationDateTime] _[Sensor]_[PixelSpacing]_[ProductVersion]_[LayerName].Ext + + Input: 1442.OPERA_L2_RTC-S1_T063-133415-IW2_20170620T001327Z_20250925T045340Z_S1A_30_v1.0_VV.tif + Returns: 1442.OPERA_L2_RTC-133415-IW2_20170620_S1A_30_v1.0_VV.tif + """ + + # ['1442.OPERA', 'L2', 'RTC-S1', 'T063-133415-IW2', '20170620T001327Z', '20250925T045340Z', 'S1A', '30', 'v1.0', 'VV.tif'] + name_parts = template_filename.split('_') + + name_parts.pop(5) # Remove Product Generation Time + name_parts.pop(3) # Remove Burst ID + + return '_'.join(name_parts) + + +def is_valid_rtc(mask_path: Path, label_path: Path) -> bool: + with rasterio.open(mask_path) as ds: + validity_mask = ds.read(1) + + with rasterio.open(label_path) as ds: + event_mask = ds.read(1) + + is_event_pixel = event_mask == 1 + # https://hyp3-docs.asf.alaska.edu/guides/opera_rtc_product_guide/#validity-mask + is_valid_pixel = np.isin(validity_mask, [0, 1]) + + total_event_pixels = is_event_pixel.sum() + valid_event_pixels = (is_event_pixel & is_valid_pixel).sum() + + pct_valid_data = 100.0 * valid_event_pixels / total_event_pixels + print(f'Percent of the event with valid data: {pct_valid_data:.1f}%') + + return pct_valid_data > 50.0 + + +def filter_rtc_chips(chips: dict[str, dict]) -> list[dict]: + good_chips = [] + + for tile_id, chip in chips.items(): + with rasterio.open(chip['BANDS']) as ds: + rtc_data = ds.read() + + with rasterio.open(chip['EVENT']) as ds: + event_mask = ds.read(1) + + has_nan_pixels = np.isnan(rtc_data).sum() > 0 + + num_pixels = event_mask.size + num_event_pixels = np.count_nonzero(event_mask > 0) + + pct_pixels_over_event = 100.0 * (num_event_pixels / num_pixels) + data_overlaps_event = pct_pixels_over_event > 1 + + if not has_nan_pixels and data_overlaps_event: + good_chips.append(chip) + + return good_chips + + +def normalize_image_array(input_array: np.ndarray, vmin: float, vmax: float) -> np.ndarray: + input_array = input_array.astype(float) + scaled_array = (input_array - vmin) / (vmax - vmin) + scaled_array[np.isnan(input_array)] = 0 + normalized_array = np.round(np.clip(scaled_array, 0, 1) * 255).astype(np.uint8) + + return normalized_array + + +def get_rtc_img(rtc_data: np.ndarray) -> np.ndarray: + vv = normalize_image_array(np.sqrt(rtc_data[0]), 0.14, 0.52) + vh = normalize_image_array(np.sqrt(rtc_data[1]), 0.05, 0.259) + + img = np.stack([vv, vh, vv], axis=-1) + + return img diff --git a/src/satchip/view.py b/src/satchip/view.py new file mode 100644 index 0000000..0de8d2f --- /dev/null +++ b/src/satchip/view.py @@ -0,0 +1,234 @@ +from pathlib import Path + +import rasterio +from rasterio.merge import merge +import numpy as np +import matplotlib.pyplot as plt +import cartopy.crs as ccrs +from shapely.geometry import box + +from satchip import models, merge_modality + + +def view_merged( + stacked_data_file: Path, + event: models.Event, + modality: models.Modality, + rgb_bands: tuple[int, int, int], + save_to_file: Path | None = None, + quite: bool = False, +): + save_to_file.parent.mkdir(exist_ok=True, parents=True) + crs_pc = ccrs.PlateCarree() + + with rasterio.open(stacked_data_file) as ds: + bounds = ds.bounds + full_extent = [bounds.left, bounds.right, bounds.bottom, bounds.top] + band_data = ds.read() + + img = get_img(band_data, modality, rgb_bands) + + # plot BANDS and geom + fig, ax = plt.subplots( + 1, + 1, + subplot_kw={'projection': crs_pc}, + figsize=(12, 12), + layout='constrained', + ) + + event_geom = event.wgs84_geometry + + ax.imshow(img, extent=full_extent, origin='upper', transform=crs_pc) + ax.add_geometries([event_geom], edgecolor='red', linewidth=2, facecolor='none', crs=crs_pc) + + if save_to_file: + plt.savefig( + save_to_file, + dpi=300, + bbox_inches='tight', + ) + + if not quite: + plt.show() + + plt.close(fig) + + +def view_chip( + chip: models.ChipStack, + modality: models.Modality, + rgb_bands: tuple[int, int, int], + save_to_file: Path | None = None, + quite: bool = False, +): + save_to_file.parent.mkdir(exist_ok=True, parents=True) + + with rasterio.open(chip.data) as ds: + bounds = ds.bounds + full_extent = [bounds.left, bounds.right, bounds.bottom, bounds.top] + band_data = ds.read() + + with rasterio.open(chip.label) as ds: + label_data = ds.read().squeeze() + + img = get_img(band_data, modality, rgb_bands) + + crs_pc = ccrs.PlateCarree() + fig, ax = plt.subplots( + 1, + 2, + subplot_kw={'projection': crs_pc}, + figsize=(12, 12), + layout='constrained', + ) + + ax[0].imshow(img, extent=full_extent, origin='upper', transform=crs_pc) + ax[1].imshow(label_data, extent=full_extent, origin='upper', transform=crs_pc) + + if save_to_file: + plt.savefig( + save_to_file, + dpi=300, + bbox_inches='tight', + ) + + if not quite: + plt.show() + + plt.close(fig) + + +def view_chips( + chips: models.ChipStack, + modality: models.Modality, + rgb_bands: tuple[int, int, int], + save_to_file: Path | None = None, + quite: bool = False, +): + save_to_file.parent.mkdir(exist_ok=True, parents=True) + + merged_data = _merge_chips([chip.data for chip in chips], Path.cwd() / 'merged_data.tif') + merged_label = _merge_chips([chip.label for chip in chips], Path.cwd() / 'merged_label.tif') + try: + with rasterio.open(merged_data) as ds: + bounds = ds.bounds + full_extent = [bounds.left, bounds.right, bounds.bottom, bounds.top] + band_data = ds.read() + + with rasterio.open(merged_label) as ds: + label_data = ds.read().squeeze() + + img = get_img(band_data, modality, rgb_bands) + + crs_pc = ccrs.PlateCarree() + # plot BANDS and geom + fig, ax = plt.subplots( + 1, + 2, + subplot_kw={'projection': crs_pc}, + figsize=(12, 12), + layout='constrained', + ) + + ax[0].imshow(img, extent=full_extent, origin='upper', transform=crs_pc) + ax[1].imshow(label_data, extent=full_extent, origin='upper', transform=crs_pc) + + if save_to_file: + plt.savefig( + save_to_file, + dpi=300, + bbox_inches='tight', + ) + + if not quite: + plt.show() + + plt.close(fig) + finally: + merged_data.unlink(missing_ok=True) + merged_label.unlink(missing_ok=True) + + +def _merge_chips( + chips: list[Path], + output_file: Path, +) -> Path: + datasets = [rasterio.open(p) for p in chips] + + try: + mosaic, transform = merge(datasets) + + out_meta = datasets[0].meta.copy() + out_meta.update( + { + 'height': mosaic.shape[1], + 'width': mosaic.shape[2], + 'transform': transform, + } + ) + + with rasterio.open(output_file, 'w', **out_meta) as dst: + if len(mosaic.shape) == 2: + dst.write(mosaic, 1) + else: + dst.write(mosaic) + finally: + for ds in datasets: + ds.close() + + return output_file + + +def get_img(band_data: np.ndarray, modality: models.Modality, rgb_bands: tuple[int, int, int]): + if 'HLS' in modality['id']: + img = get_hls_img(band_data, rgb_bands) + elif 'RTC' in modality['id']: + img = get_rtc_img(band_data, rgb_bands) + + return img + + +def get_hls_img(hls_data: np.ndarray, rgb_bands: tuple[int, int, int]) -> np.ndarray: + r_band, g_band, b_band = rgb_bands + + r = bytescale(np.sqrt(np.clip(hls_data[r_band] / 10000.0, 0, 2)), 0, 0.5) + g = bytescale(np.sqrt(np.clip(hls_data[g_band] / 10000.0, 0, 2)), 0, 0.5) + b = bytescale(np.sqrt(np.clip(hls_data[b_band] / 10000.0, 0, 2)), 0, 0.5) + + rgb = np.dstack((r, g, b)) + return rgb + + +def bytescale(arr, cmin=0, cmax=1, low=0, high=255): + # clip the data to be in the range of cmin to cmax + arr = np.clip(arr, cmin, cmax) + high = float(high) + low = float(low) + cmax = float(cmax) + cmin = float(cmin) + m = (high - low) / (cmax - cmin) # slope + b = high - (m * cmax) # intercept + arr = np.uint8((m * arr) + b) + return arr + + +def get_rtc_img(rtc_data: np.ndarray, rgb_bands: tuple[int, int, int]) -> np.ndarray: + r_band, g_band, b_band = rgb_bands + + r = normalize_image_array(np.sqrt(rtc_data[r_band]), 0.14, 0.52) + g = normalize_image_array(np.sqrt(rtc_data[g_band]), 0.05, 0.259) + b = normalize_image_array(np.sqrt(rtc_data[b_band]), 0.14, 0.52) + + img = np.stack([r, g, b], axis=-1) + + return img + + +def normalize_image_array(input_array: np.ndarray, vmin: float, vmax: float) -> np.ndarray: + input_array = input_array.astype(float) + scaled_array = (input_array - vmin) / (vmax - vmin) + scaled_array[np.isnan(input_array)] = 0 + normalized_array = np.round(np.clip(scaled_array, 0, 1) * 255).astype(np.uint8) + + return normalized_array diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..868682c --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,77 @@ +from datetime import datetime +from pathlib import Path + +import geopandas as gpd +import pandas as pd +import pytest +from shapely import wkt + +from satchip import download_data, models + + +DATA_PATH = Path(__file__).parent / 'data' + + +def pytest_addoption(parser): + parser.addoption('--download', action='store_true', default=False) + + +def pytest_configure(config): + config.addinivalue_line('markers', 'download: marks tests that download data') + + +def pytest_collection_modifyitems(config, items): + if not config.getoption('--download'): + skip = pytest.mark.skip(reason='pass --download to run') + for item in items: + if item.get_closest_marker('download'): + item.add_marker(skip) + + +@pytest.fixture +def pristine_gdf(tmp_path): + shp_path = DATA_PATH / 'hwds_pristine' + df = gpd.read_file(shp_path) + + df['SwathDate'] = pd.to_datetime(df['SwathDate'], format='%Y-%m-%d') + df['HLSDate'] = pd.to_datetime(df['HLSDate'], format='%Y-%m-%d') + + return df + + +@pytest.fixture(scope='session') +def warped_event_files(): + warped_base = DATA_PATH / 'warped' + + return [ + warped_base / '102a.HLS_S30.2019-07-09.Fmask.tif', + warped_base / '102a.MASK.tif', + warped_base / '102a.stacked.tif', + ] + + +@pytest.fixture(scope='session') +def s2_local_files(s2_event): + download_path = DATA_PATH / 'raw' + print('Downloading test data...') + + local_files = download_data.download_data(s2_event, models.HLS_S30, download_path) + + return local_files + + +@pytest.fixture(scope='session') +def s2_event() -> models.Event: + geom = wkt.loads( + 'POLYGON Z ((-97.46296108799999 41.72598898100006 0, -97.45654663199997 41.72832979400005 0, -97.45286506599996 41.735911517000034 0, -97.44692854399995 41.73806293200005 0, -97.44255668499994 41.741031193000026 0, -97.43649360699999 41.746392472000025 0, -97.43643518499994 41.74966361400004 0, -97.44477798299994 41.751269181000055 0, -97.44546922799998 41.757304729000055 0, -97.44413504799996 41.76092731800003 0, -97.43601483299994 41.761099638000076 0, -97.43493753799999 41.76400779800008 0, -97.43574582899998 41.76925407400006 0, -97.43507482199999 41.775018801000044 0, -97.42929198599995 41.77529913700005 0, -97.42434097299997 41.77189406900004 0, -97.41643891099994 41.77095250000008 0, -97.40855823399994 41.772049811000045 0, -97.40518880099995 41.769938841000055 0, -97.40212726299995 41.772478093000075 0, -97.39708935799996 41.77557899800007 0, -97.39216022699998 41.77683695500008 0, -97.38494019299998 41.779687625000065 0, -97.37802386699997 41.78005042400008 0, -97.37595922599996 41.772699628000055 0, -97.36870573799996 41.772803176000025 0, -97.36473372899997 41.776025134000065 0, -97.35972290999996 41.77556353600005 0, -97.35569283499996 41.77524485300006 0, -97.35141833499995 41.77176253300007 0, -97.35179576499996 41.769616866000035 0, -97.35247585899998 41.76537635100004 0, -97.35455601099994 41.765104164000036 0, -97.35748216699994 41.763087014000064 0, -97.35755091099998 41.760709776000056 0, -97.35659214299994 41.75815047000003 0, -97.34941309099997 41.75410074900003 0, -97.34425890099999 41.75299627900006 0, -97.33955929299998 41.75280198300004 0, -97.33427827799994 41.75626863400004 0, -97.32638406299998 41.75688761500004 0, -97.32401029299996 41.75030413500008 0, -97.32617321299995 41.74390741500008 0, -97.33893306699997 41.74281858000006 0, -97.33915551999996 41.73881615700003 0, -97.33542314499994 41.73642923700004 0, -97.33487090999995 41.73053873300006 0, -97.33494400699999 41.72165371600005 0, -97.34530133099997 41.72321647500007 0, -97.34362613199994 41.72861741600008 0, -97.34850260899998 41.72835093900005 0, -97.34904493499994 41.72262336800003 0, -97.34959717099997 41.71581247300003 0, -97.35659214299994 41.713051300000075 0, -97.35972147399997 41.70716079600004 0, -97.36818907299994 41.707897109000044 0, -97.37554466699999 41.70552724600003 0, -97.37794521999996 41.71268314300005 0, -97.38066826999994 41.713925654000036 0, -97.38197931699995 41.71100923000006 0, -97.38053398199997 41.702322166000044 0, -97.37438188299996 41.70104730000003 0, -97.37256136399998 41.69851681000006 0, -97.37141382499999 41.691215095000075 0, -97.37359654199997 41.685445995000066 0, -97.37716937399995 41.686764063000055 0, -97.37836702299995 41.692105638000044 0, -97.38678588699997 41.69234998300004 0, -97.38711030299999 41.69474924900004 0, -97.38893657299997 41.696020211000075 0, -97.39203310699997 41.69605112600004 0, -97.39708425799995 41.69660031600006 0, -97.39699799999994 41.69790345900003 0, -97.39728821299997 41.70098871600004 0, -97.39958472399996 41.70204029100006 0, -97.40051777999997 41.703388508000046 0, -97.40357986399994 41.71025483200003 0, -97.40555695699999 41.71360353400007 0, -97.41463431499994 41.71473101400005 0, -97.41599053299996 41.71041910900004 0, -97.40847919899994 41.70361729000007 0, -97.40905360899995 41.698218834000045 0, -97.42033287699996 41.69235264400004 0, -97.42815022299999 41.69435335000003 0, -97.42415654299998 41.69969951000007 0, -97.42249215599998 41.70624040400003 0, -97.42727819099997 41.70936973400006 0, -97.43331385799996 41.70941182400003 0, -97.43942735499996 41.70716079600004 0, -97.44623824999996 41.70182252600006 0, -97.45065612799999 41.700902135000035 0, -97.45614052299999 41.70244454500005 0, -97.45855998799999 41.70632093900008 0, -97.46176883199996 41.71106968300006 0, -97.46907545699997 41.71252207500004 0, -97.47778018699995 41.71142744600007 0, -97.48360613399996 41.71139459500006 0, -97.48556606699998 41.71563875600003 0, -97.47960867399996 41.71848116100006 0, -97.47027376699998 41.720714549000036 0, -97.46593462299995 41.726661585000045 0, -97.46296108799999 41.72598898100006 0))' + ) + name = '95a' + date = datetime(year=2018, month=7, day=2) + + event = models.Event( + name=name, + date=date, + wgs84_geometry=geom, + ) + + return event diff --git a/tests/test_chip_data.py b/tests/test_chip_data.py new file mode 100644 index 0000000..a5feed4 --- /dev/null +++ b/tests/test_chip_data.py @@ -0,0 +1,72 @@ +import rasterio + +from satchip import chip_data, models + + +def test_make_grid(warped_event_files): + label = [f for f in warped_event_files if 'MASK' in f.name].pop() + + grid = chip_data.make_grid_from_reference(label) + assert len(grid) == 27 + + grid = chip_data.make_grid_from_reference(label, chip_size=512) + assert len(grid) == 4 + + grids = (tuple(chip_data.make_grid_from_reference(f)) for f in warped_event_files) + assert len(set(grids)) == 1 + + +def test_chip_data(warped_event_files, tmp_path): + label = [f for f in warped_event_files if 'MASK' in f.name].pop() + grid = chip_data.make_grid_from_reference(label) + + for file in warped_event_files: + chips = chip_data.chip_data(grid, file, tmp_path / 'chips') + + chip = chips[0] + assert file.name in chip.path.name + assert chip.id in chip.path.name + + assert len(chips) == len(grid) + + shapes = set() + for chip in chips: + with rasterio.open(chip.path) as ds: + shapes.add(ds.shape) + + assert shapes == {(256, 256)} + + +def test_make_chip_stacks(warped_event_files, tmp_path): + fmask, label, data = warped_event_files + grid = chip_data.make_grid_from_reference(fmask) + + fmask_chips = chip_data.chip_data(grid, fmask, tmp_path / 'chips') + label_chips = chip_data.chip_data(grid, label, tmp_path / 'chips') + data_chips = chip_data.chip_data(grid, data, tmp_path / 'chips') + + stacks = chip_data.make_chip_stacks(data_chips, fmask_chips, label_chips, models.HLS_S30) + + assert len(stacks) == len(fmask_chips) + + for stack in stacks: + assert stack.id in stack.validation_mask.name + assert stack.id in stack.data.name + assert stack.id in stack.label.name + assert stack.modality == models.HLS_S30 + + +def test_filter_chips(warped_event_files, tmp_path): + fmask, label, data = warped_event_files + grid = chip_data.make_grid_from_reference(fmask) + + fmask_chips = chip_data.chip_data(grid, fmask, tmp_path / 'chips') + label_chips = chip_data.chip_data(grid, label, tmp_path / 'chips') + data_chips = chip_data.chip_data(grid, data, tmp_path / 'chips') + + stacks = chip_data.make_chip_stacks(data_chips, fmask_chips, label_chips, models.HLS_S30) + + filtered = chip_data.filter_chips(stacks) + + assert len(filtered) > 0 + assert len(filtered) < len(stacks) diff --git a/tests/test_chip_hyp3s1rtc.py b/tests/test_chip_hyp3s1rtc.py deleted file mode 100644 index 0216cab..0000000 --- a/tests/test_chip_hyp3s1rtc.py +++ /dev/null @@ -1,172 +0,0 @@ -import datetime -from pathlib import Path -from unittest.mock import MagicMock, patch - -import pytest -from shapely.geometry import box, mapping - -from satchip import chip_hyp3s1rtc, utils - - -def test_bounds_check(): - chip_hyp3s1rtc._check_bounds_size(utils.Bounds(0, 0, 1, 1)) - chip_hyp3s1rtc._check_bounds_size(utils.Bounds(0, 0, 2.9, 1)) - chip_hyp3s1rtc._check_bounds_size(utils.Bounds(-107.79192, 45.74287, -105.01543, 46.48598)) - - with pytest.raises(AssertionError): - chip_hyp3s1rtc._check_bounds_size(utils.Bounds(0, 0, 3, 1)) - - -def test_get_granules(): - bounds = utils.Bounds(-107.79192, 45.74287, -105.01543, 46.48598) - date_start = datetime.datetime(2020, 7, 7) - date_end = date_start + datetime.timedelta(days=14) - - mock_search_result = ['granule1', 'granule2'] - - with patch('satchip.chip_hyp3s1rtc.asf.geo_search', return_value=mock_search_result) as mock_geo_search: - results = chip_hyp3s1rtc._get_granules(bounds, date_start, date_end) - - mock_geo_search.assert_called_once() - - assert results == mock_search_result - - args, kwargs = mock_geo_search.call_args - assert ( - kwargs['intersectsWith'] - == 'POLYGON ((-105.01543 45.74287, -105.01543 46.48598, -107.79192 46.48598, -107.79192 45.74287, -105.01543 45.74287))' - ) - assert kwargs['start'] == date_start - assert kwargs['end'] == date_end + datetime.timedelta(days=1) - - -def test_get_slcs_for_each_chip_custom_intersect(): - granule1 = MagicMock() - granule1.geometry = mapping(box(0, 0, 2, 2)) - granule1.properties = {'startTime': '2025-01-01T00:00:00Z'} - - granule2 = MagicMock() - granule2.geometry = mapping(box(3, 3, 5, 5)) - granule2.properties = {'startTime': '2025-01-02T00:00:00Z'} - - granule3 = MagicMock() - granule3.geometry = mapping(box(10, 10, 15, 15)) - granule3.properties = {'startTime': '2025-01-03T00:00:00Z'} - - chip1 = MagicMock() - chip1.name = 'chip1' - chip1.bounds = [0, 0, 1, 1] - - chip2 = MagicMock() - chip2.name = 'chip2' - chip2.bounds = [1, 1, 2, 2] - - chip3 = MagicMock() - chip3.name = 'chip3' - chip3.bounds = [3, 3, 4, 4] - - chips = [chip1, chip2, chip3] - granules = [granule1, granule2, granule3] - - result = chip_hyp3s1rtc._get_slcs_for_each_chip(chips, granules, strategy='BEST') # type: ignore - - assert result['chip1'] == [granule1] - assert result['chip2'] == [granule1] - assert result['chip3'] == [granule2] - - -def test_get_slcs_for_each_chip_with_different_strategies(): - granule1 = MagicMock() - granule1.geometry = mapping(box(0, 0, 1, 1)) - granule1.properties = {'startTime': '2025-01-01T00:00:00Z'} - - granule2 = MagicMock() - granule2.geometry = mapping(box(0, 0, 5, 5)) - granule2.properties = {'startTime': '2025-01-02T00:00:00Z'} - - granule3 = MagicMock() - granule3.geometry = mapping(box(0, 0, 15, 15)) - granule3.properties = {'startTime': '2025-01-03T00:00:00Z'} - - chip1 = MagicMock() - chip1.name = 'chip1' - chip1.bounds = [0, 0, 5, 10] - - chips = [chip1] - granules = [granule1, granule2, granule3] - - result = chip_hyp3s1rtc._get_slcs_for_each_chip(chips, granules, strategy='BEST', intersection_pct=49) # type: ignore - assert result['chip1'] == [granule3] - - result = chip_hyp3s1rtc._get_slcs_for_each_chip(chips, granules, strategy='ALL', intersection_pct=49) # type: ignore - assert result['chip1'] == [granule3, granule2] - - -def test_get_slcs_for_each_chip_no_matches(): - chip = MagicMock() - chip.name = 'chip1' - chip.bounds = [0, 0, 1, 1] - - with pytest.raises(ValueError, match='No products found for chip chip1'): - chip_hyp3s1rtc._get_slcs_for_each_chip([chip], [], strategy='BEST') - - -class MockS1Product: - def __init__(self, scene_name: str): - self.properties = {'sceneName': scene_name} - - -def test_get_rtcs_for(): - slcs_for_chips = { - 'chip_001': [MockS1Product('SLC_1'), MockS1Product('SLC_2')], - 'chip_002': [MockS1Product('SLC_3'), MockS1Product('SLC_4')], - } - scratch_dir = Path('/tmp') - - mock_jobs = [] - for slc_name in ['SLC_1', 'SLC_2', 'SLC_3', 'SLC_4']: - job = MagicMock() - job.job_parameters = {'granules': [slc_name]} - mock_jobs.append(job) - - with ( - patch('satchip.chip_hyp3s1rtc._process_rtcs', return_value=mock_jobs) as mock_process_rtcs, - patch('satchip.chip_hyp3s1rtc._download_hyp3_rtc') as mock_download, - ): - - def mock_download_fn(job, scratch): - return { - 'VV': Path(f'/tmp/{job.job_parameters["granules"][0]}_rtc_VV.tif'), - 'VH': Path(f'/tmp/{job.job_parameters["granules"][0]}_rtc_VH.tif'), - } - - mock_download.side_effect = mock_download_fn - - result = chip_hyp3s1rtc._get_rtcs_for(slcs_for_chips, scratch_dir) - - expected = { - 'chip_001': [ - { - 'VV': Path('/tmp/SLC_1_rtc_VV.tif'), - 'VH': Path('/tmp/SLC_1_rtc_VH.tif'), - }, - { - 'VV': Path('/tmp/SLC_2_rtc_VV.tif'), - 'VH': Path('/tmp/SLC_2_rtc_VH.tif'), - }, - ], - 'chip_002': [ - { - 'VV': Path('/tmp/SLC_3_rtc_VV.tif'), - 'VH': Path('/tmp/SLC_3_rtc_VH.tif'), - }, - { - 'VV': Path('/tmp/SLC_4_rtc_VV.tif'), - 'VH': Path('/tmp/SLC_4_rtc_VH.tif'), - }, - ], - } - - assert result == expected - mock_process_rtcs.assert_called_once_with({'SLC_1', 'SLC_2', 'SLC_3', 'SLC_4'}) - assert mock_download.call_count == 4 diff --git a/tests/test_chip_sentinel2.py b/tests/test_chip_sentinel2.py deleted file mode 100644 index 32a81d4..0000000 --- a/tests/test_chip_sentinel2.py +++ /dev/null @@ -1,24 +0,0 @@ -from collections import namedtuple - -from satchip.chip_sentinel2 import get_latest_image_versions - - -ItemStub = namedtuple('ItemStub', ['id', 'properties']) - - -def test_get_latest_image_versions(): - items = [ - ItemStub(id='S2B_13TEG_20190623_0_L2A', properties={'s2:sequence': 0}), - ItemStub(id='S2B_13TEG_20190623_1_L2A', properties={'s2:sequence': 1}), - ItemStub(id='S2A_13TEG_20190621_0_L2A', properties={'s2:sequence': 0}), - ItemStub(id='S2A_13TEG_20190618_0_L2A', properties={'s2:sequence': 0}), - ItemStub(id='S2A_13TEG_20190618_1_L2A', properties={'s2:sequence': 1}), - ItemStub(id='S2A_13TEG_20190618_3_L2A', properties={'s2:sequence': 2}), - ] - - latest_items = get_latest_image_versions(items) # type: ignore - - assert len(latest_items) == 3 - assert any(item.id == 'S2B_13TEG_20190623_1_L2A' for item in latest_items) - assert any(item.id == 'S2A_13TEG_20190621_0_L2A' for item in latest_items) - assert any(item.id == 'S2A_13TEG_20190618_3_L2A' for item in latest_items) diff --git a/tests/test_download_data.py b/tests/test_download_data.py new file mode 100644 index 0000000..f1267e0 --- /dev/null +++ b/tests/test_download_data.py @@ -0,0 +1,12 @@ +import pytest + +from satchip import download_data, models + + +@pytest.mark.download +def test_download_data(s2_event, tmp_path): + local_files = download_data.download_data(s2_event, models.HLS_S30, tmp_path) + assert len(local_files) > 1 + + local_files = download_data.download_data(s2_event, models.HLS_L30, tmp_path) + assert len(local_files) == 0 diff --git a/tests/test_generate_labels.py b/tests/test_generate_labels.py new file mode 100644 index 0000000..95be506 --- /dev/null +++ b/tests/test_generate_labels.py @@ -0,0 +1,15 @@ +import rasterio + +from satchip import generate_labels + + +def test_generate_labales(s2_local_files, s2_event, tmp_path): + template_file = s2_local_files[0] + + output = generate_labels.binary_mask_from_template(template_file, s2_event, tmp_path / 'labels') + + assert output.name.endswith('MASK.tif') + + with rasterio.open(output) as mask_ds: + with rasterio.open(template_file) as template_ds: + assert mask_ds.shape == template_ds.shape diff --git a/tests/test_integration.py b/tests/test_integration.py deleted file mode 100644 index dba034c..0000000 --- a/tests/test_integration.py +++ /dev/null @@ -1,59 +0,0 @@ -from datetime import datetime -from pathlib import Path - -import pytest -from osgeo import gdal - -from satchip.chip_data import create_chips -from satchip.chip_label import chip_labels - - -gdal.UseExceptions() - - -def create_dataset(outpath: Path, start: tuple[int, int]) -> Path: - x, y = start - pixel_size = 10 - cols, rows = 512, 512 - driver = gdal.GetDriverByName('GTiff') - dataset = driver.Create(str(outpath), cols, rows, 1, gdal.GDT_UInt16) - dataset.SetGeoTransform((x, pixel_size, 0, y, 0, -pixel_size)) - dataset.SetProjection('EPSG:32611') - array = dataset.GetRasterBand(1).ReadAsArray() - array[:, :] = 0 - array[128:384, 128:384] = 1 - dataset.GetRasterBand(1).WriteArray(array) - dataset.FlushCache() - dataset = None - return outpath - - -def create_label_and_data(label_tif, out_dir, image_dir): - chip_labels(label_tif, datetime.fromisoformat('20240115'), out_dir) - for platform in ['S2L2A', 'HLS', 'S1RTC']: - create_chips( - list((out_dir / 'LABEL').glob('*.zarr.zip')), - platform, - datetime.fromisoformat('20240101'), - datetime.fromisoformat('20240215'), - 'BEST', - 20, - out_dir, - image_dir, - ) - - -@pytest.mark.integration -def test_integration(): - data_dir = Path('integration_test') - train_dir = data_dir / 'train' - train_dir.mkdir(parents=True, exist_ok=True) - val_dir = data_dir / 'val' - val_dir.mkdir(parents=True, exist_ok=True) - image_dir = data_dir / 'images' - image_dir.mkdir(parents=True, exist_ok=True) - - train_tif = create_dataset(data_dir / 'train.tif', (431795, 3943142)) - create_label_and_data(train_tif, train_dir, image_dir) - val_tif = create_dataset(data_dir / 'val.tif', (431795, 3943142 - 10 * 512)) - create_label_and_data(val_tif, val_dir, image_dir) diff --git a/tests/test_merge_modality.py b/tests/test_merge_modality.py new file mode 100644 index 0000000..68b0ac5 --- /dev/null +++ b/tests/test_merge_modality.py @@ -0,0 +1,95 @@ +from pathlib import Path +import rasterio + +from satchip import generate_labels, merge_modality, models, download_data + + +def test_multi_projection_event(pristine_gdf): + swath = pristine_gdf.iloc[2] + + data_path = Path(__file__).parent / 'data' + raw_path = data_path / 'raw' + merge_path = Path(__file__).parent / 'data' / 'merge' + + modality = models.HLS_S30 + + event = models.Event(name=swath['HLSID'], date=swath['SwathDate'], wgs84_geometry=swath['geometry'], buffer_m=10000) + + local_files = download_data.download_data(event, modality, raw_path) + merged_event = merge_modality.merge_modality(local_files, modality, event=event, output_path=merge_path) + + +def test_make_merge_name(s2_event): + merged_name = merge_modality._make_merge_name(s2_event.name, s2_event.date, 'BAND', 'MOD') + + assert 'BAND' in merged_name + assert 'MOD' in merged_name + assert s2_event.name in merged_name + assert merged_name.endswith('.tif') + + +def test_merge_s2_modality_empty(s2_local_files, s2_event, tmp_path): + empty_result = merge_modality.merge_modality([], models.HLS_S30, s2_event, tmp_path) + assert len(empty_result) == 0 + + +def test_merge_s2_modality(s2_local_files, s2_event, tmp_path): + merged_files = merge_modality.merge_modality(s2_local_files, models.HLS_S30, s2_event, tmp_path / 'merged') + + assert len(merged_files) == len(models.HLS_S30['bands']) + assert all('merged' in str(result) for result in merged_files) + + shapes = set() + for merged_file in merged_files: + with rasterio.open(merged_file) as ds: + shapes.add(ds.shape) + + assert len(shapes) == 1 + + +def test_stack_bands(s2_local_files, s2_event, tmp_path): + bands = tuple(b for b in models.HLS_S30['bands'] if b.id != 'Fmask') + + merged_bands = merge_modality.merge_modality(s2_local_files, models.HLS_S30, s2_event, tmp_path / 'merged', bands) + stacked_data = merge_modality.stack_bands(merged_bands, stacked_filename=tmp_path / 'stacked.tif') + + assert 'stacked.tif' in stacked_data.name + + with rasterio.open(stacked_data) as ds: + num_bands = ds.count + + assert num_bands == len(bands) + + +def test_reproject_files(s2_local_files, s2_event, tmp_path): + reprojected_files = merge_modality.reproject_files(s2_local_files, tmp_path / 'wgs84') + + assert len(reprojected_files) == len(s2_local_files) + + for f in reprojected_files: + with rasterio.open(f) as ds: + epsg_code = ds.profile['crs'].to_epsg() + assert epsg_code == 4326 + + +def test_align_to_reference(s2_local_files, s2_event, tmp_path): + merged = merge_modality.merge_modality(s2_local_files, models.HLS_S30, s2_event, tmp_path / 'merged') + reprojected = merge_modality.reproject_files(merged, tmp_path / 'wgs84') + + bands = reprojected[:-1] + fmask = reprojected[-1] + + stacked = merge_modality.stack_bands(bands, stacked_filename=tmp_path / 'stacked.tif') + label = generate_labels.binary_mask_from_template(stacked, s2_event, tmp_path) + + outputs = merge_modality.align_to_reference( + label, [stacked, fmask], tmp_path / 'aligned', s2_event.buffered_geometry().bounds + ) + + shapes = set() + for output in outputs: + with rasterio.open(output) as ds: + shapes.add(ds.shape) + + assert len(outputs) == 3 + assert len(shapes) == 1 diff --git a/tests/test_models.py b/tests/test_models.py new file mode 100644 index 0000000..49a2d83 --- /dev/null +++ b/tests/test_models.py @@ -0,0 +1,55 @@ +import pytest + +from satchip import models + + +def test_hls_model_len(): + assert len(models.HLS_S30_BANDS) == 7 + assert len(models.HLS_L30_BANDS) == 7 + + +def test_hls_model_names(): + assert models.HLS_S30_BANDS['N'].id != models.HLS_L30_BANDS['N'].id + assert models.HLS_S30_BANDS['SW1'].id != models.HLS_L30_BANDS['SW1'].id + assert models.HLS_S30_BANDS['SW2'].id != models.HLS_L30_BANDS['SW2'].id + + +@pytest.mark.parametrize( + 'filename, expected_band', + [ + ('OPERA_L2_RTC-S1_T085-181260-IW1_20190822T125312Z_20250913T222203Z_S1A_30_v1.0_VH.tif', 'VH'), + ('OPERA_L2_RTC-S1_T165-352512-IW3_20200723T000516Z_20250908T213809Z_S1B_30_v1.0_VH.tif', 'VH'), + ('OPERA_L2_RTC-S1_T085-181260-IW1_20190822T125312Z_20250913T222203Z_S1A_30_v1.0_VV.tif', 'VV'), + ('OPERA_L2_RTC-S1_T165-352512-IW3_20200723T000516Z_20250908T213809Z_S1B_30_v1.0_VV.tif', 'VV'), + ('OPERA_L2_RTC-S1_T085-181260-IW1_20190822T125312Z_20250913T222203Z_S1A_30_v1.0_mask.tif', 'mask'), + ('OPERA_L2_RTC-S1_T165-352512-IW3_20200723T000516Z_20250908T213809Z_S1B_30_v1.0_mask.tif', 'mask'), + ], +) +def test_band_id_from_filename_rtc(filename, expected_band): + assert models.band_id_from_filename(filename, 'OPERA_RTC').id == expected_band + + +@pytest.mark.parametrize( + 'filename, expected_band', + [ + ('HLS.S30.T13TFL.2019187T174919.v2.0.B12.tif', 'B12'), + ('HLS.S30.T13TGM.2018184T173901.v2.0.B11.tif', 'B11'), + ('HLS.S30.T14TLQ.2019219T173911.v2.0.B8A.tif', 'B8A'), + ('HLS.S30.T15TXG.2017167T170311.v2.0.B04.tif', 'B04'), + ('HLS.S30.T13TFL.2019187T174919.v2.0.Fmask.tif', 'Fmask'), + ('HLS.S30.T13TGM.2018184T173901.v2.0.B12.tif', 'B12'), + ('HLS.S30.T14TLQ.2019219T173911.v2.0.B11.tif', 'B11'), + ('HLS.S30.T15TXG.2017167T170311.v2.0.B8A.tif', 'B8A'), + ('HLS.S30.T13TFL.2020157T174911.v2.0.B02.tif', 'B02'), + ('HLS.S30.T13TGM.2018184T173901.v2.0.Fmask.tif', 'Fmask'), + ('HLS.S30.T13TFL.2020157T174911.v2.0.B03.tif', 'B03'), + ('HLS.S30.T14SKH.2019211T172909.v2.0.B02.tif', 'B02'), + ('HLS.S30.T14TLQ.2019219T173911.v2.0.Fmask.tif', 'Fmask'), + ('HLS.S30.T13TFL.2020157T174911.v2.0.B04.tif', 'B04'), + ('HLS.S30.T14SKH.2019211T172909.v2.0.B03.tif', 'B03'), + ('HLS.S30.T14TLQ.2020156T172859.v2.0.B02.tif', 'B02'), + ('HLS.S30.T15TXG.2017167T170311.v2.0.Fmask.tif', 'Fmask'), + ], +) +def test_band_id_from_filename_hls_s30(filename, expected_band): + assert models.band_id_from_filename(filename, 'HLS_S30').id == expected_band diff --git a/tests/test_stub.py b/tests/test_stub.py deleted file mode 100644 index 31ba320..0000000 --- a/tests/test_stub.py +++ /dev/null @@ -1,2 +0,0 @@ -def test_stub(): - assert True diff --git a/tests/test_terra_mind_grid.py b/tests/test_terra_mind_grid.py deleted file mode 100644 index 93b4d15..0000000 --- a/tests/test_terra_mind_grid.py +++ /dev/null @@ -1,23 +0,0 @@ -import numpy as np -import pytest - -from satchip.terra_mind_grid import TerraMindGrid - - -@pytest.mark.parametrize( - 'tm_name, point', - [ - ('374U_897R_3_2', [96.80658, 33.62825]), - ('735U_418L_0_1', [-92.35263, 66.07510]), - ('611U_462L_2_1', [-72.02323, 54.93515]), - ], -) -def test_terra_mind_grid(tm_name, point): - mt_name = f'{tm_name.split("_")[0]}_{tm_name.split("_")[1]}' - tmp_grid = TerraMindGrid((np.floor(point[1]), np.ceil(point[1])), (np.floor(point[0]), np.ceil(point[0]))) - # mt_chip = [x for x in tmp_grid.major_tom_chips if x.name == mt_name][0] - chips = [x for x in tmp_grid.terra_mind_chips if x.name.startswith(mt_name)] - in_lon = [x for x in chips if x.bounds[0] < point[0] < x.bounds[2]] - in_lon_lat = [x for x in in_lon if x.bounds[1] < point[1] < x.bounds[3]] - assert len(in_lon_lat) == 1 - assert in_lon_lat[0].name == tm_name diff --git a/tests/test_view.py b/tests/test_view.py new file mode 100644 index 0000000..beec136 --- /dev/null +++ b/tests/test_view.py @@ -0,0 +1,41 @@ +from satchip import view, merge_modality, chip_data, models + + +def test_view_merged(s2_local_files, s2_event, tmp_path): + bands = tuple(b for b in models.HLS_S30['bands'] if b.id != 'Fmask') + + merged_bands = merge_modality.merge_modality(s2_local_files, models.HLS_S30, s2_event, tmp_path / 'merged', bands) + stacked_data = merge_modality.stack_bands(merged_bands, stacked_filename=tmp_path / 'stacked.tif') + reprojected = merge_modality.reproject_files([stacked_data], tmp_path / 'reproj')[0] + + view.view_merged( + reprojected, s2_event, models.HLS_S30, rgb_bands=[2, 1, 0], save_to_file=tmp_path / 'output.png', quite=False + ) + + +def test_view_chip(warped_event_files, tmp_path): + fmask, label, data = warped_event_files + grid = chip_data.make_grid_from_reference(fmask) + + fmask_chips = chip_data.chip_data(grid, fmask, tmp_path / 'chips') + label_chips = chip_data.chip_data(grid, label, tmp_path / 'chips') + data_chips = chip_data.chip_data(grid, data, tmp_path / 'chips') + + stacks = chip_data.make_chip_stacks(data_chips, fmask_chips, label_chips, models.HLS_S30) + filtered = chip_data.filter_chips(stacks) + + view.view_chip(filtered[0], models.HLS_S30, [2, 1, 0]) + + +def test_view_chips(warped_event_files, tmp_path): + fmask, label, data = warped_event_files + grid = chip_data.make_grid_from_reference(fmask) + + fmask_chips = chip_data.chip_data(grid, fmask, tmp_path / 'chips') + label_chips = chip_data.chip_data(grid, label, tmp_path / 'chips') + data_chips = chip_data.chip_data(grid, data, tmp_path / 'chips') + + stacks = chip_data.make_chip_stacks(data_chips, fmask_chips, label_chips, models.HLS_S30) + filtered = chip_data.filter_chips(stacks) + + view.view_chips(filtered, models.HLS_S30, [2, 1, 0])