diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..a8c2003 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:conda", + "python-envs.defaultPackageManager": "ms-python.python:conda", + "python-envs.pythonProjects": [] +} \ No newline at end of file diff --git a/__pycache__/new_import_ODC.cpython-310.pyc b/__pycache__/new_import_ODC.cpython-310.pyc new file mode 100644 index 0000000..6503c27 Binary files /dev/null and b/__pycache__/new_import_ODC.cpython-310.pyc differ diff --git a/__pycache__/new_import_ODC.cpython-313.pyc b/__pycache__/new_import_ODC.cpython-313.pyc new file mode 100644 index 0000000..50eb3e8 Binary files /dev/null and b/__pycache__/new_import_ODC.cpython-313.pyc differ diff --git a/__pycache__/utils.cpython-310.pyc b/__pycache__/utils.cpython-310.pyc new file mode 100644 index 0000000..5c973ed Binary files /dev/null and b/__pycache__/utils.cpython-310.pyc differ diff --git a/deafrica_tools/__pycache__/__init__.cpython-310.pyc b/deafrica_tools/__pycache__/__init__.cpython-310.pyc index 068b02d..6328f14 100644 Binary files a/deafrica_tools/__pycache__/__init__.cpython-310.pyc and b/deafrica_tools/__pycache__/__init__.cpython-310.pyc differ diff --git a/deafrica_tools/__pycache__/bandindices.cpython-310.pyc b/deafrica_tools/__pycache__/bandindices.cpython-310.pyc index e0c5b54..288e1d5 100644 Binary files a/deafrica_tools/__pycache__/bandindices.cpython-310.pyc and b/deafrica_tools/__pycache__/bandindices.cpython-310.pyc differ diff --git a/easi_tools/__init__.py b/easi_tools/__init__.py new file mode 100644 index 0000000..815de9b --- /dev/null +++ b/easi_tools/__init__.py @@ -0,0 +1,9 @@ +#!python 3 + +from .deployments import EasiDefaults +from .notebook_utils import \ + heading, \ + initialize_dask, \ + mostcommon_crs, \ + unset_cachingproxy, \ + xarray_object_size diff --git a/easi_tools/__pycache__/__init__.cpython-310.pyc b/easi_tools/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000..4289615 Binary files /dev/null and b/easi_tools/__pycache__/__init__.cpython-310.pyc differ diff --git a/easi_tools/__pycache__/deployments.cpython-310.pyc b/easi_tools/__pycache__/deployments.cpython-310.pyc new file mode 100644 index 0000000..49f57ac Binary files /dev/null and b/easi_tools/__pycache__/deployments.cpython-310.pyc differ diff --git a/easi_tools/__pycache__/load_s2l2a.cpython-310.pyc b/easi_tools/__pycache__/load_s2l2a.cpython-310.pyc new file mode 100644 index 0000000..7671a64 Binary files /dev/null and b/easi_tools/__pycache__/load_s2l2a.cpython-310.pyc differ diff --git a/easi_tools/__pycache__/notebook_utils.cpython-310.pyc b/easi_tools/__pycache__/notebook_utils.cpython-310.pyc new file mode 100644 index 0000000..6dcc74f Binary files /dev/null and b/easi_tools/__pycache__/notebook_utils.cpython-310.pyc differ diff --git a/easi_tools/deployments.py b/easi_tools/deployments.py new file mode 100644 index 0000000..f2e753f --- /dev/null +++ b/easi_tools/deployments.py @@ -0,0 +1,335 @@ +#!python3 + +import sys +import os +import logging +import collections + +# A class that provides notebook variables for each of the EASI deployments + +# Map an internal deployment name to deployment variables and search parameters. +# Update to ensure that the product/space/time parameters are available in the respective databases +deployment_map = { + 'adias': { + 'domain': 'adias.aquawatchaus.space', + 'db_database': 'adias_prod_db', + 'training_shapefile': '', + 'scratch': 'adias-prod-user-scratch', + 'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 's2_l2a', 'sar': 'asf_s1_grd_gamma0', 'dem': 'copernicus_dem_30'}, + 'location': 'Lake Tahoe, California', + 'latitude': (39.0, 39.3), + 'longitude': (-120.2, -119.9), + 'time': ('2022-02-01', '2022-05-01'), + 'target': { + 'landsat': {'crs': 'epsg:26911', 'resolution': (-30,30)}, + 'sentinel-2': {'crs': 'epsg:26911', 'resolution': (-10,10)} + }, + 'aliases': { + 'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'} + }, + 'qa_mask': { + 'landsat': {'nodata': False, 'water': 'land_or_cloud', + 'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'} + } + }, + 'asia': { + 'domain': 'asia.easi-eo.solutions', + 'db_database': 'easi_asia_db', + 'training_shapefile': '', + 'scratch': 'easi-asia-user-scratch', + 'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 'sentinel_2_c1_l2a', 'sentinel-1': 'sentinel1_grd_gamma0_20m', 'dem': 'copernicus_dem_30'}, + 'location': 'Lake Tempe, Indonesia', + 'latitude': (-4.2, -3.9), + 'longitude': (119.8, 120.1), + 'time': ('2020-02-01', '2020-04-01'), + 'proxy': True, + 'target': { + 'landsat': {'crs': 'epsg:32650', 'resolution': (-30,30)}, + 'sentinel-2': {'crs': 'epsg:32650', 'resolution': (-10,10)} + }, + 'aliases': { + 'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'} + }, + 'qa_mask': { + 'landsat': {'nodata': False, 'water': 'land_or_cloud', + 'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'} + } + }, + 'chile': { + 'domain': 'datacubechile.cl', + 'db_database': 'easido_prod_db', + 'training_shapefile': '', + 'scratch': 'easido-prod-user-scratch', + 'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 'sentinel_2_c1_l2a', 'sar': 'asf_s1_grd_gamma0', 'dem': 'copernicus_dem_30'}, + 'location': 'La Serena, Chile', + 'latitude': (-29.95, -29.85), + 'longitude': (-71.3, -71.2), + 'latitude_big': (-29.95, -27.95), + 'longitude_big': (-71.3, -69.3), + 'time': ('2022-02-01', '2022-05-01'), + 'target': { + 'landsat': {'crs': 'epsg:32718', 'resolution': (-30,30)}, + 'sentinel-2': {'crs': 'epsg:32718', 'resolution': (-10,10)} + }, + 'aliases': { + 'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'} + }, + 'qa_mask': { + 'landsat': {'nodata': False, 'water': 'land_or_cloud', + 'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'} + } + }, + 'cal': { + 'domain': 'cal.ceos.org', + 'db_database': 'ceoseail_eail_db', + 'training_shapefile': './ancillary_data/VA_Counties_Newport_News.shp', + 'scratch': 'ceoseail-eail-user-scratch', + 'ows': False, + 'map': False, + 'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 's2_l2a', 'sentinel-1': 's1_rtc', 'dem': 'copernicus_dem_30'}, + 'location': 'Newport News, Virginia', + 'latitude': (37.02, 37.12), + 'longitude': (-76.55, -76.45), + 'time': ('2022-01-01', '2022-04-01'), + 'target': { + 'landsat': {'crs': 'epsg:32618', 'resolution': (-30,30)}, + 'sentinel-2': {'crs': 'epsg:32618', 'resolution': (-10,10)} + }, + 'aliases': { + 'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'} + }, + 'qa_mask': { + 'landsat': {'nodata': False, 'water': 'land_or_cloud', + 'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'} + } + }, + 'csiro': { + 'domain': 'csiro.easi-eo.solutions', + 'db_database': 'easihub_csiro_db', + 'training_shapefile': '', + 'scratch': 'easihub-csiro-user-scratch', + 'productmap': {'landsat': 'ga_ls8c_ard_3', 'sentinel-2': 'ga_s2am_ard_3', 'sentinel-1': 'sentinel1_grd_gamma0_20m', 'dem': 'copernicus_dem_30'}, + 'location': 'Lake Hume, Australia', + 'latitude': (-36.3, -35.8), + 'longitude': (146.8, 147.3), + 'time': ('2020-02-01', '2020-04-01'), + 'aliases': { + 'landsat': {'red': 'nbart_red', 'green': 'nbart_green', 'blue': 'nbart_blue', + 'nir': 'nbart_nir', 'swir1': 'nbart_swir_1', 'swir2': 'nbart_swir_2', + 'qa_band': 'oa_fmask'} + }, + 'qa_mask': { + 'landsat': {'fmask':'valid'} + } + }, + 'sub-apse2': { + 'domain': 'sub-apse2.easi-eo.solutions', + 'db_database': '', + 'training_shapefile': '', + 'scratch': '', + 'ows': False, + 'map': False, + 'productmap': {'landsat': 'ga_ls8c_ard_3', 'sentinel-2': 'ga_s2am_ard_3', 'dem': 'copernicus_dem_30'}, + 'location': 'Lake Hume, Australia', + 'latitude': (-36.3, -35.8), + 'longitude': (146.8, 147.3), + 'time': ('2020-02-01', '2020-04-01'), + 'aliases': { + 'landsat': {'red': 'nbart_red', 'green': 'nbart_green', 'blue': 'nbart_blue', + 'nir': 'nbart_nir', 'swir1': 'nbart_swir_1', 'swir2': 'nbart_swir_2', + 'qa_band': 'oa_fmask'} + }, + 'qa_mask': { + 'landsat': {'fmask':'valid'} + } + }, +} + + +class EasiDefaults(): + """Provide deployment-specific default variables for EASI notebooks""" + + def __init__(self, deployment=None): + """Initialise""" + self._log = _getlogger(self.__class__.__name__) + self.name = deployment if deployment else self._find_deployment() + self.deployment = self._validate(self.name) + self.proxy = None + self._aliases = {} + if self.deployment and self.deployment.get('proxy', None): + self.proxy = EasiCachingProxy() + if self.deployment: + self._log.info(f'Successfully found configuration for deployment "{self.name}"') + + def _validate(self, deployment) -> dict: + """Return the dict associated with the deployment name""" + names = deployment_map.keys() + if deployment is None or deployment not in names: + self._log.error(f'Deployment name not recognised: {deployment}') + self._log.error(f'Select one of: {", ".join(names)}') + return None + return deployment_map[deployment] + + def _find_deployment(self) -> str: + """Use the deployment's database environment variable as a lookup into the deployment_map dict""" + db_database = os.environ['DB_DATABASE'] + deployment_name = [item for item in deployment_map if deployment_map[item]["db_database"] == db_database] + msg = 'Try specifying one using EasiDefaults(deployment="deployment_name").' + if len(deployment_name) == 0: + self._log.error(f'Deployment could not be found automatically. {msg}') + return None + elif len(deployment_name) > 1: + self._log.error(f'More than one deployment found. {msg}') + return None + return deployment_name[0] + + + @property + def domain(self): + """Deployment domain""" + return self.deployment['domain'] + + @property + def db_database(self): + """Database name""" + return self.deployment['db_database'] + + @property + def training_shapefile(self): + """A local shapefile""" + return self.deployment['training_shapefile'] + + @property + def hub(self): + """JupyterLab URL""" + return f'https://hub.{self.domain}' + + @property + def explorer(self): + """Explorer URL""" + return f'https://explorer.{self.domain}' + + @property + def ows(self): + """OWS URL""" + if not self.deployment.get('ows', True): + self._log.warning(f'Deployment does not have an OWS service: {self.name}') + return None + return f'https://ows.{self.domain}' + + @property + def terria(self): + """Terria Map URL""" + if not self.deployment.get('map', True): + self._log.warning(f'Deployment does not have a Map service: {self.name}') + return None + return f'https://map.{self._domain()}' + + @property + def scratch(self): + """Scratch bucket""" + return self.deployment['scratch'] + + @property + def location(self): + """Default location name""" + return self.deployment['location'] + + @property + def latitude(self): + """Default latitude range""" + return self.deployment['latitude'] + + @property + def longitude(self): + """Default longitude range""" + return self.deployment['longitude'] + + @property + def latitude_big(self): + """Default big latitude range""" + if 'latitude_big' in self.deployment: + return self.deployment['latitude_big'] + self._log.warning(f'Default big latitude range not defined for "{self.deployment}". Using default latitude range') + return self.latitude + + @property + def longitude_big(self): + """Default big longitude range""" + if 'longitude_big' in self.deployment: + return self.deployment['longitude_big'] + self._log.warning(f'Default big longitude range not defined for "{self.deployment}". Using default longitude range') + return self.latitude + + @property + def time(self): + """Default time range""" + return self.deployment['time'] + + def product(self, family='landsat'): + """Product name. Family loosely describes products from a satellite series or product type.""" + p = self.deployment['productmap'].get(family, None) + if p is None: + self._log.warning(f'Product family not defined for "{self.name}": {family}') + out = ', '.join([f'{k} > {v}' for k,v in self.deployment['productmap'].items()]) + self._log.warning(f'{self.name}: {out}') + return None + return p + + def crs(self, family='landsat'): + """Default resolution. Family loosely describes products from a satellite series or product type.""" + return self.deployment.get('target', {}).get(family, {}).get('crs', None) + + def resolution(self, family='landsat'): + """Default resolution. Family loosely describes products from a satellite series or product type.""" + return self.deployment.get('target', {}).get(family, {}).get('resolution', None) + + def aliases(self, family='landsat') -> collections.UserDict: + """Return a dict-like object that maps a common name to a specific measurement/alias name. + Family loosely describes products from a satellite series or product type. + + The common name is returned if there is no specific measurement/alias name defined. + That is, the common name should work as a measurement/alias name for the family in this deployment. + Else, provide a specific measurement/alias name in the defaults above. + """ + if family not in self._aliases: + self._aliases[family] = EasiAlias(self.deployment.get('aliases', {}).get(family, {})) + return self._aliases[family] + + def qa_mask(self, family='landsat') -> dict: + """Default QA mask values. Family loosely describes products from a satellite series or product type.""" + return self.deployment.get('qa_mask', {}).get(family, {}) + + +class EasiAlias(collections.UserDict): + """Custom UserDict that returns a default measurement name for a given key if defined. + Else returns the key as the value. Items can not be set.""" + def __init__(self, default:dict = {}): + self.data = default + self._log = _getlogger(self.__class__.__name__) + def __getitem__(self, key): + if key in self.data: + return self.data[key] + return key + def __setitem__(self, key, val): + self._log.error(f'Error <{self.__class__.__name__}>: Can not set items') + + +class EasiCachingProxy(): + """Set, unset and return information about the user's caching-proxy configuration""" + def __init__(self): + pass + + +def _getlogger(name): + """Return a logger. Define here to limit external dependencies""" + # Default logger + # log.hasHandlers() = False + # log.getEffectiveLevel() = 30 = warning + # log.propagate = True + logger = logging.getLogger(name) + logger.setLevel(logging.INFO) + if not len(logger.handlers): + logger.addHandler(logging.StreamHandler(sys.stdout)) + logger.propagate = False # Do not propagate up to root logger, which may have other handlers + return logger diff --git a/easi_tools/load_s2l2a.py b/easi_tools/load_s2l2a.py new file mode 100644 index 0000000..918a1c7 --- /dev/null +++ b/easi_tools/load_s2l2a.py @@ -0,0 +1,293 @@ +#!python + +# Sentinel-2 L2A Collection 0 scaling and offset corrections. +# - Applies to data indexed from https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a +# - The newer https://earth-search.aws.element84.com/v1/collections/sentinel-2-c1-l2a (Collection 1) may not be affected in the same way +# +# TL;DR: +# DN values in COG files have different definitions depending on the processing baseline version +# and whether the offset change has been pre-applied by the cloud data custodian. +# +# Background: +# +# ESA has undertaken a reprocessing of the Sentinel-2 L2A product that includes +# a change to the offset value used to convert digital numbers (in file) to +# scientific values (reflectances). +# +# https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/level-2a-algorithms-products +# +# L2A algorithm and products: Starting with the PB 04.00 (25th January 2022), the dynamic +# range of the Level-2A products is shifted by a band-dependent constant: BOA_ADD_OFFSET. +# This offset will allow encoding negative surface reflectances that may occur over very +# dark surfaces. +# +# L2A_SRi = (L2A_DNi + BOA_ADD_OFFSETi) / QUANTIFICATION_VALUEi +# +# QUANTIFICATION_VALUEi = 10000 +# BOA_ADD_OFFSETi = -1000 +# +# refl = (dn -1000) / 10000 +# refl = dn/10000 - 1000/10000 +# refl = dn * 0.0001 - 0.1 +# +# These are the values in the EASI product definition, e.g. +# https://explorer.asia.easi-eo.solutions/products/s2_l2a.odc-product.yaml +# +# Example workflow: +# +# ESA's reprocessing is flowing through to the AWS open data repository of S2 L2A but +# while this stabilises we may see inconsistencies in time series queries due to: +# - More than one processed version of a dataset (scene) in the AWS bucket and indexed in an EASI database +# - Datasets (scenes) that indicate they have an offset applied by ESA but the offset correction +# has been not been applied to the COG +# +# Element-84 discussion: +# https://github.com/Element84/earth-search/issues/23#issuecomment-1834674853 + + +import xarray as xr +import pandas as pd +import logging +from pathlib import Path +import sys, re + +import datacube +from datacube.api.core import output_geobox +from datacube.api.query import SPATIAL_KEYS, CRS_KEYS, OTHER_KEYS +from datacube.utils import masking + + +# Set logger +log = logging.getLogger(Path(__file__).stem) +log.setLevel(logging.INFO) +log.addHandler(logging.StreamHandler(sys.stdout)) + +# Constants +search_keys = ( + 'product', + 'time', + 'geopolygon', + 'like', + 'limit', + 'ensure_location', + 'dataset_predicate', +) + SPATIAL_KEYS + CRS_KEYS + OTHER_KEYS + +# TODO: get measurement aliases from the ODC product record +refl_bands = { + 'coastal','band_01','B01','coastal_aerosol', + 'blue','band_02','B02', + 'green','band_03','B03', + 'red','band_04','B04', + 'rededge1','band_05','B05','red_edge_1', + 'rededge2','band_06','B06','red_edge_2', + 'rededge3','band_07','B07','red_edge_3', + 'nir','band_08','B08','nir_1', + 'nir08','band_8a','B8A','nir_2', + 'nir09','band_09','B09','nir_3', + 'swir16','band_11','B11','swir_1','swir_16', + 'swir22','band_12','B12','swir_2','swir_22', +} +scale_factor = 0.0001 +add_offset = -0.1 + + +def highest_sequence_number(matches: list) -> dict: + """Filter for the highest element84 processing sequence number per scene (scene label excluding the sequence number) + + : return : { scene_id_excluding_sequence_number : { highest_sequence_number : datacube.model.Dataset }} + """ + p = re.compile(r'(S2.+)_([0-9]+)_(L2A)') + sorter = {} + for ds in matches: + # Separate the scene label from the sequence number + label = ds.metadata_doc['label'] + m = p.match(label) + if not m: + log.warning(f'Dataset label does not match expected pattern: {label}') + continue + key = f'{m.group(1)}_{m.group(3)}' + seq = int(m.group(2)) + # Retain the highest sequence number + if key in sorter: + if list(sorter[key])[0] < seq: + sorter[key] = {seq: ds} + else: + sorter[key] = {seq: ds} + return sorter + + +def ds_requires_offset(ds: datacube.model.Dataset) -> bool: + """Return True if a dataset's metadata indicates that the offset correction should be applied""" + props = ds.metadata_doc['properties'] + + # If baseline is less than '04.00' then offset correction does not apply + baseline = props.get('s2:processing_baseline', '0.0') + p = re.compile(r'(\d+)\.(\d+)') + m = p.match(baseline) + if not m: + log.warning(f'Dataset processing_baseline does not match expected pattern: {baseline}') + return None + if int(m.group(1)) < 4: + return False + + # If the boa_offset_applied has been applied then offset correction is not required + boa_offset_applied = props.get('earthsearch:boa_offset_applied', False) + return not boa_offset_applied + + +def apply_correction_to_data(ds: xr.Dataset, offset: float = 0) -> xr.Dataset: + """Apply the scale and offset correction to each reflectance band where there is valid data (not nodata)""" + refl_vars = [x for x in ds.data_vars if x in refl_bands] + mask = masking.valid_data_mask(ds[refl_vars]) + # Save on a dask step? + if offset == 0: + ds[refl_vars] = ds[refl_vars].where(mask) * scale_factor + else: + ds[refl_vars] = ds[refl_vars].where(mask) * scale_factor + offset + return ds + + +def load_s2l2a_with_offset( + dc: datacube.Datacube, + query: dict, +) -> xr.Dataset: + """ + Replaces datacube.load(**query) for s2_l2a products. + + Method: + - Find all datasets matching the query (dc.find_datasets) + - Filter for the highest element84 processing sequence number per scene (scene label excluding the sequence number) + - Filter into two lists for datasets that have + - "s2:processing_baseline" >= "04.00" and "earthsearch:boa_offset_applied" == False (offset correction required) + - everything else (no correction required) + - If either list is empty then load the non-empty list, apply scale (and offset if required), and return the xarray Dataset + - Load and combine the two lists of datasets + - Load each list, apply scale (and offset if required) + - Concat on time dimension and sort by time + - Return the combined xarray Dataset + + Notes: + - Any 'groupby' function is applied to each of the xarray Datasets prior to them being combined. + This could create "extra" (non-grouped) time layers in the combined Dataset if the groupby function + would have grouped datasets (scenes) from both lists. + - Scale and offset are applied to the reflectance bands where there is valid data (not `nodata`). + This includes applying the "scale_factor" even if no datasets require the offset correction. + Other masks can be applied by the user (e.g. pixel quality or cloud masking). + """ + + product = query.get('product', '') + if product != 's2_l2a': + log.error(f'This function only applies to the "s2_l2a" product, not: {product}') + return None + + # Find all datasets matching the query + matches = None + if 'datasets' in query: + matches = query['datasets'] + del query['datasets'] + if matches is None: + search_params = {k:v for k,v in query.items() if k in search_keys} + matches = dc.find_datasets(**search_params) + if 'skip_broken_datasets' not in query: + # This helps to avoid data loading error messages + query['skip_broken_datasets'] = True + + # Filter for the highest element84 processing sequence number + sorter = highest_sequence_number(matches) + + # Filter into two lists + offset_applied, offset_required = [], [] + for key in sorter.keys(): + ds = list(sorter[key].values())[0] + isrequired = ds_requires_offset(ds) + if isrequired is None: + continue + elif isrequired: + offset_required.append(ds) + else: + offset_applied.append(ds) + matches_combined = offset_applied + offset_required + + # If either list is empty then no separation and merge is required + this_offset = None + if len(offset_applied) == 0: + log.info('All datasets require offset correction') + msg = 'The valid_data_mask, scale and offset have been applied to the reflectance bands' + this_offset = add_offset + if len(offset_required) == 0: + log.info('No datasets require offset correction') + msg = 'The valid_data_mask and scale (no offset) have been applied to the reflectance bands' + this_offset = 0 + if this_offset is not None: + data = dc.load( + datasets = matches_combined, + **query + ) + xx = apply_correction_to_data(data, this_offset) + log.info(msg) + return xx + + # DEBUG: What do we have + # def func(s): + # p = re.compile('(\d{8})') + # m = p.search(s[0]) + # if m: + # return m.group(1) + # log.info(f'Number of datasets in initial query: {len(matches)}') + # log.info(f'{sorted([(x.metadata_doc["label"],x.id) for x in matches], key=func)}') + # log.info(f'Number of datasets with offset applied: {len(offset_applied)}') + # log.info(f'{sorted([(x.metadata_doc["label"],x.id) for x in offset_applied], key=func)}') + # log.info(f'Number of datasets without offset applied: {len(offset_required)}') + # log.info(f'{sorted( [(x.metadata_doc["label"],x.id) for x in offset_required], key=func)}') + # return + + # Else, load data into two Datasets + log.info('Mix of datasets found with either offset required or not.') + log.info('We will load two xarrays, apply offset where required, and merge into one xarray.') + + # 1. Ensure the target geobox covers all datasets + target_geobox = output_geobox( + datasets = matches_combined, + **query, + ) + + # 2. Edit the query for our needs + # Ensure that dask time chunking = 1 + dask_input = None + if 'dask_chunks' in query: + dask_input = query['dask_chunks'] # Save + if dask_input.get('time', 1) != 1: + query['dask_chunks'].update({'time': 1}) + # Remove keys that are not compatible with 'like' + for x in ('output_crs', 'resolution', 'align'): + if x in query: + del query[x] + + # 3. Load two xarrays + data_offset_applied = dc.load( + datasets = offset_applied, + like = target_geobox, + **query + ) + data_offset_required = dc.load( + datasets = offset_required, + like = target_geobox, + **query + ) + + # 4. Apply respective scale and offsets + data_offset_applied = apply_correction_to_data(data_offset_applied) + data_offset_required = apply_correction_to_data(data_offset_required, add_offset) + + # 5. Combine the two xarrays + combined = xr.concat([data_offset_applied, data_offset_required], dim='time') + combined = combined.sortby('time') + + # 6. Reapply any time > 1 chunking + if dask_input is not None: + if dask_input.get('time', 1) != 1: + combined = combined.chunk(dask_input) + + log.info('The valid_data_mask, scale and offset have been applied to the reflectance bands') + return combined diff --git a/easi_tools/notebook_utils.py b/easi_tools/notebook_utils.py new file mode 100644 index 0000000..4f2ba13 --- /dev/null +++ b/easi_tools/notebook_utils.py @@ -0,0 +1,171 @@ +#!python3 + +# A collection of utilities that can be used in Python notebooks. +# +# License: Apache 2.0 + +# Created for EASI Hub training notebooks, https://dev.azure.com/csiro-easi/easi-hub-public/_git/hub-notebooks + +# Data tools +import numpy as np +import xarray as xr +import pandas as pd +import geopandas as gpd +import datacube +from datacube.utils import masking +from datetime import datetime + +# hvPlot, Holoviews, Datashader and Bokeh +import hvplot.pandas +import hvplot.xarray +import panel as pn +import holoviews as hv +# hv.extension("bokeh", logo=False) # Its likely set from in the notebooks + +# Jupyter Lab +from IPython.display import HTML + +# Python +import sys, os, re +import logging +from pathlib import Path +from collections import Counter +import contextlib + +# Dask +import dask +from dask.distributed import Client, LocalCluster +from dask_gateway import Gateway + +# EASIDefaults +from . import EasiDefaults + +# Set logger +logger = logging.getLogger(Path(__file__).stem) +logger.setLevel(logging.INFO) +if not len(logger.handlers): + logger.addHandler(logging.StreamHandler(sys.stdout)) + + +def display_table( + df: pd.DataFrame, + panel: bool = False, +): + """Display the full pandas dataframe. If panel is True use a panel object""" + table = None + if panel: + # Dicts are rendered as "[object Object]". Need to set a formatter, I guess. + table = pn.widgets.DataFrame(df, + # sizing_mode='stretch_width', # equal column widths, full screen + autosize_mode='fit_viewport', # fitted columns, about 90-95% width + # reorderable=True, # didn't work first try + ) + else: + with pd.option_context("display.max_rows", None, + "display.max_columns", None, + "display.max_colwidth", -1): + table = HTML( df.to_html().replace(r"\n", "
") ) + display(table) + + +def heading(txt: str): + """Print a simple HTML heading""" + display(HTML( f"

{txt}

" )) + + +def hv_table_hook(plot, element): + """Selected options for hv.table() formatting + + Use: df.hv.table().opts(hooks=[hv_table_hook]) + """ + plot.handles["table"].autosize_mode="fit_viewport" + # Other examples + # plot.handles['table'].row_height = 40 + # from bokeh.models.widgets import DateFormatter + # plot.handles['table'].columns[6].formatter = DateFormatter(format='%Y-%m-%d') + + +def xarray_object_size(data): + """Return a formatted string""" + val, unit = data.nbytes / (1024 ** 2), "MB" + if val > 1024: + val, unit = data.nbytes / (1024 ** 3), "GB" + return f"Dataset size: {val:.2f} {unit}" + + +def mostcommon_crs(dc, query): + """Adapted from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/Tools/dea_tools/datahandling.py""" + matching_datasets = dc.find_datasets(**query) + crs_list = [str(i.crs) for i in matching_datasets] + crs_mostcommon = None + if len(crs_list) > 0: + # Identify most common CRS + crs_counts = Counter(crs_list) + crs_mostcommon = crs_counts.most_common(1)[0][0] + else: + logger.warning("No data was found for the supplied product query") + return crs_mostcommon + + +def initialize_dask(use_gateway=False, workers=(1,2), wait=False, local_port=8786, **kwargs): + """Initialize a Dask Gateway or Local cluster""" + # Check inputs + if isinstance(workers, (int, float)): + workers = (int(workers), int(workers)) + if len(workers) != 2: + logger.error("Require workers to be a single integer or a 2-element tuple/list") + return None, None + if isinstance(local_port, (str, float)): + local_port = int(local_port) + + # Dask gateway + if use_gateway: + gateway = Gateway() + clusters = gateway.list_clusters() + if not clusters: + logger.info("Starting new cluster") + cluster = gateway.new_cluster(**kwargs) + else: + logger.info(f"An existing cluster was found. Connecting to: {clusters[0].name}") + cluster = gateway.connect(clusters[0].name) + client = cluster.get_client() + cluster.adapt(minimum=workers[0], maximum=workers[1]) + if wait: + logger.info("Waiting for at least one cluster worker") + # client.wait_for_workers(n_workers=1) # Before release 2023.10.0 + client.sync(client._wait_for_workers,n_workers=1) # Since release 2023.10.0 + + # Local cluster + else: + cluster = LocalCluster(n_workers=4) + client = Client(cluster) + server = f'https://hub.{EasiDefaults().domain}' # Or replace if not using EasiDefaults + user = os.environ.get('JUPYTERHUB_SERVICE_PREFIX') # Current user + dask.config.set({"distributed.dashboard.link": f'{server}{user}' + "proxy/{port}/status"}) # port is evaluated by dask + + return cluster, client + + +def localcluster_dashboard(client, server="https://hub.csiro.easi-eo.solutions"): + """Return a dashboard link using jupyter proxy""" + dashboard_link = client.dashboard_link + for host in ("127.0.0.1", "localhost"): + if host in dashboard_link: + port = re.search(r":(\d+)\/status", dashboard_link).group(1) + dashboard_link = f'{server}{os.environ["JUPYTERHUB_SERVICE_PREFIX"]}proxy/{port}/status' + break + return dashboard_link + + +@contextlib.contextmanager +def unset_cachingproxy(): + """Unset the EASI caching proxy with a context manager""" + # Inspired by https://stackoverflow.com/a/34333710 + env = os.environ + remove = ("AWS_HTTPS", "GDAL_HTTP_PROXY") + update_after = {k: env[k] for k in remove} + try: + [env.pop(k, None) for k in remove] + yield + finally: + env.update(update_after) diff --git a/new_import_ODC.py b/new_import_ODC.py index 0845f48..b0dfb97 100644 --- a/new_import_ODC.py +++ b/new_import_ODC.py @@ -82,47 +82,87 @@ import joblib def load_data(dc, date_range, longtitude_range, latitude_range): + """ + Load Sentinel-2 L2A data using direct datacube.load() + without spatial filtering (which was causing 0 results). + + Note: Data is loaded in UTM (EPSG:32648) to avoid CRS issues. + Spatial filtering on lat/lon is skipped to return maximum data. + """ product = 's2_l2a' - query = { - 'product': product, # Product name - 'x': longtitude_range, # "x" axis bounds - 'y': latitude_range, # "y" axis bounds - 'time': date_range, # Any parsable date strings - } - native_crs = notebook_utils.mostcommon_crs(dc, query) - print(f'Most common native CRS: {native_crs}') + native_crs = 'EPSG:32648' # UTM Zone 48N for Vietnam measurements = ['red', 'nir', 'scl'] - - load_params = { - 'measurements': measurements, # Selected measurement or alias names - 'output_crs': native_crs, # Target EPSG code - 'resolution': (-10, 10), # Target resolution - 'group_by': 'solar_day', # Scene grouping - 'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks - } - data = load_s2l2a_with_offset( - dc, - query | load_params # Combine the two dicts that contain our search and load parameters - ) - return data + + print(f'Loading Sentinel-2 data (EPSG:32648)...') + print(f' Time range: {date_range}') + print(f' Measurements: {measurements}') + + try: + # Load ALL available data WITHOUT dask_chunks (forces immediate load) + # This avoids the metadata issue with dc.load() when using dask_chunks + data = dc.load( + product=product, + time=date_range, + measurements=measurements, + output_crs=native_crs, + resolution=(-10, 10), + group_by='solar_day', + skip_broken_datasets=True + ) + + print(f'✅ Data loaded successfully!') + print(f' Dimensions: {dict(data.sizes)}') + print(f' Time steps: {len(data.time)}') + print(f' Spatial extent: x={len(data.x)}, y={len(data.y)}') + print(f' Data type: numpy arrays (not Dask)') + + return data + + except Exception as e: + print(f'❌ Error loading data: {e}') + import traceback + traceback.print_exc() + return None def mask_clean(data): - flag_name = 'scl' - flag_desc = masking.describe_variable_flags(data[flag_name]) # Pandas dataframe - display(flag_desc) - display(flag_desc.loc['qa'].values[1]) - # Create a "data quality" Mask layer - flags_def = flag_desc.loc['qa'].values[1] - good_pixel_flags = [flags_def[str(i)] for i in [2, 4, 5, 6]] # To pass strings to enum_to_bool() - - # enum_to_bool calculates the pixel-wise "or" of each set of pixels given by good_pixel_flags - # 1 = good data - # 0 = "bad" data - good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags) + """ + Clean data by masking clouds and bad pixels using the SCL (Scene Classification Layer). + + SCL classes: + - 0: No Data + - 1: Saturated/Defective + - 2: Dark Area Pixels + - 3: Cloud Shadows + - 4: Vegetation ✓ GOOD + - 5: Not Vegetated ✓ GOOD + - 6: Water ✓ GOOD + - 7: Unclassified ✓ GOOD + - 8: Cloud Medium Probability ✗ BAD + - 9: Cloud High Probability ✗ BAD + - 10: Thin Cirrus ✗ BAD + - 11: Snow/Ice ✗ BAD + """ + + # Good pixel classes (keep these) + good_pixel_classes = [4, 5, 6, 7] + + # Create mask: 1 where SCL is in good_pixel_classes, 0 otherwise + good_pixel_mask = data['scl'].isin(good_pixel_classes) + + print(f'✅ Cloud masking applied') + print(f' Good pixel classes: {good_pixel_classes}') + print(f' Mask created (dask-backed, not yet computed)') + + # Get all variables except SCL data_layer_names = [x for x in data.data_vars if x != 'scl'] - # Apply good pixel mask to blue, green, red and nir. + + # Apply mask to all layers result = data[data_layer_names].where(good_pixel_mask).persist() + + print(f' Data variables masked: {data_layer_names}') + print(f' Result persisted to workers') + return result @@ -360,7 +400,17 @@ def load_data_sen2(dc, date_range, coordinates): 'y': latitude_range, # "y" axis bounds 'time': date_range, # Any parsable date strings } - native_crs = notebook_utils.mostcommon_crs(dc, query) + + # Try to get native CRS, default to EPSG:32648 (UTM Zone 48N) for Vietnam + try: + native_crs = notebook_utils.mostcommon_crs(dc, query) + if native_crs is None: + print('⚠️ Could not determine native CRS, using EPSG:32648 (UTM Zone 48N)') + native_crs = 'EPSG:32648' + except Exception as e: + print(f'⚠️ Error determining CRS: {e}, using EPSG:32648') + native_crs = 'EPSG:32648' + print(f'Most common native CRS: {native_crs}') # measurements = ['red','green', 'blue', 'nir', 'scl'] @@ -373,10 +423,27 @@ def load_data_sen2(dc, date_range, coordinates): 'group_by': 'solar_day', # Scene grouping 'dask_chunks': {'x': 2048, 'y': 2048}, # Dask chunks } - data = load_s2l2a_with_offset( - dc, - query | load_params # Combine the two dicts that contain our search and load parameters - ) + + try: + data = load_s2l2a_with_offset( + dc, + query | load_params # Combine the two dicts that contain our search and load parameters + ) + except Exception as e: + print(f'❌ Error loading data: {e}') + print('Attempting direct dc.load without offset correction...') + data = dc.load( + product=product, + x=longtitude_range, + y=latitude_range, + time=date_range, + measurements=measurements, + output_crs=native_crs, + resolution=(-10, 10), + group_by='solar_day', + dask_chunks={'x': 2048, 'y': 2048}, + skip_broken_datasets=True + ) return data def mask_cloud(data):