basketballcantho/remote-sensing
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

thêm chức năng train trên odc predict trên planetary

Browse Source
This commit is contained in:
Victor Phan
2026-03-04 23:03:04 +07:00
parent ebb8e6e4b3
commit 8a1e7bb22e
39 changed files with 15297 additions and 711 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/data_products/novasar-1-gamma0.ipynb
+952
View File
@@ -0,0 +1,952 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "dae926b2-fa06-4f00-b18d-7f74e92ce676",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"## NovaSAR-1 data products <img align=\"right\" src=\"../../resources/csiro_easi_logo.png\">\n",
"\n",
"#### Index\n",
"- [Overview](#Overview)\n",
" - [Tips and tricks](#Tips-and-tricks)\n",
" - [Background on SAR data](#Background-on-SAR-data)\n",
"- [Setup (imports, defaults, dask, odc)](#Setup)\n",
"- [Example query](#Example-query)\n",
" - [Review scene metadata for the query](#Review-scene-metadata-for-the-query)\n",
" - [Filter scenes based on metadata only](#Filter-scenes-based-on-metadata-only)\n",
"- [Load data into a virtual dask array](#Load-data-into-a-virtual-dask-array)\n",
"- [Conversion and helper functions](#Conversion-and-helper-functions)\n",
" - [Filter scenes based on valid data](#Filter-scenes-based-on-valid-data)\n",
" - [Add dB and amplitude values to the dataset](#Add-dB-and-amplitude-values-to-the-dataset)\n",
"- [Persist the data into the dask workers](#Persist-the-data-into-the-dask-workers)\n",
" - [Plot the data](#Plot-the-data)\n",
" - [Plot histograms of the dB data](#Plot-histograms-of-the-dB-data)\n",
" - [Make an RGB image](#Make-an-RGB-image)\n",
"- [Export to Geotiffs](#Export-to-Geotiffs)"
]
},
{
"cell_type": "markdown",
"id": "23474415-05a4-4065-b28a-3acbe7129f7f",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"This notebook demonstrates how to load and use NovaSAR-1 ARD _gamma-0 backscatter_ data products from the [CSIRO NovaSAR Facility](https://research.csiro.au/cceo/novasar/).\n",
"\n",
"The [NovaSAR-1 data products](https://research.csiro.au/cceo/novasar/about/novasar-1-user-guide/#products) include:\n",
"- Level-1 product types Multi-Look Detected (GRD, SCD, SRD) and Single Look Complex (SLC)\n",
"- Level-2 analysis ready data (ARD) gamma-0 radiometric terrain corrected backscatter\n",
"\n",
"The **Level-1** data products are stored in swath (line, row) raster grids, with a \"fake\" WGS-84 bounding box grid for ODC indexing. This means they are not suited to `datacube.load()` directly. Instructions TBC.\n",
"\n",
"The **Level-2 ARD** data products are remapped to WGS-84 grids (various resolutions) as part of the gamma-0 backscatter processing. These data can be opened directly with `datacube.load()`."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "35b91cdb-2c5a-48c0-a0bf-7de2c2637bc5",
"metadata": {},
"source": [
"## Tips and tricks\n",
"\n",
"#### 1. NovaSAR-1 product names\n",
"\n",
"The ODC currently requires a “dataset” to contain all bands that are defined for its parent product. In practice, each of the NovaSAR-1 product types (GRD, SCD, SRD, SLC, ARD) can have different combinations of available polarization bands across all scene acquisitions. This means we need multiple ODC products; one for each combination of polarizations.\n",
"\n",
"> WIP: relax this constraint in ODC so that we can have \"optional\" bands defined, which will allow us to aggregate polarization combinations for a single product type.\n",
"\n",
"View the coverage of indexed NovaSAR-1 products and datasets in the [CSIRO EASI Explorer \"SAR products group\"](https://explorer.csiro.easi-eo.solutions/products#synthetic-aperture-radar-sar-group).\n",
"\n",
"#### 2. Load to target grid or crs/resolution\n",
"\n",
"Different scenes (ODC datasets) in a product may have different spatial resolutions depending on the acquisition mode, so its usually appropriate to provide a target grid. For example, use either\n",
"\n",
"- `datacube.load(..., like=target_geobox)`. See [odc-geo](https://github.com/opendatacube/odc-geo) to help create a target \"geobox\", or use an existing (ODC-compatible) xarray object\n",
"- `datacube.load(..., output_crs=target_crs, resolution=target_res, align=target_res/2)`\n",
"\n",
"For convenience and general applicability, default load parameters are defined for all NovaSAR-1 L2 ARD products as:\n",
"\n",
"- `CRS = \"epsg:4326\" (WGS-84)`, `Pixel_size = 0.0002 deg` (20 m) and `Align = 0.0001` (10 m). \n",
"\n",
"This means a `datacube.load()` call without specific output grid parameters (no `like`, `output_crs`, `resolution`) will load data into a _WGS-84, 20 m, aligned to pixel centres_ grid with extents given by the `latitude/longitude/x/y/crs` parameters.\n",
"\n",
"> Native NovaSAR-1 L2 ARD products are aligned to pixel centres (AREA_OR_POINT=Point). You should ensure this is taken into account if mapping pixels to a target grid that is approximately the same resolution as the native scene resolution (otherwise you may be shifting by half-pixel). If mapping to a coarser resolution if may not matter as much.\n",
"\n",
"#### 3. Select or filter scenes before loading\n",
"\n",
"Scenes can be pre-filtered using the scene metadata and the `datacube.find_datasets()` function before datasets are loaded into an xarray object. There are two easy (and similar) ways to achieve this:\n",
"\n",
"1. Use `datacube.find_datasets()` with _product, time, space_ search parameters. This returns a list of dataset items that can be filtered manually, e.g. by scene name or scene metadata. Pass your filtered list to `datacube.load(datasets=my_datasets_list, like=target_geobox)`.\n",
" - This is demonstrated below by parsing likely useful metadata for a list of datasets items into a pandas table, for viewing.\n",
"1. Create a metadata predicate function that returns `True` if a dataset's metadata satisifies the conditions specified in your function. This predicate function can be passed to `datacube.load(..., dataset_predicate=my_filter_function`), which is then applied internally in a call to `find_datasets`. See [datacube.load() help](https://opendatacube.readthedocs.io/en/latest/api/indexed-data/generate/datacube.Datacube.load.html) for an example.\n",
"\n",
"#### 4. Units and conversions\n",
"The `novasar_l2ard_*` data are given in _Digital numbers_ (DN). DNs can be converted to _decibel (dB)_ or _linear amplitude_, and vice-versa, with the following equations. Practical _Xarray_ examples are given below.\n",
"\n",
"Amplitude to/from dB:\n",
"```\n",
"dB = 20 * log10(DN) + K\n",
"DN = 10^((dB-K)/20)\n",
"\n",
"where K is a calibration factor, which for NovaSAR-1 is -83 dB.\n",
"```\n",
"\n",
"Digital numbers to/from Amplitude:\n",
"```\n",
"amplitude = DN / 14125.3754\n",
"DN = amplitude * 14125.3754\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "5e5120c6-e943-46a9-a4d8-c5b759586e62",
"metadata": {},
"source": [
"## Background on SAR data\n",
"\n",
"An excellent introduction and overview to using SAR data is provided in the [CEOS Laymans SAR Interpretation Guide](https://ceos.org/ard/files/Laymans_SAR_Interpretation_Guide_3.0.pdf). This guide has also been converted to a set of set of Jupyter notebooks that you can download from https://github.com/AMA-Labs/cal-notebooks/tree/main/examples/SAR.\n",
"\n",
"The SAR instrument on the NovaSAR-1 satellite operates in the ***S-band at approximately 9.4 cm wavelength***. This means that it can \"see\" objects of about this size and larger, and smaller objects are relatively transparent. Compared to *Sentinel-1 (C-band, 5.6 cm)*, NovaSAR-1 S-band is better able to penetrate through vegetation.\n",
"\n",
"> The SAR signal responds to the orientation and scattering from surface features of comparable size or larger than the wavelength.\n",
"> - A bright backscatter value typically means the surface was orientated perpendicular to the signal incidence angle and most of the signal was reflected back to the satellite (direct backscatter)\n",
"> - A dark backscatter value means most of the signal was reflected away from the satellite (forward scattering) and typically responds to a smooth surface (relative to the wavelength) such as calm water or bare soil\n",
"> - Rough surfaces (relative to the wavelength) result in diffuse scattering where some of the signal is returned to the satellite.\n",
"> - Complex surfaces may result in volume scattering (scattering within a tree canopy) or double-bounce scattering (perpendicular objects such as buildings and structures)\n",
"> - The relative backscatter values of co-polarisation (HH, VV) and cross-polarisation (HV) measurements can provide information on the scattering characteristics of the surface features.\n",
"\n",
"Using NovaSAR-1 backscatter data requires interpretation of the data for different surface features, including as these features change spatially or in time. It may also be necessary to carefully consider the incidence angle of the SAR signal relative to the surface features using the *incidence_angle* band or the satellite direction metadata (descending = north to south; ascending = south to north)."
]
},
{
"cell_type": "markdown",
"id": "dd5c0ff4-e1a2-4adf-98d1-1ca5b9f3126d",
"metadata": {},
"source": [
"## Set up"
]
},
{
"cell_type": "markdown",
"id": "4cb7c012-719e-4f13-a206-36f271c6a91c",
"metadata": {},
"source": [
"#### Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c848c17-2e7c-4b68-9682-b8994644f521",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Common imports and settings\n",
"import os, sys, re\n",
"from pathlib import Path\n",
"from IPython.display import Markdown\n",
"import pandas as pd\n",
"pd.set_option(\"display.max_rows\", None)\n",
"import xarray as xr\n",
"import numpy as np\n",
"\n",
"# Datacube\n",
"import datacube\n",
"from datacube.utils.aws import configure_s3_access\n",
"import odc.geo.xr # https://github.com/opendatacube/odc-geo\n",
"from datacube.utils import masking # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py\n",
"from dea_tools.plotting import display_map # https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools\n",
"\n",
"# Dask\n",
"import dask\n",
"from dask.distributed import Client, LocalCluster\n",
"\n",
"# Basic plots\n",
"%matplotlib inline\n",
"# import matplotlib.pyplot as plt\n",
"# plt.rcParams['figure.figsize'] = [12, 8]\n",
"\n",
"# Holoviews\n",
"# https://holoviz.org/tutorial/Composing_Plots.html\n",
"# https://holoviews.org/user_guide/Composing_Elements.html\n",
"import hvplot.xarray\n",
"import holoviews as hv\n",
"import panel as pn\n",
"from datashader import reductions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bde8313b-887d-4529-bd57-89f740ff6380",
"metadata": {},
"outputs": [],
"source": [
"# EASI defaults and convenience functions\n",
"# easi-tools is in the repo https://github.com/csiro-easi/easi-notebooks\n",
"# These are convenience functions for using these notebooks in EASI; comment-out if not required\n",
"\n",
"repo = Path.home() / 'easi-notebooks' # Change path as necessary\n",
"if str(repo) not in sys.path:\n",
" sys.path.append(str(repo))\n",
"\n",
"from easi_tools import initialize_dask, xarray_object_size # Useful"
]
},
{
"cell_type": "markdown",
"id": "bb66366a-f2d5-461a-a34f-7eed5a5296fb",
"metadata": {
"tags": []
},
"source": [
"#### Dask cluster\n",
"\n",
"Using a local _Dask_ cluster is a good habit to get into. It can simplify loading and processing of data in many cases, and it provides a dashboard that shows the loading/processing progress.\n",
"\n",
"To learn more about _Dask_ see the set of [dask notebooks](https://github.com/csiro-easi/easi-notebooks/tree/main/html#dask-tutorials)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "02f9e778-755c-4ff2-a99f-ad5c48438665",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Local cluster - good for small exploratory and testing work\n",
"cluster, client = initialize_dask(workers=8)\n",
"display(client)\n",
"\n",
"# Or use Dask Gateway - this may take a few minutes\n",
"# cluster, client = initialize_dask(use_gateway=True, workers=4)\n",
"# display(client)"
]
},
{
"cell_type": "markdown",
"id": "a746582e-be79-448c-a9c4-7a2ae2d703ce",
"metadata": {},
"source": [
"#### ODC database\n",
"\n",
"Connect to the ODC database. Configure the environment and low-level tools to read from AWS buckets."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0e2232ee-f0a5-43fc-a6c6-d422b437cd69",
"metadata": {},
"outputs": [],
"source": [
"dc = datacube.Datacube()\n",
"\n",
"# Access AWS \"requester-pays\" buckets\n",
"# This is necessary for reading data from most third-party AWS S3 buckets such as for Landsat and Sentinel-2\n",
"configure_s3_access(aws_unsigned=False, requester_pays=True, client=client);"
]
},
{
"cell_type": "markdown",
"id": "1801db2b-e24e-45fa-85d4-e153ebcc8f6f",
"metadata": {},
"source": [
"## Example query\n",
"\n",
"Change any of the parameters in the `query` object below to adjust the location, time, projection, or spatial resolution of the returned datasets.\n",
"\n",
"Use the Explorer interface to check the temporal and spatial coverage for each product."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a04d875-20d4-4d72-a16d-e76b7f48161c",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Set your own latitude / longitude\n",
"\n",
"# Menindee Lakes, Australia\n",
"latitude_range = (-32.1, -33.6)\n",
"longitude_range = (141.5, 143)\n",
"\n",
"# Great Western Woodlands, Australia\n",
"# latitude_range = (-33, -32.6)\n",
"# longitude_range = (120.5, 121)\n",
"\n",
"time_range = None # All available times\n",
"\n",
"available_ard_products = [\n",
" 'novasar_l2ard_hh', # Moderate coverage, https://explorer.csiro.easi-eo.solutions/products/novasar_l2ard_hh\n",
" 'novasar_l2ard_hh_hv', # High coverage, https://explorer.csiro.easi-eo.solutions/products/novasar_l2ard_hh_hv\n",
" 'novasar_l2ard_hv', # Low coverage, https://explorer.csiro.easi-eo.solutions/products/novasar_l2ard_hv\n",
" 'novasar_l2ard_vv', # Moderate coverage, https://explorer.csiro.easi-eo.solutions/products/novasar_l2ard_vv\n",
" 'novasar_l2ard_vv_hh', # Moderate coverage, https://explorer.csiro.easi-eo.solutions/products/novasar_l2ard_vv_hh\n",
" 'novasar_l2ard_vv_hh_hv', # High coverage, https://explorer.csiro.easi-eo.solutions/products/novasar_l2ard_vv_hh_hv\n",
"]\n",
"\n",
"# Select a product\n",
"product_name = 'novasar_l2ard_hh_hv'\n",
"polarizations = re.findall('_([vh]+)', product_name) # The set of polarization bands for the selected product\n",
"\n",
"query = {\n",
" 'product': product_name, # Product name\n",
" 'measurements': polarizations + ['angle', 'mask', 'scatteringarea', 'gammatosigmaratio'], # All bands for testing\n",
" 'x': longitude_range, # \"x\" axis bounds\n",
" 'y': latitude_range, # \"y\" axis bounds\n",
" 'time': time_range, # Any parsable date strings\n",
"}\n",
"\n",
"# Convenience function to display the selected area of interest\n",
"display_map(longitude_range, latitude_range)"
]
},
{
"cell_type": "markdown",
"id": "b399fe85-ba78-4a8e-8a33-3dcdfa95b7f2",
"metadata": {},
"source": [
"## Review scene metadata for the query\n",
"\n",
"Run `datacube.find_datasets()` for the query and summarise resulting scene information in a pandas table"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4d84334a-1749-412c-b10c-e7bbc53c43d8",
"metadata": {},
"outputs": [],
"source": [
"datasets = dc.find_datasets(**query)\n",
"\n",
"# Examples\n",
"# datasets[0].metadata_doc --> dict of all metadata\n",
"# datasets[0].measurements --> dict of measurements\n",
"# datasets[0].grids --> dict of grids (shape and transform) defined for measurements\n",
"# datasets[0].crs --> CRS object\n",
"# datasets[0].extent --> geometry object\n",
"# datasets[0].time --> start and end time objects\n",
"# datasets[0].properties --> dict of properties\n",
"\n",
"# Select your own fields\n",
"fields = {\n",
" 'name': None,\n",
" 'time': None,\n",
" 'gsd': None,\n",
" 'polarizations': None,\n",
" 'pass_direction': None,\n",
" 'antenna_pointing': None,\n",
" 'platform_heading': None,\n",
" 'incident_angle_near_range': None,\n",
" 'incident_angle_far_range': None,\n",
" 'source_product_level': None,\n",
" 'operational_mode_name': None,\n",
" 'filter_applied': None,\n",
" 'noise_removal_applied': None,\n",
"}\n",
"\n",
"# Parse into a pandas DataFrame\n",
"df = []\n",
"for ds in datasets:\n",
" ff = fields.copy()\n",
" # Non-properties fields\n",
" ff['name'] = ds.metadata.label\n",
" ff['time'] = ds.time[0] # time->range(start, end) or center_time->(end-start)/2\n",
" # Properties fields\n",
" for xx in fields:\n",
" if xx in ('name', 'time'):\n",
" continue\n",
" pre = 'eo' if xx == 'gsd' else 'novasar'\n",
" ff[xx] = ds.properties[f\"{pre}:{xx}\"]\n",
" df.append(ff)\n",
"df = pd.DataFrame(df)\n",
"display(df.sort_values('time') if not df.empty else Markdown(\"**Dataframe is empty**. Choose different product, space, time parameters\"))"
]
},
{
"cell_type": "markdown",
"id": "1f263ddb-f156-4841-8a01-9ced702cc820",
"metadata": {},
"source": [
"## Filter scenes based on metadata only\n",
"\n",
"This step is optional and configurable. Here we show a selection of simple filters:\n",
"\n",
"- Select descending passes only\n",
"- Select the N most recent scenes"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "93b759fa-21d1-4d4b-929a-38927e6e0448",
"metadata": {},
"outputs": [],
"source": [
"selected_datasets = datasets # Start with the full list\n",
"\n",
"# Select descending passes only\n",
"selected_datasets = [\n",
" xx for xx in selected_datasets if xx.properties[\"novasar:pass_direction\"] == \"DESCENDING\"\n",
"]\n",
"\n",
"# Select the N most recent scenes\n",
"selected_datasets = sorted(selected_datasets, key=lambda ds: ds.time[0])[-10:]\n",
"\n",
"display(Markdown(f\"**Number of selected scenes**: {len(selected_datasets)}\"))\n",
"selected_datasets"
]
},
{
"cell_type": "markdown",
"id": "9fefeef8-0c21-49c2-b86f-c68e5178a7ac",
"metadata": {
"tags": []
},
"source": [
"## Load data into a virtual dask array"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2f4c724f-636a-4833-a165-adb053acaa07",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Target xarray parameters\n",
"# - Target GeoBox or output CRS and resolution\n",
"# - Usually we group input scenes on the same day to a single time layer (groupby)\n",
"# - Select a reasonable Dask chunk size. This should be adjusted depending on the spatial extents and target grid parameters you choose\n",
"\n",
"# Default target grid is WGS-84, 20 m pixels for the given spatial extents\n",
"# - Create or reuse a GeoBox for precise mapping of input pixels to a target grid\n",
"# - Or provide output_crs and resolution, which will create best-fit GeoBox\n",
"\n",
"# Example similar Geobox constructor\n",
"# geobox = odc.geo.geobox.GeoBox.from_bbox(\n",
"# [query['x'][0], query['y'][1], query['x'][1], query['y'][0]],\n",
"# crs='epsg:4326', resolution=0.0002, anchor=0.0001\n",
"# )\n",
"\n",
"load_params = {\n",
" 'datasets': selected_datasets,\n",
" 'group_by': 'solar_day', # Scene grouping\n",
" 'dask_chunks': {'latitude':2048, 'longitude':2048}, # Dask chunks\n",
"}\n",
"\n",
"# Load data\n",
"data = dc.load(**(query | load_params))\n",
"\n",
"display(xarray_object_size(data))\n",
"display(data)"
]
},
{
"cell_type": "markdown",
"id": "e3d86506-39ba-4f02-ac6e-2d1d27119668",
"metadata": {},
"source": [
"## Conversion and helper functions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "99eb8a2c-fe00-44b7-bd39-250ab3d323cc",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# These functions use numpy, which should be satisfactory for most notebooks.\n",
"# Calculations for larger or more complex arrays may require Xarray's \"ufunc\" capability.\n",
"# https://docs.xarray.dev/en/stable/examples/apply_ufunc_vectorize_1d.html\n",
"#\n",
"# Apply numpy.log10 to the DataArray\n",
"# log10_data = xr.apply_ufunc(np.log10, data)\n",
"\n",
"def dn_to_decibel(da: 'xr.DataArray', K=0):\n",
" \"\"\"Return an array converted to dB values\"\"\"\n",
" xx = da.where(da > 0, np.nan) # Set values <= 0 to NaN\n",
" xx = 20*np.log10(xx) + K\n",
" xx.attrs.update({\"units\": \"dB\"})\n",
" return xx\n",
"\n",
"def decibel_to_dn(da: 'xr.DataArray', K=0):\n",
" \"\"\"Return an array converted to digital number values\"\"\"\n",
" xx = np.power(10, (da-K)/20.0)\n",
" xx.attrs.update({\"units\": \"DN\"})\n",
" return xx\n",
"\n",
"def dn_to_amplitude(da: 'xr.DataArray'):\n",
" \"\"\"Return an array converted to linear amplitude values\"\"\"\n",
" scale_factor = da.attrs['scale_factor'] if 'scale_factor' in da.attrs else 1\n",
" # print(f\"DN to amplitude scale_factor: 1/{1/scale_factor} = {scale_factor}\")\n",
" xx = da.where(da > 0, np.nan) # Set values <= 0 to NaN\n",
" xx = xx * scale_factor\n",
" xx.attrs.update({\"units\": \"amplitude\"})\n",
" return xx\n",
"\n",
"def select_valid_time_layers(ds: 'xarray', percent: float = 5):\n",
" \"\"\"Select time layers that have at least a given percentage of valid data (e.g., >=5%)\n",
"\n",
" Example usage:\n",
" selected = select_valid_time_layers(ds, percent=5)\n",
" filtered == ds.sel(time=selected)\n",
" \"\"\"\n",
" spatial_dims = ds.odc.spatial_dims\n",
" nelements = ds.sizes[spatial_dims[0]] * ds.sizes[spatial_dims[1]]\n",
" if ds.dtype =='bool':\n",
" return ds.sum(dim=spatial_dims).values / nelements >= (percent/100.0)\n",
" return ds.count(dim=spatial_dims).values / nelements >= (percent/100.0)\n",
"\n",
"# Examples to check that the intensity to/from dB functions work as expected\n",
"# xx = data.vv.isel(time=0,latitude=np.arange(0, 5),longitude=np.arange(0, 5))\n",
"# xx[0] = 0 # manually change some values\n",
"# xx[1] = -0.001 # manually change some values\n",
"# display(\"digital numbers:\", xx.values)\n",
"# yy = dn_to_decibel(xx, K=-83)\n",
"# display(\"decibels:\", yy.values)\n",
"# zz = decibel_to_dn(yy, K=-83)\n",
"# display(\"digital numbers:\", zz.values)\n",
"# aa = dn_to_amplitude(xx)\n",
"# display(\"amplitude:\", aa.values)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57bcf24c-99bb-4986-9afd-3a6e63b9a8cc",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# hvPlot convenience functions\n",
"def make_image(ds: 'xarray', frame_height=300, **kwargs):\n",
" \"\"\"Return a Holoviews DynamicMap (image) object that can be displayed or combined\"\"\"\n",
" spatial_dims = ds.odc.spatial_dims\n",
" defaults = dict(\n",
" cmap=\"Greys_r\",\n",
" y = spatial_dims[0], x = spatial_dims[1],\n",
" groupby = 'time',\n",
" rasterize = True,\n",
" geo = True,\n",
" robust = True,\n",
" frame_height = frame_height,\n",
" clabel = ds.attrs.get('units', None),\n",
" )\n",
" defaults.update(**kwargs)\n",
" return ds.hvplot.image(**defaults)\n",
"\n",
"def rgb_image(ds: 'xarray', frame_height=300, **kwargs):\n",
" \"\"\"Return a Holoviews DynamicMap (RBG image) object that can be displayed or combined\"\"\"\n",
" spatial_dims = ds.odc.spatial_dims\n",
" defaults = dict(\n",
" bands='band',\n",
" y = spatial_dims[0], x = spatial_dims[1],\n",
" groupby = 'time',\n",
" rasterize = True,\n",
" geo = True,\n",
" robust = True,\n",
" frame_height = frame_height,\n",
" )\n",
" defaults.update(**kwargs)\n",
" return ds.hvplot.rgb(**defaults)\n",
"\n",
"def mask_image(ds: 'xarray', frame_height=300, **kwargs):\n",
" \"\"\"Return a Holoviews DynamicMap (mask image) object that can be displayed or combined\"\"\"\n",
" # NovaSAR ARD mask\n",
" color_def = [\n",
" (0, '#ff0004', 'InvalidData'), # red\n",
" (1, '#eeeeee', 'ValidData'), # light grey\n",
" (5, '#ff52ff', 'Layover'), # cyan\n",
" (17, '#774c0b', 'Shadow'), # brown\n",
" (18, 'black', 'upper-limit'),\n",
" ]\n",
" cvals = [x[0] for x in color_def] # values, including upper-limit\n",
" cmap = [x[1] for x in color_def[0:-1]] # colors, excluding upper-limit\n",
" cticks = [(x[0], f\"[{x[0]}] {x[2]}\") for x in color_def[0:-1]] # labels, excluding upper-limit\n",
"\n",
" # Image options\n",
" defaults = {\n",
" 'aggregator': reductions.mode(),\n",
" 'cmap': cmap,\n",
" 'clim': (cvals[0], cvals[-1]),\n",
" 'colorbar': True,\n",
" 'frame_height': frame_height,\n",
" }\n",
" # Colorbar options for categories\n",
" extra_opts = {\n",
" 'color_levels': cvals,\n",
" 'cticks': cticks,\n",
" }\n",
" defaults.update(**kwargs)\n",
" return make_image(ds, **defaults).options(hv.opts.Image(**extra_opts)) #.hist(bins=bin_edges)"
]
},
{
"cell_type": "markdown",
"id": "d0592571-8f7e-471b-acf4-57eb93563b25",
"metadata": {},
"source": [
"## Filter scenes based on valid data\n",
"\n",
"This step is optional and configurable. Here we show another simple filter:\n",
"\n",
"- Exclude time layers with less than XX% valid data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "91c113df-101a-49d2-911e-5a83e77ad14a",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Exclude time layers with less than XX% valid data\n",
"\n",
"data['mask'] = data.mask.persist()\n",
"valid_data_mask = masking.valid_data_mask(data.mask) # mask != mask.nodata -> bool\n",
"selected = select_valid_time_layers(valid_data_mask, 20) # Exclude time layers with less than 20% valid data\n",
"data = data.sel(time=selected)\n",
"\n",
"display(xarray_object_size(data))\n",
"display(data)"
]
},
{
"cell_type": "markdown",
"id": "8d244be5-8895-4d71-b9af-4ce026e14e01",
"metadata": {},
"source": [
"## Add dB and amplitude values to the dataset"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "961295e1-675b-4dda-92ba-c44b15506de8",
"metadata": {},
"outputs": [],
"source": [
"# Convenience loop for available polarization bands\n",
"\n",
"for pp in polarizations:\n",
" data[f\"{pp}_db\"] = dn_to_decibel(data[pp], K=-83).astype('float32')\n",
" data[f\"{pp}_amp\"] = dn_to_amplitude(data[pp]).astype('float32')\n",
"\n",
"display(xarray_object_size(data))\n",
"display(data)"
]
},
{
"cell_type": "markdown",
"id": "e0f94783-ad8d-4c28-b145-091bd3a5c9ae",
"metadata": {},
"source": [
"## Persist the data into the dask workers\n",
"\n",
"This will begin the data loading and calculations in the background.\n",
"\n",
"View progress in the dask dashboard (link given above where the dask cluster is defined)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f4715d63-b9ee-410b-814d-f7b94730c516",
"metadata": {},
"outputs": [],
"source": [
"# Refine for selected variables if not plotting/analysing everything\n",
"data = data.persist()"
]
},
{
"cell_type": "markdown",
"id": "5fc56259-2e20-48f8-bab5-13ef56993bf9",
"metadata": {
"tags": []
},
"source": [
"## Plot the data\n",
"\n",
"- Stronger co-polarization (VV) indicates direct backscatter while stronger cross-polarization (VH) may indicate a complex surface or volume scattering.\n",
"- Amplitude data are linear-scaled so can tend to disciminate across a range of backscatter returns.\n",
"- Decibel data are log-scaled so can tend to discriminate high and low backscatter returns.\n",
"\n",
"> Note the different data ranges for plotting (`clim`) between `vv`, `vh`, _amplitude_ and _dB_."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2b306030-ecb7-47b4-a63d-27fdb491c4f9",
"metadata": {},
"outputs": [],
"source": [
"# VV, HH and HV (amplitude and dB) and Angle hvPlots\n",
"\n",
"clim = {\n",
" 'vv_amp': (0, 0.5),\n",
" 'hh_amp': (0, 0.5),\n",
" 'hv_amp': (0, 0.25),\n",
" 'vv_db': (-20, -5),\n",
" 'hh_db': (-20, -5),\n",
" 'hv_db': (-35, -10),\n",
"}\n",
"\n",
"bands_plot_list = []\n",
"for units in ('amplitude', 'dB'):\n",
" for pp in polarizations:\n",
" varname = f\"{pp}_{units[0:3].lower()}\"\n",
" title = f\"{pp.upper()} ({units})\"\n",
" bands_plot_list.append(\n",
" make_image(data[varname], title=title, clim=clim[varname])\n",
" )\n",
"\n",
"# Which is essentially doing this for all polarizations\n",
"# vv_plot = make_image(data.vv_amp, title='VV (amplitude)', clim=(0, 0.5))\n",
"# vv_db_plot = make_image(data.vv_db, title='VV (dB)', clim=(-20, -5))\n",
"\n",
"# # Add plots for the non-polarization bands\n",
"bands_plot_list.append(make_image(data.angle, title='Incidence angle'))\n",
"bands_plot_list.append(make_image(data.scatteringarea, title='Scattering area'))\n",
"bands_plot_list.append(make_image(data.gammatosigmaratio, title='Gamma to sigma ratio'))\n",
"bands_plot_list.append(mask_image(data.mask, title='Mask'))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6531e1f9-d8c0-45bd-aeaa-f7e1c95bbc7d",
"metadata": {},
"outputs": [],
"source": [
"# Arrange plots with linked axes and time slider. Adjust browser window width if required.\n",
"\n",
"num_plot_cols = 2 if len(polarizations) <= 2 else len(polarizations)\n",
"layout = pn.panel(\n",
" hv.Layout(bands_plot_list).cols(num_plot_cols),\n",
" widget_location='top',\n",
")\n",
"\n",
"print(layout) # Helpful to see how the hvplot is constructed\n",
"display(layout)"
]
},
{
"cell_type": "markdown",
"id": "77fbc6c5-4e0e-40c2-8c8f-82757ae0909b",
"metadata": {},
"source": [
"## Plot histograms of the dB data\n",
"\n",
"A histogram can help separate water from land features. Here we show a histogram for the _dB_ channels for all time layers.\n",
"- If the histogram shows two clear peaks then a value between the peaks could be used as a water / land threshold\n",
"- If not then try selected time layers, a different area of interest, or other channels or combinations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a122cfe-f2bd-4b63-983c-c625e6849f5e",
"metadata": {},
"outputs": [],
"source": [
"# vals, bins, hist_plot = data.hv_db.plot.hist(bins=np.arange(-30, 0, 1), color='red') # Matplotlib\n",
"\n",
"hist_plot_list = []\n",
"for pp in polarizations:\n",
" varname = f\"{pp}_db\"\n",
" title = f\"{pp.upper()} (dB), combined times\"\n",
" hist_plot_list.append(\n",
" data[varname].hvplot.hist(bins=np.arange(-30, 0, 1), color='red', title=title, height=300)\n",
" )\n",
"\n",
"layout = hv.Layout(hist_plot_list).cols(len(polarizations))\n",
"\n",
"print(layout) # Helpful to see how the hvplot is constructed\n",
"display(layout)"
]
},
{
"cell_type": "markdown",
"id": "1bbb8282-3b77-43ce-a774-88b031fb94a3",
"metadata": {},
"source": [
"## Make an RGB image\n",
"\n",
"A common strategy to create an RGB colour composite image for SAR data from two channels is to use the ratio of the channels to represent the third colour.\n",
"\n",
"Here we choose:\n",
"\n",
"- For a tri-pol (3 polarization bands) product we plot each polarization as red, green and blue.\n",
"- For a dual-pol (2 polarization bands) product we plot each polarization as red and green and the ratio as blue.\n",
"- For a single-pol (1 polarization band) product we can not make a useful RGB image\n",
"\n",
"Recall that:\n",
"- VV or HH ... direct scattering\n",
"- VH or HV ... complex scattering\n",
"- ratio ... relatively more of one than the other"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "082eeed7-dc16-4098-8bad-5dc748b6a79d",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Select bands combinations for RGB\n",
"if len(polarizations) == 3:\n",
" rgb_bands = polarizations\n",
"elif len(polarizations) == 2:\n",
" combo = \"_\".join(polarizations)\n",
" data[f\"{combo}\"] = data[polarizations[0]] / data[polarizations[1]]\n",
" rgb_bands = polarizations + [combo]\n",
"else:\n",
" display(Markdown(\"**Can not make a useful RGB image from one polarization band**. Choose a dual- or tri-pol product\"))\n",
" rgb_bands = []"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "63a142d3-f278-49c8-a635-bb219675f239",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"if rgb_bands:\n",
" # Scale RGB bands by their median so they have a similar range for visualization\n",
" spatial_dims = data.odc.spatial_dims\n",
" for bb in rgb_bands:\n",
" data[f\"{bb}_scaled\"] = (data[bb] / data[bb].median(dim=spatial_dims))\n",
" rgb_bands = [f\"{bb}_scaled\" for bb in rgb_bands]\n",
"\n",
" # odc-geo function\n",
" rgb_data = data.odc.to_rgba(bands=rgb_bands, vmin=0, vmax=2).persist()\n",
"\n",
" # As subplots\n",
" # rgb_plot = rgb_image(\n",
" # rgb_data,\n",
" # ).layout().cols(4)\n",
"\n",
" # As movie. Select \"loop\" and use \"-\" button to adjust the speed to allow for rendering. After a few cycles the images should play reasonably well.\n",
" rgb_plot = rgb_image(\n",
" rgb_data,\n",
" precompute = True,\n",
" widget_type='scrubber', widget_location='bottom',\n",
" frame_height = 500,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bd2da45d-a1db-4d75-b725-706d0fa252c9",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"if rgb_bands:\n",
" print(rgb_plot) # Helpful to see how the hvplot is constructed\n",
" display(rgb_plot)"
]
},
{
"cell_type": "markdown",
"id": "fd1df15e-31e8-4ad8-9940-e28bd6822441",
"metadata": {},
"source": [
"## Export to Geotiffs\n",
"\n",
"Recall that to write a dask dataset to a file requires the dataset to be `.compute()`ed. This may result in a large memory increase on your JupyterLab node if the area of interest is large enough, which in turn may kill the kernel. If so then skip this step, choose a smaller area or find a different way to export data."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1347d080-74d5-491a-bb45-45252fc385ea",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Make a directory to save outputs to\n",
"target = Path.home() / 'output'\n",
"if not target.exists(): target.mkdir()\n",
"\n",
"def write_band(ds, varname):\n",
" \"\"\"Write the variable name of the xarray dataset to a Geotiff file for each time layer\"\"\"\n",
" for i in range(len(ds.time)):\n",
" date = ds[varname].isel(time=i).time.dt.strftime('%Y%m%d').data\n",
" fname = f'{target}/example_novasar-1_{varname}_{date}.tif'\n",
" single = ds[varname].isel(time=i).compute()\n",
" single.odc.write_cog(\n",
" fname=fname,\n",
" overwrite=True,\n",
" )\n",
" print(f'Wrote: {fname}')\n",
"\n",
"for pp in polarizations:\n",
" write_band(data, pp)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef7c0753-3ba2-4736-acf0-3f52ca36bd17",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
notebooks/data_products/sentinel-1-gamma0.ipynb
+716
View File
@@ -0,0 +1,716 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "dae926b2-fa06-4f00-b18d-7f74e92ce676",
"metadata": {},
"source": [
"## Sentinel-1 RTC Gamma0 data <img align=\"right\" src=\"../../resources/csiro_easi_logo.png\">\n",
"\n",
"#### Index\n",
"- [Overview](#Overview)\n",
"- [Setup (imports, defaults, dask, odc)](#Setup)\n",
"- [Example query](#Example-query)\n",
"- [Product definition](#Product-definition)\n",
"- [Quality layer](#Quality-layer)\n",
"- [Create and apply a good quality pixel mask](#Create-and-apply-a-good-quality-pixel-mask)\n",
"- [Plot and browse the data](#Plot-and-browse-the-data)"
]
},
{
"cell_type": "markdown",
"id": "35b91cdb-2c5a-48c0-a0bf-7de2c2637bc5",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"This notebook demonstrates how to load and use Sentinel-1 Radiometric Terrain Corrected (RTC) Gamma0 data generated in EASI.\n",
"\n",
"These _analysis ready data_ S1 gamma-0 backscatter data are processed from Sentinel-1 GRD scenes using the [SNAP-10 Toolbox](https://step.esa.int/main/download/snap-download/) with Graph Processing Tool (GPT) xml receipes. See the [RTC Gamma0 product variants](#RTC-Gamma0-product-variants) section for further details.\n",
"\n",
"For most uses we recommend the smoothed 20 m product (`sentinel1_grd_gamma0_20m`).\n",
"We can process the 10 m products (`sentinel1_grd_gamma0_10m`, `sentinel1_grd_gamma0_10m_unsmooth`) and other variants on request."
]
},
{
"cell_type": "markdown",
"id": "5e5120c6-e943-46a9-a4d8-c5b759586e62",
"metadata": {},
"source": [
"#### Using Sentinel-1 backscatter data\n",
"\n",
"An excellent introduction and overview to using SAR data is provided in the [CEOS Laymans SAR Interpretation Guide](https://ceos.org/ard/files/Laymans_SAR_Interpretation_Guide_3.0.pdf). This guide has also been converted to a set of set of Jupyter notebooks that you can download from https://github.com/AMA-Labs/cal-notebooks/tree/main/examples/SAR.\n",
"\n",
"Synthetic Aperture Radar operates in the microwave range of the electromagnetic spectrum as an active pulse sent by the satellite and scattered by features on the Earth's surface. The return signal from the surface is measured at the satellite in terms of the signal intensity, phase and polarisation compared to the signal that was sent.\n",
"\n",
"The SAR instrument on the Sentinel-1 satellites operate in the C-band at approximately 5.6 cm wavelength. This means that it can \"see\" objects of about this size and larger, and smaller objects are relatively transparent. This makes Sentinel-1 more sensitive to tree canopies, sparse and low biomass vegetation, and surface water (smooth and wind affected).\n",
"\n",
"> The SAR signal responds to the orientation and scattering from surface features of comparable size or larger than the wavelength.\n",
"> - A bright backscatter value typically means the surface was orientated perpendicular to the signal incidence angle and most of the signal was reflected back to the satellite (direct backscatter)\n",
"> - A dark backscatter value means most of the signal was reflected away from the satellite (forward scattering) and typically responds to a smooth surface (relative to the wavelength) such as calm water or bare soil\n",
"> - Rough surfaces (relative to the wavelength) result in diffuse scattering where some of the signal is returned to the satellite.\n",
"> - Complex surfaces may result in volume scattering (scattering within a tree canopy) or double-bounce scattering (perpendicular objects such as buildings and structures)\n",
"> - The relative backscatter values of co-polarisation (VV) and cross-polarisation (VH) measurements can provide information on the scattering characteristics of the surface features.\n",
"\n",
"Using Sentinel-1 backscatter data requires interpretation of the data for different surface features, including as these features change spatially or in time. It may also be necessary to carefully consider the incidence angle of the SAR signal relative to the surface features using the _incidence_angle_ band or the satellite direction metadata (descending = north to south; ascending = south to north)."
]
},
{
"cell_type": "markdown",
"id": "69f88789-9098-4669-b783-57d8982e9762",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"#### Units and conversions\n",
"The `sentinel1_grd_gamma0_*` data are given in _Intensity_ (or backscatter _power_) units. Intensity can be converted to _decibel (dB)_ or _amplitude_, and vice-versa, with the following equations. Practical _Xarray_ examples are given below.\n",
"\n",
"Intensity to/from dB:\n",
"```\n",
" dB = 10 * log10(intensity) + K\n",
"intensity = 10^((dB-K)/10)\n",
"\n",
"where K is a calibration factor, which for Sentinel-1 is 0 dB.\n",
"```\n",
"\n",
"Intensity to/from Amplitude:\n",
"```\n",
"intensity = amplitude * amplitude\n",
"amplitude = sqrt(intensity)\n",
"```\n",
"\n",
"Additional reference: https://forum.step.esa.int/t/what-stage-of-processing-requires-the-linear-to-from-db-command"
]
},
{
"cell_type": "markdown",
"id": "dd5c0ff4-e1a2-4adf-98d1-1ca5b9f3126d",
"metadata": {},
"source": [
"## Set up"
]
},
{
"cell_type": "markdown",
"id": "4cb7c012-719e-4f13-a206-36f271c6a91c",
"metadata": {},
"source": [
"#### Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c848c17-2e7c-4b68-9682-b8994644f521",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Common imports and settings\n",
"import os, sys, re\n",
"from pathlib import Path\n",
"from IPython.display import Markdown\n",
"import pandas as pd\n",
"pd.set_option(\"display.max_rows\", None)\n",
"import xarray as xr\n",
"import numpy as np\n",
"\n",
"# Datacube\n",
"import datacube\n",
"from datacube.utils.aws import configure_s3_access\n",
"import odc.geo.xr # https://github.com/opendatacube/odc-geo\n",
"from datacube.utils import masking # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py\n",
"from dea_tools.plotting import display_map # https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools\n",
"\n",
"# Basic plots\n",
"%matplotlib inline\n",
"# import matplotlib.pyplot as plt\n",
"# plt.rcParams['figure.figsize'] = [12, 8]\n",
"\n",
"# Holoviews\n",
"# https://holoviz.org/tutorial/Composing_Plots.html\n",
"# https://holoviews.org/user_guide/Composing_Elements.html\n",
"import hvplot.xarray\n",
"import panel as pn"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bde8313b-887d-4529-bd57-89f740ff6380",
"metadata": {},
"outputs": [],
"source": [
"# EASI defaults\n",
"# These are convenience functions so that the notebooks in this repository work in all EASI deployments\n",
"\n",
"# The `git.Repo()` part returns the local directory that easi-notebooks has been cloned into\n",
"# If using the `easi-tools` functions from another path, replace `repo` with your local path to `easi-notebooks` directory\n",
"try:\n",
" import git\n",
" repo = git.Repo('.', search_parent_directories=True).working_tree_dir # Path to this cloned local directory\n",
"except (ImportError, git.InvalidGitRepositoryError):\n",
" repo = Path.home() / 'easi-notebooks' # Reasonable default\n",
" if not repo.is_dir():\n",
" raise RuntimeError('To use `easi-tools` please provide the local path to `https://github.com/csiro-easi/easi-notebooks`')\n",
"if repo not in sys.path:\n",
" sys.path.append(str(repo)) # Add the local path to `easi-notebooks` to python\n",
"\n",
"from easi_tools import EasiDefaults\n",
"from easi_tools import initialize_dask, xarray_object_size, heading"
]
},
{
"cell_type": "markdown",
"id": "d2916064-8680-44a2-b537-0ecdd0b3e05b",
"metadata": {
"tags": []
},
"source": [
"#### EASI defaults\n",
"\n",
"These default values are configured for each EASI instance. They help us to use the same training notebooks in each EASI instance. You may find some of the functions convenient for your work or you can easily override the values in your copy of this notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "44d6583e-d2f5-4587-b035-6814d1313740",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"easi = EasiDefaults()\n",
"\n",
"family = 'sentinel-1'\n",
"product = easi.product(family) # 'sentinel1_grd_gamma0_20m'\n",
"display(Markdown(f'Default {family} product for \"{easi.name}\": [{product}]({easi.explorer}/products/{product})'))"
]
},
{
"cell_type": "markdown",
"id": "bb66366a-f2d5-461a-a34f-7eed5a5296fb",
"metadata": {
"tags": []
},
"source": [
"#### Dask cluster\n",
"\n",
"Using a local _Dask_ cluster is a good habit to get into. It can simplify loading and processing of data in many cases, and it provides a dashboard that shows the loading/processing progress.\n",
"\n",
"To learn more about _Dask_ see the set of [dask notebooks](https://github.com/csiro-easi/easi-notebooks/tree/main/html#dask-tutorials)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "02f9e778-755c-4ff2-a99f-ad5c48438665",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Local cluster\n",
"cluster, client = initialize_dask(workers=4)\n",
"display(client)\n",
"\n",
"# Or use Dask Gateway - this may take a few minutes\n",
"# cluster, client = initialize_dask(use_gateway=True, workers=4)\n",
"# display(client)"
]
},
{
"cell_type": "markdown",
"id": "a746582e-be79-448c-a9c4-7a2ae2d703ce",
"metadata": {},
"source": [
"#### ODC database\n",
"\n",
"Connect to the ODC database. Configure the environment and low-level tools to read from AWS buckets."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0e2232ee-f0a5-43fc-a6c6-d422b437cd69",
"metadata": {},
"outputs": [],
"source": [
"dc = datacube.Datacube()\n",
"\n",
"# Access AWS \"requester-pays\" buckets\n",
"# This is necessary for reading data from most third-party AWS S3 buckets such as for Landsat and Sentinel-2\n",
"configure_s3_access(aws_unsigned=False, requester_pays=True, client=client);"
]
},
{
"cell_type": "markdown",
"id": "1801db2b-e24e-45fa-85d4-e153ebcc8f6f",
"metadata": {},
"source": [
"## Example query\n",
"\n",
"Change any of the parameters in the `query` object below to adjust the location, time, projection, or spatial resolution of the returned datasets.\n",
"\n",
"Use the Explorer interface to check the temporal and spatial coverage for each product."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a04d875-20d4-4d72-a16d-e76b7f48161c",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Explorer link\n",
"display(Markdown(f'See: {easi.explorer}/products/{product}'))\n",
"\n",
"# EASI defaults\n",
"display(Markdown(f'#### Location: {easi.location}'))\n",
"latitude_range = easi.latitude\n",
"longitude_range = easi.longitude\n",
"time_range = easi.time\n",
"\n",
"# Or set your own latitude / longitude\n",
"# Australia GWW\n",
"# latitude_range = (-33, -32.6)\n",
"# longitude_range = (120.5, 121)\n",
"# time_range = ('2020-01-01', '2020-01-31')\n",
"\n",
"query = {\n",
" 'product': product, # Product name\n",
" 'x': longitude_range, # \"x\" axis bounds\n",
" 'y': latitude_range, # \"y\" axis bounds\n",
" 'time': time_range, # Any parsable date strings\n",
"}\n",
"\n",
"# Convenience function to display the selected area of interest\n",
"display_map(longitude_range, latitude_range)"
]
},
{
"cell_type": "markdown",
"id": "9fefeef8-0c21-49c2-b86f-c68e5178a7ac",
"metadata": {
"tags": []
},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2f4c724f-636a-4833-a165-adb053acaa07",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Target xarray parameters\n",
"# - Select a set of measurements to load\n",
"# - output CRS and resolution\n",
"# - Usually we group input scenes on the same day to a single time layer (groupby)\n",
"# - Select a reasonable Dask chunk size (this should be adjusted depending on the\n",
"# spatial and resolution parameters you choose\n",
"load_params = {\n",
" 'group_by': 'solar_day', # Scene grouping\n",
" 'dask_chunks': {'latitude':2048, 'longitude':2048}, # Dask chunks\n",
"}\n",
"\n",
"# Load data\n",
"data = dc.load(**(query | load_params))\n",
"display(xarray_object_size(data))\n",
"display(data)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "be474a1a-b5f1-40e0-b35a-5b9c9f1a15bd",
"metadata": {},
"outputs": [],
"source": [
"# When happy with the shape and size of chunks, persist() the result\n",
"data = data.persist()"
]
},
{
"cell_type": "markdown",
"id": "e3d86506-39ba-4f02-ac6e-2d1d27119668",
"metadata": {},
"source": [
"## Conversion and helper functions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "99eb8a2c-fe00-44b7-bd39-250ab3d323cc",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# These functions use numpy, which should be satisfactory for most notebooks.\n",
"# Calculations for larger or more complex arrays may require Xarray's \"ufunc\" capability.\n",
"# https://docs.xarray.dev/en/stable/examples/apply_ufunc_vectorize_1d.html\n",
"#\n",
"# Apply numpy.log10 to the DataArray\n",
"# log10_data = xr.apply_ufunc(np.log10, data)\n",
"\n",
"def intensity_to_db(da: 'xr.DataArray', K=0):\n",
" \"\"\"Return an array converted to dB values\"\"\"\n",
" xx = da.where(da > 0, np.nan) # Set values <= 0 to NaN\n",
" xx = 10*np.log10(xx) + K\n",
" xx.attrs.update({\"units\": \"dB\"})\n",
" return xx\n",
"\n",
"def db_to_intensity(da: 'xr.DataArray', K=0):\n",
" \"\"\"Return an array converted to intensity values\"\"\"\n",
" xx = np.power(10, (da-K)/10.0)\n",
" xx.attrs.update({\"units\": \"intensity\"})\n",
" return xx\n",
"\n",
"def select_valid_time_layers(ds: 'xarray', percent: float = 5):\n",
" \"\"\"Select time layers that have at least a given percentage of valid data (e.g., >=5%)\n",
"\n",
" Example usage:\n",
" selected = select_valid_time_layers(ds, percent=5)\n",
" filtered == ds.sel(time=selected)\n",
" \"\"\"\n",
" spatial_dims = ds.odc.spatial_dims\n",
" return ds.count(dim=spatial_dims).values / (ds.sizes[spatial_dims[0]]*ds.sizes[spatial_dims[1]]) >= (percent/100.0)\n",
"\n",
"# Examples to check that the intensity to/from dB functions work as expected\n",
"# xx = data.vv.isel(time=0,latitude=np.arange(0, 5),longitude=np.arange(0, 5))\n",
"# xx[0] = 0\n",
"# xx[1] = -0.001\n",
"# display(xx.values)\n",
"# yy = intensity_to_db(xx)\n",
"# display(yy.values)\n",
"# zz = db_to_intensity(yy)\n",
"# display(zz.values)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57bcf24c-99bb-4986-9afd-3a6e63b9a8cc",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# hvPlot convenience functions\n",
"def make_image(ds: 'xarray', frame_height=300, **kwargs):\n",
" \"\"\"Return a Holoviews DynamicMap (image) object that can be displayed or combined\"\"\"\n",
" spatial_dims = ds.odc.spatial_dims\n",
" defaults = dict(\n",
" cmap=\"Greys_r\",\n",
" y = spatial_dims[0], x = spatial_dims[1],\n",
" groupby = 'time',\n",
" rasterize = True,\n",
" geo = True,\n",
" robust = True,\n",
" frame_height = frame_height,\n",
" clabel = ds.attrs.get('units', None),\n",
" )\n",
" defaults.update(**kwargs)\n",
" return ds.hvplot.image(**defaults)\n",
"\n",
"def rgb_image(ds: 'xarray', frame_height=300, **kwargs):\n",
" \"\"\"Return a Holoviews DynamicMap (RBG image) object that can be displayed or combined\"\"\"\n",
" spatial_dims = ds.odc.spatial_dims\n",
" defaults = dict(\n",
" bands='band',\n",
" y = spatial_dims[0], x = spatial_dims[1],\n",
" groupby = 'time',\n",
" rasterize = True,\n",
" geo = True,\n",
" robust = True,\n",
" frame_height = frame_height,\n",
" )\n",
" defaults.update(**kwargs)\n",
" return ds.hvplot.rgb(**defaults)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "91c113df-101a-49d2-911e-5a83e77ad14a",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Optional time layer filter\n",
"\n",
"selected = select_valid_time_layers(data.vv, 10) # Exclude time layers with less than 10% valid data\n",
"data = data.sel(time=selected).persist()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "961295e1-675b-4dda-92ba-c44b15506de8",
"metadata": {},
"outputs": [],
"source": [
"# Add db values to the dataset\n",
"\n",
"data['vh_db'] = intensity_to_db(data.vh).persist()\n",
"data['vv_db'] = intensity_to_db(data.vv).persist()"
]
},
{
"cell_type": "markdown",
"id": "5fc56259-2e20-48f8-bab5-13ef56993bf9",
"metadata": {
"tags": []
},
"source": [
"## Plot the data\n",
"\n",
"> Note the different data ranges for plotting (`clim`) between `vv`, `vh`, _intensity_ and _dB_.\n",
"\n",
"- Stronger co-polarisation (VV) indicates direct backscatter while stronger cross-polarisation (VH) may indicate a complex surface or volume scattering.\n",
"- Intensity data are linear-scaled so can tend to disciminate across a range of backscatter returns.\n",
"- Decibel data are log-scaled so can tend to discriminate high and low backscatter returns."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a7630dec-2f78-41b5-b1b4-0eb32206c120",
"metadata": {},
"outputs": [],
"source": [
"# VV and VH (intensity and dB) and Angle hvPlots\n",
"\n",
"vv_plot = make_image(data.vv, clim=(0, 0.5), title='VV (intensity)')\n",
"vh_plot = make_image(data.vh, clim=(0, 0.1), title='VH (intensity)')\n",
"ia_plot = make_image(data.angle, title='Incidence angle')\n",
"\n",
"vv_db_plot = make_image(data.vv_db, clim=(-30, -3), title='VV (dB)')\n",
"vh_db_plot = make_image(data.vh_db, clim=(-30, -1), title='VH (dB)')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6531e1f9-d8c0-45bd-aeaa-f7e1c95bbc7d",
"metadata": {},
"outputs": [],
"source": [
"# Arrange plots with linked axes and time slider. Adjust browser window width if required.\n",
"\n",
"layout = pn.panel(\n",
" (vv_plot + vh_plot + ia_plot + vv_db_plot + vh_db_plot).cols(3),\n",
" widget_location='top',\n",
")\n",
"print(layout) # Helpful to see how the hvplot is constructed\n",
"layout"
]
},
{
"cell_type": "markdown",
"id": "77fbc6c5-4e0e-40c2-8c8f-82757ae0909b",
"metadata": {},
"source": [
"## Plot a histogram of the dB data\n",
"\n",
"A histogram can help separate water from land features. Here we show a histogram for the _VH (db)_ channel for all time layers.\n",
"- If the histogram shows two clear peaks then a value between the peaks could be used as a water / land threshold\n",
"- If not then try selected time layers, a different area of interest, or other channels or combinations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a122cfe-f2bd-4b63-983c-c625e6849f5e",
"metadata": {},
"outputs": [],
"source": [
"# vals, bins, hist_plot = data.vh_db.plot.hist(bins=np.arange(-30, 0, 1), color='red') # Matplotlib\n",
"hist_plot = data.vh_db.hvplot.hist(bins=np.arange(-30, 0, 1), color='red', title='Combined times', height=400) # hvPlot\n",
"\n",
"print(hist_plot) # Helpful to see how the hvplot is constructed\n",
"hist_plot"
]
},
{
"cell_type": "markdown",
"id": "1bbb8282-3b77-43ce-a774-88b031fb94a3",
"metadata": {},
"source": [
"## Make an RGB image\n",
"\n",
"A common strategy to create an RGB colour composite image for SAR data from two channels is to use the ratio of the channels to represent the third colour. Here we choose\n",
"\n",
"To create an RGB colour composite image we can use the ratio of VH and VV to represent a third channel. Here we choose\n",
"- Red = VH ... complex scattering\n",
"- Green = VV ... direct scattering\n",
"- Blue = VH/VV ... relatively more complex than direct"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bd2da45d-a1db-4d75-b725-706d0fa252c9",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Add the vh/vv band to represent 'blue'\n",
"data['vh_vv'] = data.vh / data.vv\n",
"\n",
"# Scale the measurements by their median so they have a similar range for visualization\n",
"spatial_dims = data.odc.spatial_dims\n",
"data['vh_scaled'] = data.vh / data.vh.median(dim=spatial_dims).persist()\n",
"data['vv_scaled'] = data.vv / data.vv.median(dim=spatial_dims).persist()\n",
"data['vh_vv_scaled'] = data.vh_vv / data.vh_vv.median(dim=spatial_dims).persist()\n",
"\n",
"# odc-geo function\n",
"rgb_data = data.odc.to_rgba(bands=['vh_scaled','vv_scaled','vh_vv_scaled'], vmin=0, vmax=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8834880-fb2d-42ec-bc01-f1b37ac97038",
"metadata": {},
"outputs": [],
"source": [
"# As subplots\n",
"# rgb_plot = rgb_image(\n",
"# rgb_data,\n",
"# ).layout().cols(4)\n",
"\n",
"# As movie. Select \"loop\" and use \"-\" button to adjust the speed to allow for rendering. After a few cycles the images should play reasonably well.\n",
"rgb_plot = rgb_image(\n",
" rgb_data,\n",
" precompute = True,\n",
" widget_type='scrubber', widget_location='bottom',\n",
" frame_height = 500,\n",
")\n",
"\n",
"print(rgb_plot) # Helpful to see how the hvplot is constructed\n",
"rgb_plot"
]
},
{
"cell_type": "markdown",
"id": "fd1df15e-31e8-4ad8-9940-e28bd6822441",
"metadata": {},
"source": [
"## Export to Geotiffs\n",
"\n",
"Recall that to write a dask dataset to a file requires the dataset to be `.compute()`ed. This may result in a large memory increase on your JupyterLab node if the area of interest is large enough, which in turn may kill the kernel. If so then skip this step, choose a smaller area or find a different way to export data."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1347d080-74d5-491a-bb45-45252fc385ea",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Make a directory to save outputs to\n",
"target = Path.home() / 'output'\n",
"if not target.exists(): target.mkdir()\n",
"\n",
"def write_band(ds, varname):\n",
" \"\"\"Write the variable name of the xarray dataset to a Geotiff file for each time layer\"\"\"\n",
" for i in range(len(ds.time)):\n",
" date = ds[varname].isel(time=i).time.dt.strftime('%Y%m%d').data\n",
" fname = f'{target}/example_sentinel-1_{varname}_{date}.tif'\n",
" single = ds[varname].isel(time=i).compute()\n",
" single.odc.write_cog(\n",
" fname=fname,\n",
" overwrite=True,\n",
" )\n",
" print(f'Wrote: {fname}')\n",
" \n",
"write_band(data, 'vv')\n",
"write_band(data, 'vh')"
]
},
{
"cell_type": "markdown",
"id": "7ed2add5-12ad-4327-801b-9e60b3baf031",
"metadata": {},
"source": [
"## Appendix"
]
},
{
"cell_type": "markdown",
"id": "776027ff-c4ef-4d34-962f-b2a33ba09a60",
"metadata": {},
"source": [
"### RTC Gamma0 product variants\n",
"\n",
"The set of products listed here differ by the selection and configuration of processing steps and options. The set of SNAP operators conform with [CEOS Analysis Ready Data](https://ceos.org/ard/) specifications for _normalised radar backscatter_.\n",
"\n",
"S1 gamma-0 backscatter data are processed from Sentinel-1 GRD scenes using the [SNAP-10 Toolbox](https://step.esa.int/main/download/snap-download/) with Graph Processing Tool (GPT) xml receipes (available on request).\n",
"\n",
"| | sentinel1_grd_gamma0_20m | sentinel1_grd_gamma0_10m | sentinel1_grd_gamma0_10m_unsmooth |\n",
"|--|--|--|--|\n",
"| **DEM** | | | |\n",
"| copernicus_dem_30 | Y | Y | Y |\n",
"| Scene to DEM extent multiplier| 3.0 | 3.0 | 3.0 |\n",
"| **SNAP operator** | | | |\n",
"| Apply-Orbit-File | Y | Y | Y |\n",
"| ThermalNoiseRemoval | Y | Y | Y |\n",
"| Remove-GRD-Border-Noise | Y | Y | Y |\n",
"| Calibration | Y | Y | Y |\n",
"| SetNoDataValue | Y | Y | Y |\n",
"| Terrain-Flattening | Y | Y | Y |\n",
"| Speckle-Filter | Y | Y | N |\n",
"| Multilook | Y | Y | N |\n",
"| Terrain-Correction | Y | Y | Y |\n",
"| **Output** | | | |\n",
"| Projection | WGS84, epsg:4326 | WGS84, epsg:4326 | WGS84, epsg:4326 |\n",
"| Pixel resolution | 20 m | 10 m | 10 m |\n",
"| Pixel alignment</br>_PixelIsArea = top-left_ | PixelIsArea | PixelIsArea | PixelIsArea |"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef7c0753-3ba2-4736-acf0-3f52ca36bd17",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "folium",
"language": "python",
"name": "folium"
},
"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.12.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
notebooks/data_products/sentinel-2-c1-l2a.ipynb
+632
View File
@@ -0,0 +1,632 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"# Sentinel-2 Collection 1, L2A <img align=\"right\" src=\"../../resources/csiro_easi_logo.png\">\n",
"\n",
"#### Index\n",
"- [Overview](#Overview)\n",
"- [Setup (imports, defaults, dask, odc)](#Setup)\n",
"- [Example query](#Example-query)\n",
"- [Product definition](#Product-definition)\n",
"- [Quality layer](#Quality-layer)\n",
"- [Create and apply a good quality pixel mask](#Create-and-apply-a-good-quality-pixel-mask)\n",
"- [Plot and browse the data](#Plot-and-browse-the-data)"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Overview\n",
"\n",
"Sentinel-2 is an Earth observation mission from the EU Copernicus Programme that systematically acquires optical imagery at high spatial resolution (up to 10 m for some bands). The mission is a constellation of two identical satellites in the same orbit, 180° apart for optimal coverage and data delivery. Together, they cover all Earth's land surfaces, large islands, inland and coastal waters every 3-5 days.\n",
"\n",
"Sentinel-2A was launched on 23 June 2015 and Sentinel-2B followed on 7 March 2017.\n",
"Both of the Sentinel-2 satellites carry a wide swath high-resolution multispectral imager with 13 spectral bands.\n",
"For more information see:\n",
"- [ESA Sentinel missions](https://www.esa.int/Applications/Observing_the_Earth/Copernicus/The_Sentinel_missions)\n",
"- [Sentinel-2 technical specifications](https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi)\n",
"\n",
"_Selected text adapted from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/DEA_datasets/Sentinel_2.ipynb_\n",
"\n",
"#### Data source and documentation\n",
"\n",
"ESA produces a surface reflectance \"S2 Collection 1 L2A\" product using their [sen2cor](https://step.esa.int/main/snap-supported-plugins/sen2cor/) software. [Element84](https://github.com/Element84/earth-search) convert these data to Cloud-Optimized Geotiff format and makes them publicly available at their [Earth Search STAC API](https://earth-search.aws.element84.com/v1) endpoint for programmatic access.\n",
"\n",
"EASI uses its STAC indexing tools to index datasets into our ODC databases.\n",
"\n",
"| Name | Product | Source | Information | Index\n",
"|--|--|--|--|--|\n",
"| Sentinel-2 C1 L2A COGs | `sentinel_2_c1_l2a` | [Earth Search STAC](https://earth-search.aws.element84.com/v1/collections/sentinel-2-c1-l2a) | Use for global (land) surface reflectance | Select COGS via STAC and convert to ODC metadata\n",
"\n",
"#### Collection 1 baseline processing\n",
"\n",
"ESA is reprocessing Sentinel-2 to a [\"Collection 1\"](https://sentinels.copernicus.eu/web/sentinel/sentinel-data-access/sentinel-products/sentinel-2-data-products/collection-1-level-2a) product using processing baseline >=5.00.\n",
"\n",
"Element-84 are kindly processing these data to COGs. Their processing documentation and advice is expected to be made available once the full repreocessing is completed by ESA and Element84 (see [issues](https://github.com/Element84/earth-search/issues)).\n",
"\n",
"In the meantime this notebook will show how we think the data should be loaded and used."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Common imports and settings\n",
"import os, sys\n",
"from pathlib import Path\n",
"# os.environ['USE_PYGEOS'] = '0'\n",
"from IPython.display import Markdown\n",
"import pandas as pd\n",
"pd.set_option(\"display.max_rows\", None)\n",
"import xarray as xr\n",
"\n",
"# Datacube\n",
"import datacube\n",
"from datacube.utils.aws import configure_s3_access\n",
"import odc.geo.xr # https://github.com/opendatacube/odc-geo\n",
"from datacube.utils import masking # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py\n",
"from odc.algo import enum_to_bool # https://github.com/opendatacube/odc-tools/blob/develop/libs/algo/odc/algo/_masking.py\n",
"from dea_tools.plotting import display_map, rgb # https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools\n",
"\n",
"# Basic plots\n",
"%matplotlib inline\n",
"# import matplotlib.pyplot as plt\n",
"# plt.rcParams['figure.figsize'] = [12, 8]\n",
"\n",
"# Holoviews\n",
"# https://holoviz.org/tutorial/Composing_Plots.html\n",
"# https://holoviews.org/user_guide/Composing_Elements.html\n",
"import hvplot.xarray\n",
"import panel as pn\n",
"import colorcet as cc\n",
"import cartopy.crs as ccrs\n",
"from datashader import reductions\n",
"from holoviews import opts"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# EASI defaults\n",
"# These are convenience functions so that the notebooks in this repository work in all EASI deployments\n",
"\n",
"# The `git.Repo()` part returns the local directory that easi-notebooks has been cloned into\n",
"# If using the `easi-tools` functions from another path, replace `repo` with your local path to `easi-notebooks` directory\n",
"try:\n",
" import git\n",
" repo = git.Repo('.', search_parent_directories=True).working_tree_dir # Path to this cloned local directory\n",
"except (ImportError, git.InvalidGitRepositoryError):\n",
" repo = Path.home() / 'easi-notebooks' # Reasonable default\n",
" if not repo.is_dir():\n",
" raise RuntimeError('To use `easi-tools` please provide the local path to `https://github.com/csiro-easi/easi-notebooks`')\n",
"if repo not in sys.path:\n",
" sys.path.append(str(repo)) # Add the local path to `easi-notebooks` to python\n",
"\n",
"from easi_tools import EasiDefaults\n",
"from easi_tools import initialize_dask, xarray_object_size, mostcommon_crs, heading"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"#### EASI defaults\n",
"\n",
"These default values are configured for each EASI instance. They help us to use the same training notebooks in each EASI instance. You may find some of the functions convenient for your work or you can easily override the values in your copy of this notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"easi = EasiDefaults()\n",
"\n",
"family = 'sentinel-2'\n",
"product = 'sentinel_2_c1_l2a' # Sentinel-2 collection 1, L2A\n",
"display(Markdown(f'Default {family} product for \"{easi.name}\": [{product}]({easi.explorer}/products/{product})'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Dask cluster\n",
"\n",
"Using a local _Dask_ cluster is a good habit to get into. It can simplify loading and processing of data in many cases, and it provides a dashboard that shows the loading/processing progress.\n",
"\n",
"To learn more about _Dask_ see the set of [dask notebooks](https://github.com/csiro-easi/easi-notebooks/tree/main/html#dask-tutorials)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Local cluster\n",
"cluster, client = initialize_dask(workers=4)\n",
"display(client)\n",
"\n",
"# Or use Dask Gateway - this may take a few minutes\n",
"# cluster, client = initialize_dask(use_gateway=True, workers=4)\n",
"# display(client)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### ODC database\n",
"\n",
"Connect to the ODC database. Configure the environment and low-level tools to read from AWS buckets."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dc = datacube.Datacube()\n",
"\n",
"# Access AWS \"requester-pays\" buckets\n",
"# This is necessary for reading data from most third-party AWS S3 buckets such as for Landsat and Sentinel-2\n",
"configure_s3_access(aws_unsigned=False, requester_pays=True, client=client);"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Example query\n",
"\n",
"Change any of the parameters in the `query` object below to adjust the location, time, projection, or spatial resolution of the returned datasets.\n",
"\n",
"Use the Explorer interface to check the temporal and spatial coverage for each product."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Explorer link\n",
"display(Markdown(f'See: {easi.explorer}/products/{product}'))\n",
"\n",
"# EASI defaults\n",
"display(Markdown(f'#### Location: {easi.location}'))\n",
"latitude_range = easi.latitude\n",
"longitude_range = easi.longitude\n",
"time_range = easi.time\n",
"\n",
"# Or set your own latitude / longitude\n",
"# latitude_range = (21.5, 23.5)\n",
"# longitude_range = (88, 90.8)\n",
"# time_range = ('2022-01-01', '2022-03-01')\n",
"\n",
"query = {\n",
" 'product': product, # Product name\n",
" 'x': longitude_range, # \"x\" axis bounds\n",
" 'y': latitude_range, # \"y\" axis bounds\n",
" 'time': time_range, # Any parsable date strings\n",
"}\n",
"\n",
"# Convenience function to display the selected area of interest\n",
"display_map(longitude_range, latitude_range)"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"#### Most common CRS\n",
"\n",
"Sentinel-2 datasets are stored with different coordinate reference systems (CRS), corresponding to the multiple UTM zones that are used for S2 L1B tiling. S2 measurement bands also have different resolutions (10 m, 20 m and 60 m). As such S2 queries need to include the following two query parameters:\n",
"\n",
"* `output_crs` - This sets a consistent CRS that all Sentinel-2 data will be reprojected to, irrespective of the UTM zone the individual image is stored in.\n",
"* `resolution` - This sets the resolution that all Sentinel-2 images will be resampled to. \n",
"\n",
"Use `mostcommon_crs()` to select a CRS. Adapted from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/Tools/dea_tools/datahandling.py"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Most common CRS\n",
"native_crs = mostcommon_crs(dc, query)\n",
"print(f'Most common native CRS: {native_crs}')\n",
"\n",
"# Target xarray parameters\n",
"# - Select a set of measurements to load\n",
"# - output CRS and resolution\n",
"# - Usually we group input scenes on the same day to a single time layer (groupby)\n",
"# - Select a reasonable Dask chunk size (this should be adjusted depending on the\n",
"# spatial and resolution parameters you choose\n",
"load_params = {\n",
" 'measurements': ['blue', 'red', 'green', 'nir', 'scl'], # Selected measurement or alias names\n",
" 'output_crs': native_crs, # Target EPSG code\n",
" 'resolution': (-20, 20), # Target resolution\n",
" 'group_by': 'solar_day', # Scene grouping\n",
" 'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Load data\n",
"\n",
"data = dc.load(**(query | load_params))\n",
"display(xarray_object_size(data))\n",
"display(data)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# When happy with the shape and size of chunks, persist() the result\n",
"data = data.persist()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Optional\n",
"\n",
"# Create a simple plot to verify that the data look reasonable\n",
"# This will load and create images from the data, which may take a few minutes\n",
"# Here we limit this plot to the first few time layers.\n",
"\n",
"data.isel(time=slice(0,4)).red.plot.imshow(col=\"time\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Product definition\n",
"\n",
"The product definition contains details on the measurements and quality layers available in the product. Datacube provides convenience functions that return this information in `pandas DataFrames`.\n",
"\n",
"Use `list_measurements` to show the details for a product, and `masking.describe_variable_flags` to show the flag definitions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Measurement definitions for the selected product\n",
"measurement_info = dc.list_measurements().loc[query['product']]\n",
"heading(f'Measurement table for product: {query[\"product\"]}')\n",
"display(measurement_info)\n",
"\n",
"# Flag definitions\n",
"flag_name = 'scl'\n",
"heading(f'Flag definition table for flag name: {flag_name}')\n",
"display(masking.describe_variable_flags(data[flag_name]))\n",
"\n",
"flags_def = masking.describe_variable_flags(data[flag_name]).loc['qa']['values']\n",
"display(flags_def)"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"#### Apply the correct _offset_ to the source data\n",
"\n",
"ESA introduced a change to their [L1C processing](#Collection-1-baseline-processing) that encodes their L1C and L2A products with _scale_ and _offset_ value such that\n",
"`phyiscal_value = encoded_value * scale_factor + offset`. The scale and offset details per band are available in the product definition."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Apply scale and offset\n",
"\n",
"for vv in data.data_vars:\n",
" scale = measurement_info.loc[vv, 'scale_factor']\n",
" offset = measurement_info.loc[vv, 'add_offset']\n",
" if not pd.isnull(scale) and not pd.isnull(offset):\n",
" data[vv] = data[vv] * scale + offset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Quality layer\n",
"\n",
"To visualise the **SCL** layer we create a custom color map following the colors used by ESA.\n",
"\n",
"Here we use `hvplot` to create a dynamic (zoom, scroll) image with an attached histogram. This example shows how a custom color map can be used with `hvplot`, as well as the [datashader aggregator](https://datashader.org/getting_started/Pipeline.html) `reductions.mode()`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Make SCL image\n",
"# https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm\n",
"# https://www.sentinel-hub.com/faq/how-get-s2a-scene-classification-sentinel-2/\n",
"\n",
"from bokeh.models.tickers import FixedTicker\n",
"\n",
"color_def = [\n",
" (0, '#000000', 'No data'), # black\n",
" (1, '#ff0004', 'Saturated or defective'), # red\n",
" (2, '#868686', 'Topographic and casted shadow'), # gray\n",
" (3, '#774c0b', 'Cloud shadows'), # brown\n",
" (4, '#10d32d', 'Vegetation'), # green\n",
" (5, '#ffff53', 'Not vegetated'), # yellow\n",
" (6, '#0000ff', 'Water'), # blue\n",
" (7, '#818181', 'Unclassified'), # medium gray\n",
" (8, '#c0c0c0', 'Cloud medium probability'), # light gray\n",
" (9, '#f2f2f2', 'Cloud high probability'), # very light gray\n",
" (10, '#53fff9', 'Thin cirrus'), # light blue/purple\n",
" (11, '#ff52ff', 'Snow or ice'), # cyan\n",
"]\n",
"color_val = [x[0] for x in color_def]\n",
"color_hex = [x[1] for x in color_def]\n",
"color_txt = [f'{x[0]:2d}: {x[2]}' for x in color_def]\n",
"color_lim = (min(color_val), max(color_val) + 1)\n",
"bin_edges = color_val + [max(color_val) + 1]\n",
"bin_range = (color_val[0] + 0.5, color_val[-1] + 0.5) # No idea why (0.5,11.5) works and (0,11) or (0,12) do not\n",
"\n",
"# These options manipulate the color map and colorbar to show the categories for this product\n",
"options = {\n",
" 'title': f'Flag data for: {query[\"product\"]} ({flag_name})',\n",
" 'cmap': color_hex,\n",
" 'clim': color_lim,\n",
" 'color_levels': bin_edges,\n",
" 'colorbar': True,\n",
" 'width': 800,\n",
" 'height': 450,\n",
" 'aspect': 'equal',\n",
" 'tools': ['hover'],\n",
" 'colorbar_opts': {\n",
" 'major_label_overrides': dict(zip(color_val, color_txt)),\n",
" 'major_label_text_align': 'left',\n",
" 'ticker': FixedTicker(ticks=color_val),\n",
" },\n",
"}\n",
"\n",
"# Set the dataset CRS, if using hvplot's projection and coastlines options\n",
"# plot_crs = native_crs\n",
"# if plot_crs == 'epsg:4326':\n",
"# plot_crs = ccrs.PlateCarree()\n",
"\n",
"# Native data and coastline overlay:\n",
"# - Comment `crs`, `projection`, `coastline` to plot in native_crs coords\n",
"# TODO: Update the axis labels to 'longitude', 'latitude' if `coastline` is used\n",
"\n",
"quality_plot = data[flag_name].hvplot.image(\n",
" x = 'x', y = 'y', # Dataset x,y dimension names\n",
" rasterize = True, # Use Datashader\n",
" aggregator = reductions.mode(), # Datashader selects mode value, requires 'hv.Image'\n",
" precompute = True, # Datashader precomputes what it can\n",
" # crs = plot_crs, # Datset crs\n",
" # projection = ccrs.PlateCarree(), # Output projection (ccrs.PlateCarree() when coastline=True)\n",
" # coastline = '10m', # Coastline = '10m'/'50m'/'110m'\n",
").options(opts.Image(**options)).hist(bin_range = bin_range)\n",
"\n",
"# display(quality_plot)\n",
"# Optional: Change the default time slider to a dropdown list, https://stackoverflow.com/a/54912917\n",
"fig = pn.panel(quality_plot, widgets={'time': pn.widgets.Select}) # widget_location='top_left'\n",
"fig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create and apply a good quality pixel mask\n",
"\n",
"Select a set of flag values that represent \"good quality\" for your application. Here we select \"vegetation\", \"not vegetated\" and \"water\"; that is we exclude clouds and low-quality features.\n",
"\n",
"The **SCL** layer uses distinct integer values to represent each class. The datacube `enum_to_bool()` function creates a boolean mask layer corresponding to a set of category values (string names).\n",
"\n",
"Recall that *scale* and *offset* (if required) have already been applied by the `load_s2l2a_with_offset()` function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Create Mask layer\n",
"\n",
"good_pixel_flags = [flags_def[str(i)] for i in [4, 5, 6]]\n",
"\n",
"good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags)\n",
"display(good_pixel_mask) # -> DataArray. Type: bool\n",
"\n",
"# Apply good pixel mask (multiple layers)\n",
"good_data = data.where(good_pixel_mask).persist()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate NDVI"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create an NDVI layer, as useful way to visualise the data for differences in vegetation and land cover\n",
"# ndvi = (nir - red) / (nir + red)\n",
"\n",
"ndvi = (good_data.nir - good_data.red) / (good_data.nir + good_data.red)\n",
"ndvi.persist()"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Plot and browse the data\n",
"\n",
"There are numerous tools we can use to plot and interact with the data. Here we use `hvplot` again because it works well with dask and allows us to zoom and scroll quite efficiently. `Hvplot` uses [Datashader](https://datashader.org/getting_started/Pipeline.html) to process and render only the pixels that are required for the viewport.\n",
"\n",
"Various options can be changed such as the data layer, colour map and colour range."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Generate a plot\n",
"\n",
"options = {\n",
" 'title': f'{query[\"product\"]}',\n",
" 'width': 800,\n",
" 'height': 450,\n",
" 'aspect': 'equal',\n",
" 'cmap': cc.rainbow,\n",
" 'clim': (0, 1), # Limit the color range depending on the layer_name\n",
" 'colorbar': True,\n",
" 'tools': ['hover'],\n",
"}\n",
"\n",
"# Set the dataset CRS, if using hvplot's projection and coastlines options\n",
"# plot_crs = native_crs\n",
"# if plot_crs == 'epsg:4326':\n",
"# plot_crs = ccrs.PlateCarree()\n",
"\n",
"# Native data and coastline overlay:\n",
"# - Comment `crs`, `projection`, `coastline` to plot in native_crs coords\n",
"# TODO: Update the axis labels to 'longitude', 'latitude' if `coastline` is used\n",
"\n",
"layer_plot = ndvi.hvplot.image(\n",
" x = 'x', y = 'y', # Dataset x,y dimension names\n",
" rasterize = True, # Use Datashader\n",
" aggregator = reductions.mean(), # Datashader selects mean value\n",
" precompute = True, # Datashader precomputes what it can\n",
" # crs = plot_crs, # Dataset crs\n",
" # projection = ccrs.PlateCarree(), # Output projection (use ccrs.PlateCarree() when coastline=True)\n",
" # coastline='10m', # Coastline = '10m'/'50m'/'110m'\n",
").options(opts.Image(**options)).hist(bin_range = options['clim'])\n",
"\n",
"# display(layer_plot)\n",
"# Optional: Change the default time slider to a dropdown list, https://stackoverflow.com/a/54912917\n",
"fig = pn.panel(layer_plot, widgets={'time': pn.widgets.Select}) # widget_location='top_left'\n",
"display(fig)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
notebooks/data_products/sentinel-2-l2a.ipynb
+622
View File
@@ -0,0 +1,622 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"# Sentinel-2 L2A <img align=\"right\" src=\"../../resources/csiro_easi_logo.png\">\n",
"\n",
"#### Index\n",
"- [Overview](#Overview)\n",
"- [Setup (imports, defaults, dask, odc)](#Setup)\n",
"- [Example query](#Example-query)\n",
"- [Product definition](#Product-definition)\n",
"- [Quality layer](#Quality-layer)\n",
"- [Create and apply a good quality pixel mask](#Create-and-apply-a-good-quality-pixel-mask)\n",
"- [Plot and browse the data](#Plot-and-browse-the-data)"
]
},
{
"cell_type": "markdown",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"## Overview\n",
"\n",
"Sentinel-2 is an Earth observation mission from the EU Copernicus Programme that systematically acquires optical imagery at high spatial resolution (up to 10 m for some bands). The mission is a constellation of two identical satellites in the same orbit, 180° apart for optimal coverage and data delivery. Together, they cover all Earth's land surfaces, large islands, inland and coastal waters every 3-5 days.\n",
"\n",
"Sentinel-2A was launched on 23 June 2015 and Sentinel-2B followed on 7 March 2017.\n",
"Both of the Sentinel-2 satellites carry a wide swath high-resolution multispectral imager with 13 spectral bands.\n",
"For more information see:\n",
"- [ESA Sentinel missions](https://www.esa.int/Applications/Observing_the_Earth/Copernicus/The_Sentinel_missions)\n",
"- [Sentinel-2 technical specifications](https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi)\n",
"\n",
"_Selected text adapted from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/DEA_datasets/Sentinel_2.ipynb_\n",
"\n",
"#### Data source and documentation\n",
"\n",
"ESA produces a surface reflectance \"S2 L2A\" product using their [sen2cor](https://step.esa.int/main/snap-supported-plugins/sen2cor/) software. [Element84](https://github.com/Element84/earth-search) convert these data to Cloud-Optimized Geotiff format and makes them publicly available at their [Earth Search STAC API](https://earth-search.aws.element84.com/v1) endpoint, and [AWS open data](https://registry.opendata.aws/sentinel-2-l2a-cogs/), for programmatic access.\n",
"\n",
"EASI uses its STAC indexing tools to index datasets into our ODC databases.\n",
"\n",
"| Name | Product | Source | Information | Index\n",
"|--|--|--|--|--|\n",
"| Sentinel-2 COGs | `s2_l2a` | [Earth Search STAC](https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a) | Use for global (land) surface reflectance | Select COGS via STAC and convert to ODC metadata\n",
"\n",
"#### Baseline processing and offset consideration\n",
"\n",
"ESA introduced a change to their L1C processing (processing baseline >=04.00, for data since 25 January 2022) that encodes their L1C and L2A products with an _offset_ value such that\n",
"`phyiscal_value = encoded_value * scale_factor + offset`\n",
"\n",
"Element84 additionally [may have pre-applied the offset or not](https://github.com/Element84/earth-search/issues/23#issuecomment-1834674853) to the data that we index and load. This introduces an inconsistency in the S2A series that we need to account for.\n",
"\n",
"In this notebook we [show how to](#Apply-the-correct-offset-to-the-source-data) correctly load `s2_l2a` data with the offset applied to scenes where required. We use a `load_s2l2a_with_offset` function provided in `easi_tools` that identifies which scenes within a query should and should not have the offset applied, loads the data and applies the offset, and combines the non-offsetted with the offsetted scenes into a single xarray Dataset. This method also applies the scale_factor to all data in the query."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Common imports and settings\n",
"import os, sys, re\n",
"from pathlib import Path\n",
"from IPython.display import Markdown\n",
"import pandas as pd\n",
"pd.set_option(\"display.max_rows\", None)\n",
"import xarray as xr\n",
"\n",
"# Datacube\n",
"import datacube\n",
"from datacube.utils.aws import configure_s3_access\n",
"import odc.geo.xr # https://github.com/opendatacube/odc-geo\n",
"from datacube.utils import masking # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py\n",
"from odc.algo import enum_to_bool # https://github.com/opendatacube/odc-tools/blob/develop/libs/algo/odc/algo/_masking.py\n",
"from dea_tools.plotting import display_map, rgb # https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools\n",
"\n",
"# Basic plots\n",
"%matplotlib inline\n",
"# import matplotlib.pyplot as plt\n",
"# plt.rcParams['figure.figsize'] = [12, 8]\n",
"\n",
"# Holoviews\n",
"# https://holoviz.org/tutorial/Composing_Plots.html\n",
"# https://holoviews.org/user_guide/Composing_Elements.html\n",
"import hvplot.pandas\n",
"import hvplot.xarray\n",
"import panel as pn\n",
"import colorcet as cc\n",
"import cartopy.crs as ccrs\n",
"from datashader import reductions\n",
"from holoviews import opts\n",
"# hv.extension('bokeh', logo=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# EASI defaults\n",
"# These are convenience functions so that the notebooks in this repository work in all EASI deployments\n",
"\n",
"# The `git.Repo()` part returns the local directory that easi-notebooks has been cloned into\n",
"# If using the `easi-tools` functions from another path, replace `repo` with your local path to `easi-notebooks` directory\n",
"try:\n",
" import git\n",
" repo = git.Repo('.', search_parent_directories=True).working_tree_dir # Path to this cloned local directory\n",
"except (ImportError, git.InvalidGitRepositoryError):\n",
" repo = Path.home() / 'easi-notebooks' # Reasonable default\n",
" if not repo.is_dir():\n",
" raise RuntimeError('To use `easi-tools` please provide the local path to `https://github.com/csiro-easi/easi-notebooks`')\n",
"if repo not in sys.path:\n",
" sys.path.append(str(repo)) # Add the local path to `easi-notebooks` to python\n",
"\n",
"from easi_tools import EasiDefaults\n",
"from easi_tools import initialize_dask, xarray_object_size, mostcommon_crs, heading\n",
"from easi_tools.load_s2l2a import load_s2l2a_with_offset"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"#### EASI defaults\n",
"\n",
"These default values are configured for each EASI instance. They help us to use the same training notebooks in each EASI instance. You may find some of the functions convenient for your work or you can easily override the values in your copy of this notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"easi = EasiDefaults()\n",
"\n",
"family = 'sentinel-2'\n",
"product = 's2_l2a' # Sentinel-2 collection 0, L2A\n",
"display(Markdown(f'Default {family} product for \"{easi.name}\": [{product}]({easi.explorer}/products/{product})'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Dask cluster\n",
"\n",
"Using a local _Dask_ cluster is a good habit to get into. It can simplify loading and processing of data in many cases, and it provides a dashboard that shows the loading/processing progress.\n",
"\n",
"To learn more about _Dask_ see the set of [dask notebooks](https://github.com/csiro-easi/easi-notebooks/tree/main/html#dask-tutorials)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Local cluster\n",
"cluster, client = initialize_dask(workers=4)\n",
"display(client)\n",
"\n",
"# Or use Dask Gateway - this may take a few minutes\n",
"# cluster, client = initialize_dask(use_gateway=True, workers=4)\n",
"# display(client)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### ODC database\n",
"\n",
"Connect to the ODC database. Configure the environment and low-level tools to read from AWS buckets."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dc = datacube.Datacube()\n",
"\n",
"# Access AWS \"requester-pays\" buckets\n",
"# This is necessary for reading data from most third-party AWS S3 buckets such as for Landsat and Sentinel-2\n",
"configure_s3_access(aws_unsigned=False, requester_pays=True, client=client);"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Example query\n",
"\n",
"Change any of the parameters in the `query` object below to adjust the location, time, projection, or spatial resolution of the returned datasets.\n",
"\n",
"Use the Explorer interface to check the temporal and spatial coverage for each product."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Explorer link\n",
"display(Markdown(f'See: {easi.explorer}/products/{product}'))\n",
"\n",
"# EASI defaults\n",
"display(Markdown(f'#### Location: {easi.location}'))\n",
"latitude_range = easi.latitude\n",
"longitude_range = easi.longitude\n",
"time_range = easi.time\n",
"\n",
"# Or set your own latitude / longitude\n",
"# latitude_range = (-36.3, -35.8)\n",
"# longitude_range = (146.8, 147.3)\n",
"# time_range = ('2022-01-01', '2022-03-01')\n",
"\n",
"query = {\n",
" 'product': product, # Product name\n",
" 'x': longitude_range, # \"x\" axis bounds\n",
" 'y': latitude_range, # \"y\" axis bounds\n",
" 'time': time_range, # Any parsable date strings\n",
"}\n",
"\n",
"# Convenience function to display the selected area of interest\n",
"display_map(longitude_range, latitude_range)"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"#### Most common CRS\n",
"\n",
"Sentinel-2 datasets are stored with different coordinate reference systems (CRS), corresponding to the multiple UTM zones that are used for S2 L1B tiling. S2 measurement bands also have different resolutions (10 m, 20 m and 60 m). As such S2 queries need to include the following two query parameters:\n",
"\n",
"* `output_crs` - This sets a consistent CRS that all Sentinel-2 data will be reprojected to, irrespective of the UTM zone the individual image is stored in.\n",
"* `resolution` - This sets the resolution that all Sentinel-2 images will be resampled to. \n",
"\n",
"Use `mostcommon_crs()` to select a CRS. Adapted from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/Tools/dea_tools/datahandling.py"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Most common CRS\n",
"native_crs = mostcommon_crs(dc, query)\n",
"print(f'Most common native CRS: {native_crs}')\n",
"\n",
"# Target xarray parameters\n",
"# - Select a set of measurements to load\n",
"# - output CRS and resolution\n",
"# - Usually we group input scenes on the same day to a single time layer (groupby)\n",
"# - Select a reasonable Dask chunk size (this should be adjusted depending on the\n",
"# spatial and resolution parameters you choose\n",
"load_params = {\n",
" 'measurements': ['blue', 'red', 'green', 'nir', 'scl'], # Selected measurement or alias names\n",
" 'output_crs': native_crs, # Target EPSG code\n",
" 'resolution': (-20, 20), # Target resolution\n",
" 'group_by': 'solar_day', # Scene grouping\n",
" 'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"#### Apply the correct _offset_ to the source data\n",
"\n",
"We provide a convenience function to load `s2_l2a` data, apply the scale (and offset if required) and return an `xarray.Dataset`. This replaces `datacube.load()` for this product."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# The usual dc.load() may give an incorrect result for this product due to inconsistent upstream handling of any offset value\n",
"# data = dc.load(query | load_params)\n",
"\n",
"# Use the replacement \"dc.load()\" function for this product\n",
"data = load_s2l2a_with_offset(dc, query | load_params) # Combine the two dicts that contain our search and load parameters\n",
"display(xarray_object_size(data))\n",
"display(data)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# When happy with the shape and size of chunks, persist() the result\n",
"data = data.persist()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Optional\n",
"\n",
"# Create a simple plot to verify that the data look reasonable\n",
"# This will load and create images from the data, which may take a few minutes\n",
"# Here we limit this plot to the first few time layers.\n",
"\n",
"data.isel(time=slice(0, 4)).red.plot.imshow(col=\"time\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Product definition\n",
"\n",
"The product definition contains details on the measurements and quality layers available in the product. Datacube provides convenience functions that return this information in `pandas DataFrames`.\n",
"\n",
"Use `list_measurements` to show the details for a product, and `masking.describe_variable_flags` to show the flag definitions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Measurement definitions for the selected product\n",
"measurement_info = dc.list_measurements().loc[query['product']]\n",
"heading(f'Measurement table for product: {query[\"product\"]}')\n",
"display(measurement_info)\n",
"\n",
"# Flag definitions\n",
"flag_name = 'scl'\n",
"heading(f'Flag definition table for flag name: {flag_name}')\n",
"display(masking.describe_variable_flags(data[flag_name]))\n",
"\n",
"flags_def = masking.describe_variable_flags(data[flag_name]).loc['qa']['values']\n",
"display(flags_def)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Quality layer\n",
"\n",
"To visualise the **SCL** layer we create a custom color map following the colors used by ESA.\n",
"\n",
"Here we use `hvplot` to create a dynamic (zoom, scroll) image with an attached histogram. This example shows how a custom color map can be used with `hvplot`, as well as the [datashader aggregator](https://datashader.org/getting_started/Pipeline.html) `reductions.mode()`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Make SCL image\n",
"# https://sentiwiki.copernicus.eu/web/s2-processing\n",
"# https://www.sentinel-hub.com/faq/how-get-s2a-scene-classification-sentinel-2/\n",
"\n",
"from bokeh.models.tickers import FixedTicker\n",
"\n",
"color_def = [\n",
" (0, '#000000', 'No data'), # black\n",
" (1, '#ff0004', 'Saturated or defective'), # red\n",
" (2, '#868686', 'Dark features or shadows'), # gray\n",
" (3, '#774c0b', 'Cloud shadows'), # brown\n",
" (4, '#10d32d', 'Vegetation'), # green\n",
" (5, '#ffff53', 'Not vegetated'), # yellow\n",
" (6, '#0000ff', 'Water'), # blue\n",
" (7, '#818181', 'Unclassified'), # medium gray\n",
" (8, '#c0c0c0', 'Cloud medium probability'), # light gray\n",
" (9, '#f2f2f2', 'cloud high probability'), # very light gray\n",
" (10, '#53fff9', 'Thin cirrus'), # cyan\n",
" (11, '#ff52ff', 'Snow or ice'), # pink\n",
"]\n",
"color_val = [x[0] for x in color_def]\n",
"color_hex = [x[1] for x in color_def]\n",
"color_txt = [f'{x[0]:2d}: {x[2]}' for x in color_def]\n",
"color_lim = (min(color_val), max(color_val) + 1)\n",
"bin_edges = color_val + [max(color_val) + 1]\n",
"bin_range = (color_val[0] + 0.5, color_val[-1] + 0.5) # No idea why (0.5,11.5) works and (0,11) or (0,12) do not\n",
"\n",
"# These options manipulate the color map and colorbar to show the categories for this product\n",
"options = {\n",
" 'title': f'Flag data for: {query[\"product\"]} ({flag_name})',\n",
" 'cmap': color_hex,\n",
" 'clim': color_lim,\n",
" 'color_levels': bin_edges,\n",
" 'colorbar': True,\n",
" 'width': 800,\n",
" 'height': 450,\n",
" 'aspect': 'equal',\n",
" 'tools': ['hover'],\n",
" 'colorbar_opts': {\n",
" 'major_label_overrides': dict(zip(color_val, color_txt)),\n",
" 'major_label_text_align': 'left',\n",
" 'ticker': FixedTicker(ticks=color_val),\n",
" },\n",
"}\n",
"\n",
"# Set the dataset CRS, if using hvplot's projection and coastlines options\n",
"# plot_crs = native_crs\n",
"# if plot_crs == 'epsg:4326':\n",
"# plot_crs = ccrs.PlateCarree()\n",
"\n",
"# Native data and coastline overlay:\n",
"# - Comment `crs`, `projection`, `coastline` to plot in native_crs coords\n",
"# TODO: Update the axis labels to 'longitude', 'latitude' if `coastline` is used\n",
"\n",
"quality_plot = data[flag_name].hvplot.image(\n",
" x = 'x', y = 'y', # Dataset x,y dimension names\n",
" rasterize = True, # Use Datashader\n",
" aggregator = reductions.mode(), # Datashader selects mode value, requires 'hv.Image'\n",
" precompute = True, # Datashader precomputes what it can\n",
" # crs = plot_crs, # Datset crs\n",
" # projection = ccrs.PlateCarree(), # Output projection (ccrs.PlateCarree() when coastline=True)\n",
" # coastline = '10m', # Coastline = '10m'/'50m'/'110m'\n",
").options(opts.Image(**options)).hist(bin_range = bin_range)\n",
"\n",
"# display(quality_plot)\n",
"# Optional: Change the default time slider to a dropdown list, https://stackoverflow.com/a/54912917\n",
"fig = pn.panel(quality_plot, widgets={'time': pn.widgets.Select}) # widget_location='top_left'\n",
"fig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create and apply a good quality pixel mask\n",
"\n",
"Select a set of flag values that represent \"good quality\" for your application. Here we select \"vegetation\", \"not vegetated\" and \"water\"; that is we exclude clouds and low-quality features.\n",
"\n",
"The **SCL** layer uses distinct integer values to represent each class. The datacube `enum_to_bool()` function creates a boolean mask layer corresponding to a set of category values (string names).\n",
"\n",
"Recall that *scale* and *offset* (if required) have already been applied by the `load_s2l2a_with_offset()` function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Create Mask layer\n",
"\n",
"good_pixel_flags = [flags_def[str(i)] for i in [4, 5, 6]]\n",
"\n",
"good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags)\n",
"display(good_pixel_mask) # -> DataArray. Type: bool\n",
"\n",
"# Apply good pixel mask (multiple layers)\n",
"good_data = data.where(good_pixel_mask).persist()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate NDVI"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create an NDVI layer, as useful way to visualise the data for differences in vegetation and land cover\n",
"# ndvi = (nir - red) / (nir + red)\n",
"\n",
"ndvi = (good_data.nir - good_data.red) / (good_data.nir + good_data.red)\n",
"ndvi.persist()"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## Plot and browse the data\n",
"\n",
"There are numerous tools we can use to plot and interact with the data. Here we use `hvplot` again because it works well with dask and allows us to zoom and scroll quite efficiently. `Hvplot` uses [Datashader](https://datashader.org/getting_started/Pipeline.html) to process and render only the pixels that are required for the viewport.\n",
"\n",
"Various options can be changed such as the data layer, colour map and colour range."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Generate a plot\n",
"\n",
"options = {\n",
" 'title': f'{query[\"product\"]}',\n",
" 'width': 800,\n",
" 'height': 450,\n",
" 'aspect': 'equal',\n",
" 'cmap': cc.rainbow,\n",
" 'clim': (0, 1), # Limit the color range depending on the layer_name\n",
" 'colorbar': True,\n",
" 'tools': ['hover'],\n",
"}\n",
"\n",
"# Set the dataset CRS, if using hvplot's projection and coastlines options\n",
"# plot_crs = native_crs\n",
"# if plot_crs == 'epsg:4326':\n",
"# plot_crs = ccrs.PlateCarree()\n",
"\n",
"# Native data and coastline overlay:\n",
"# - Comment `crs`, `projection`, `coastline` to plot in native_crs coords\n",
"# TODO: Update the axis labels to 'longitude', 'latitude' if `coastline` is used\n",
"\n",
"layer_plot = ndvi.hvplot.image(\n",
" x = 'x', y = 'y', # Dataset x,y dimension names\n",
" rasterize = True, # Use Datashader\n",
" aggregator = reductions.mean(), # Datashader selects mean value\n",
" precompute = True, # Datashader precomputes what it can\n",
" # crs = plot_crs, # Dataset crs\n",
" # projection = ccrs.PlateCarree(), # Output projection (use ccrs.PlateCarree() when coastline=True)\n",
" # coastline='10m', # Coastline = '10m'/'50m'/'110m'\n",
").options(opts.Image(**options)).hist(bin_range = options['clim'])\n",
"\n",
"# display(layer_plot)\n",
"# Optional: Change the default time slider to a dropdown list, https://stackoverflow.com/a/54912917\n",
"fig = pn.panel(layer_plot, widgets={'time': pn.widgets.Select}) # widget_location='top_left'\n",
"display(fig)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}