Source code for dea_tools.datahandling

# dea_datahandling.py
"""
Loading and manipulating Digital Earth Australia products and data
using the Open Data Cube and xarray.

License: The code in this notebook is licensed under the Apache License,
Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0). Digital Earth
Australia data is licensed under the Creative Commons by Attribution 4.0
license (https://creativecommons.org/licenses/by/4.0/).

Contact: If you need assistance, please post a question on the Open Data
Cube Slack channel (http://slack.opendatacube.org/) or on the GIS Stack
Exchange (https://gis.stackexchange.com/questions/ask?tags=open-data-cube)
using the `open-data-cube` tag (you can view previously asked questions
here: https://gis.stackexchange.com/questions/tagged/open-data-cube).

If you would like to report an issue with this script, you can file one
on Github (https://github.com/GeoscienceAustralia/dea-notebooks/issues/new).

Last modified: Feb 2024
"""

import datetime

# Import required packages
import os
import warnings
import zipfile
import requests
from collections import Counter

import rioxarray
import numpy as np
import pandas as pd
import xarray as xr
import sklearn.decomposition
from scipy.ndimage import binary_dilation
from skimage.color import hsv2rgb, rgb2hsv
from skimage.exposure import match_histograms

import odc.geo.xr
import odc.algo
from odc.algo import mask_cleanup
from datacube.utils.dates import normalise_dt


def _dc_query_only(**kw):
    """
    Remove load-only datacube parameters, the rest can be
    passed to Query/dc.find_datasets.

    Returns
    -------
    dict of query parameters
    """

    def _impl(
        measurements=None,
        output_crs=None,
        resolution=None,
        resampling=None,
        skip_broken_datasets=None,
        dask_chunks=None,
        fuse_func=None,
        align=None,
        datasets=None,
        progress_cbk=None,
        group_by=None,
        **query,
    ):
        return query

    return _impl(**kw)


def _common_bands(dc, products):
    """
    Takes a list of products and returns a list of measurements/bands
    that are present in all products

    Returns
    -------
    List of band names
    """
    common = None
    bands = None

    for p in products:
        p = dc.index.products.get_by_name(p)
        if common is None:
            common = set(p.measurements)
            bands = list(p.measurements)
        else:
            common = common.intersection(set(p.measurements))
    return [band for band in bands if band in common]


[docs] def load_ard( dc, products=None, cloud_mask="fmask", min_gooddata=0.00, mask_pixel_quality=True, mask_filters=None, mask_contiguity=False, fmask_categories=["valid", "snow", "water"], s2cloudless_categories=["valid"], ls7_slc_off=True, dtype="auto", predicate=None, **kwargs, ): """ Load multiple Geoscience Australia Landsat or Sentinel 2 Collection 3 products (e.g. Landsat 5, 7, 8, 9; Sentinel 2A and 2B), optionally apply pixel quality/cloud masking and contiguity masks, and drop time steps that contain greater than a minimum proportion of good quality (e.g. non-cloudy or shadowed) pixels. The function supports loading the following Landsat products: * ga_ls5t_ard_3 * ga_ls7e_ard_3 * ga_ls8c_ard_3 * ga_ls9c_ard_3 And Sentinel-2 products: * ga_s2am_ard_3 * ga_s2bm_ard_3 Cloud masking can be performed using the Fmask (Function of Mask) cloud mask for Landsat and Sentinel-2, and the s2cloudless (Sentinel Hub cloud detector for Sentinel-2 imagery) cloud mask for Sentinel-2. Last modified: June 2023 Parameters ---------- dc : datacube Datacube object The Datacube to connect to, i.e. ``dc = datacube.Datacube()``. This allows you to also use development datacubes if required. products : list A list of product names to load. Valid options are ['ga_ls5t_ard_3', 'ga_ls7e_ard_3', 'ga_ls8c_ard_3', 'ga_ls9c_ard_3'] for Landsat, ['ga_s2am_ard_3', 'ga_s2bm_ard_3'] for Sentinel 2. cloud_mask : string, optional The cloud mask used by the function. This is used for both masking out poor quality pixels (e.g. clouds) if ``mask_pixel_quality=True``, and for calculating the ``min_gooddata`` percentage when dropping cloudy or low quality satellite observations. Two cloud masks are supported: * 'fmask': (default; available for Landsat, Sentinel-2) * 's2cloudless' (Sentinel-2 only) min_gooddata : float, optional The minimum percentage of good quality pixels required for a satellite observation to be loaded. Defaults to 0.00 which will return all observations regardless of pixel quality (set to e.g. 0.99 to return only observations with more than 99% good quality pixels). mask_pixel_quality : str or bool, optional Whether to mask out poor quality (e.g. cloudy) pixels by setting them as nodata. Depending on the choice of cloud mask, the function will identify good quality pixels using the categories passed to the ``fmask_categories`` or ``s2cloudless_categories`` params. Set to False to turn off pixel quality masking completely. Poor quality pixels will be set to NaN (and convert all data to `float32`) if ``dtype='auto'``, or be set to the data's native nodata value (usually -999) if ``dtype='native'`` (see 'dtype' below for more details). mask_filters : iterable of tuples, optional Iterable tuples of morphological operations - ("<operation>", <radius>) to apply to the inverted pixel quality mask, where: operation: string; one of these morphological operations: * ``'dilation'`` = Expands poor quality pixels/clouds outwards * ``'erosion'`` = Shrinks poor quality pixels/clouds inwards * ``'closing'`` = Remove small holes in clouds by expanding then shrinking poor quality pixels * ``'opening'`` = Remove small or narrow clouds by shrinking then expanding poor quality pixels radius: int e.g. ``mask_filters=[('erosion', 5), ("opening", 2), ("dilation", 2)]`` mask_contiguity : str or bool, optional Whether to mask out pixels that are missing data in any band (i.e. "non-contiguous" pixels). This can be important for generating clean composite datasets. The default of False will not apply any contiguity mask. If loading NBART data, set: * ``mask_contiguity='nbart'`` (or ``mask_contiguity=True``) If loading NBAR data, specify: * ``mask_contiguity='nbar'`` Non-contiguous pixels will be set to NaN if `dtype='auto'`, or set to the data's native nodata value if `dtype='native'` (see 'dtype' below). fmask_categories : list, optional A list of Fmask cloud mask categories to consider as good quality pixels when calculating `min_gooddata` and when masking data by pixel quality if ``mask_pixel_quality=True``. The default is ``['valid', 'snow', 'water']``; all other Fmask categories ('cloud', 'shadow', 'nodata') will be treated as low quality pixels. Choose from: 'nodata', 'valid', 'cloud', 'shadow', 'snow', and 'water'. s2cloudless_categories : list, optional A list of s2cloudless cloud mask categories to consider as good quality pixels when calculating `min_gooddata` and when masking data by pixel quality if ``mask_pixel_quality=True``. The default is `['valid']`; all other s2cloudless categories ('cloud', 'nodata') will be treated as low quality pixels. Choose from: 'nodata', 'valid', or 'cloud'. ls7_slc_off : bool, optional An optional boolean indicating whether to include data from after the Landsat 7 SLC failure (i.e. SLC-off). Defaults to True, which keeps all Landsat 7 observations > May 31 2003. dtype : string, optional Controls the data type/dtype that layers are coerced to after loading. Valid values: 'native', 'auto', and 'float{16|32|64}'. When 'auto' is used, the data will be converted to `float32` if masking is used, otherwise data will be returned in the native data type of the data. Be aware that if data is loaded in its native dtype, nodata and masked pixels will be returned with the data's native nodata value (typically -999), not NaN. predicate : function, optional DEPRECATED: Please use `dataset_predicate` instead. An optional function that can be passed in to restrict the datasets that are loaded. A predicate function should take a `datacube.model.Dataset` object as an input (i.e. as returned from `dc.find_datasets`), and return a boolean. For example, a predicate function could be used to return True for only datasets acquired in January: `dataset.time.begin.month == 1` **kwargs : A set of keyword arguments to `dc.load` that define the spatiotemporal query and load parameters used to extract data. Keyword arguments can either be listed directly in the `load_ard` call like any other parameter (e.g. `measurements=['nbart_red']`), or by passing in a query kwarg dictionary (e.g. `**query`). Keywords can include `measurements`, `x`, `y`, `time`, `resolution`, `resampling`, `group_by`, `crs`; see the `dc.load` documentation for all possible options: https://datacube-core.readthedocs.io/en/latest/api/indexed-data/generate/datacube.Datacube.load.html Returns ------- combined_ds : xarray.Dataset An xarray.Dataset containing only satellite observations with a proportion of good quality pixels greater than `min_gooddata`. Notes ----- The `load_ard` function builds on the Open Data Cube's native `dc.load` function by adding the ability to load multiple satellite data products at once, and automatically apply cloud masking and filtering. For loading non-satellite data products (e.g. DEA Water Observations), use `dc.load` instead. """ ######### # Setup # ######### # Verify that products were provided if not products: raise ValueError( "Please provide a list of product names to load data from. " "Valid options are: ['ga_ls5t_ard_3', 'ga_ls7e_ard_3', " "'ga_ls8c_ard_3', 'ga_ls9c_ard_3'] for Landsat, and " "['ga_s2am_ard_3', 'ga_s2bm_ard_3'] for Sentinel 2." ) # Determine whether products are all Landsat, all S2, or mixed elif all(["ls" in product for product in products]): product_type = "ls" elif all(["s2" in product for product in products]): product_type = "s2" else: product_type = "mixed" warnings.warn( "You have selected a combination of Landsat and Sentinel-2 " "products. This can produce unexpected results as these " "products use the same names for different spectral bands " "(e.g. Landsat and Sentinel-2's 'nbart_swir_2'); use with " "caution." ) # Set contiguity band depending on `mask_contiguity`; # "oa_nbart_contiguity" if True, False or "nbart", # "oa_nbar_contiguity" if "nbar" if mask_contiguity in (True, False, "nbart"): contiguity_band = "oa_nbart_contiguity" elif mask_contiguity == "nbar": contiguity_band = "oa_nbar_contiguity" else: raise ValueError( f"Unsupported value '{mask_contiguity}' passed to " "`mask_contiguity`. Please provide either 'nbart', 'nbar', " "True, or False." ) # Set pixel quality (PQ) band depending on `cloud_mask` if cloud_mask == "fmask": pq_band = "oa_fmask" pq_categories = fmask_categories elif cloud_mask == "s2cloudless": pq_band = "oa_s2cloudless_mask" pq_categories = s2cloudless_categories # Raise error if s2cloudless is requested for Landsat products if product_type in ["ls", "mixed"]: raise ValueError( "The 's2cloudless' cloud mask is not available for " "Landsat products. Please set `mask_pixel_quality` to " "'fmask' or False." ) else: raise ValueError( f"Unsupported value '{cloud_mask}' passed to " "`cloud_mask`. Please provide either 'fmask', " "'s2cloudless', True, or False." ) # To ensure that the categorical PQ/contiguity masking bands are # loaded using nearest neighbour resampling, we need to add these to # the resampling kwarg if it exists and is not "nearest". # This only applies if a string resampling method is supplied; # if a resampling dictionary (e.g. `resampling={'*': 'bilinear', # 'oa_fmask': 'mode'}` is passed instead we assume the user wants # to select custom resampling methods for each of their bands. resampling = kwargs.get("resampling", None) if isinstance(resampling, str) and resampling not in (None, "nearest"): kwargs["resampling"] = { "*": resampling, pq_band: "nearest", contiguity_band: "nearest", } # We extract and deal with `dask_chunks` separately as every # function call uses dask internally regardless of whether the user # sets `dask_chunks` themselves dask_chunks = kwargs.pop("dask_chunks", None) # Create a list of requested measurements so that we can eventually # return only the measurements the user orignally asked for requested_measurements = kwargs.pop("measurements", None) # Copy our measurements list so we can temporarily append extra PQ # and/or contiguity masking bands when loading our data measurements = requested_measurements.copy() if requested_measurements else None # Deal with "load all" case: pick a set of bands that are common # across requested products if measurements is None: measurements = _common_bands(dc, products) # Deal with edge case where user supplies alias for PQ/contiguity # by stripping PQ/contiguity masks of their "oa_" prefix else: contiguity_band = ( contiguity_band.replace("oa_", "") if contiguity_band.replace("oa_", "") in measurements else contiguity_band ) pq_band = ( pq_band.replace("oa_", "") if pq_band.replace("oa_", "") in measurements else pq_band ) # If `measurements` are specified but do not include PQ or # contiguity variables, add these to `measurements` if pq_band not in measurements: measurements.append(pq_band) if mask_contiguity and contiguity_band not in measurements: measurements.append(contiguity_band) # Get list of data and mask bands so that we can later exclude # mask bands from being masked themselves data_bands = [ band for band in measurements if band not in (pq_band, contiguity_band) ] mask_bands = [band for band in measurements if band not in data_bands] ################# # Find datasets # ################# # Pull out query params only to pass to dc.find_datasets query = _dc_query_only(**kwargs) # If predicate is specified, use this function to filter the list # of datasets prior to load if predicate: print( "The 'predicate' parameter will be deprecated in future " "versions of this function as this functionality has now " "been added to Datacube itself. Please use " "`dataset_predicate=...` instead." ) query["dataset_predicate"] = predicate # Extract list of datasets for each product using query params dataset_list = [] # Get list of datasets for each product print("Finding datasets") for product in products: # Obtain list of datasets for product print( f" {product} (ignoring SLC-off observations)" if not ls7_slc_off and product == "ga_ls7e_ard_3" else f" {product}" ) datasets = dc.find_datasets(product=product, **query) # Remove Landsat 7 SLC-off observations if ls7_slc_off=False if not ls7_slc_off and product == "ga_ls7e_ard_3": datasets = [ i for i in datasets if normalise_dt(i.time.begin) < datetime.datetime(2003, 5, 31) ] # Add any returned datasets to list dataset_list.extend(datasets) # Raise exception if no datasets are returned if len(dataset_list) == 0: raise ValueError( "No data available for query: ensure that " "the products specified have data for the " "time and location requested" ) ############# # Load data # ############# # Note we always load using dask here so that we can lazy load data # before filtering by `min_gooddata` ds = dc.load( datasets=dataset_list, measurements=measurements, dask_chunks={} if dask_chunks is None else dask_chunks, **kwargs, ) #################### # Filter good data # #################### # Calculate pixel quality mask pq_mask = odc.algo.fmask_to_bool(ds[pq_band], categories=pq_categories) # The good data percentage calculation has to load all pixel quality # data, which can be slow. If the user has chosen no filtering # by using the default `min_gooddata = 0`, we can skip this step # completely to save processing time if min_gooddata > 0.0: # Compute good data for each observation as % of total pixels print(f"Counting good quality pixels for each time step using {cloud_mask}") data_perc = pq_mask.sum(axis=[1, 2], dtype="int32") / ( pq_mask.shape[1] * pq_mask.shape[2] ) keep = (data_perc >= min_gooddata).persist() # Filter by `min_gooddata` to drop low quality observations total_obs = len(ds.time) ds = ds.sel(time=keep) pq_mask = pq_mask.sel(time=keep) print( f"Filtering to {len(ds.time)} out of {total_obs} " f"time steps with at least {min_gooddata:.1%} " f"good quality pixels" ) # Morphological filtering on cloud masks if (mask_filters is not None) & (mask_pixel_quality != False): print(f"Applying morphological filters to pixel quality mask: {mask_filters}") pq_mask = ~mask_cleanup(~pq_mask, mask_filters=mask_filters) warnings.warn( "As of `dea_tools` v0.3.0, pixel quality masks are " "inverted before being passed to `mask_filters` (i.e. so " "that good quality/clear pixels are False and poor quality " "pixels/clouds are True). This means that 'dilation' will " "now expand cloudy pixels, rather than shrink them as in " "previous versions." ) ############### # Apply masks # ############### # Create a combined mask to hold both pixel quality and contiguity. # This is more efficient than creating multiple dask tasks for # similar masking operations. mask = None # Add pixel quality mask to combined mask if mask_pixel_quality: print(f"Applying {cloud_mask} pixel quality/cloud mask") mask = pq_mask # Add contiguity mask to combined mask if mask_contiguity: print(f"Applying contiguity mask ({contiguity_band})") cont_mask = ds[contiguity_band] == 1 # If mask already has data if mask_pixel_quality == True, # multiply with cont_mask to perform a logical 'or' operation # (keeping only pixels good in both) mask = cont_mask if mask is None else mask * cont_mask # Split into data/masks bands, as conversion to float and masking # should only be applied to data bands ds_data = ds[data_bands] ds_masks = ds[mask_bands] # Mask data if either of the above masks were generated if mask is not None: ds_data = odc.algo.keep_good_only(ds_data, where=mask) # Automatically set dtype to either native or float32 depending # on whether masking was requested if dtype == "auto": dtype = "native" if mask is None else "float32" # Set nodata values using odc.algo tools to reduce peak memory # use when converting data dtype if dtype != "native": ds_data = odc.algo.to_float(ds_data, dtype=dtype) # Put data and mask bands back together attrs = ds.attrs ds = xr.merge([ds_data, ds_masks]) ds.attrs.update(attrs) ############### # Return data # ############### # Drop bands not originally requested by user if requested_measurements: ds = ds[requested_measurements] # If user supplied `dask_chunks`, return data as a dask array # without actually loading it into memory if dask_chunks is not None: print(f"Returning {len(ds.time)} time steps as a dask array") return ds else: print(f"Loading {len(ds.time)} time steps") return ds.compute()
[docs] def mostcommon_crs(dc, product, query): """ Takes a given query and returns the most common CRS for observations returned for that spatial extent. This can be useful when your study area lies on the boundary of two UTM zones, forcing you to decide which CRS to use for your `output_crs` in `dc.load`. Parameters ---------- dc : datacube Datacube object The Datacube to connect to, i.e. `dc = datacube.Datacube()`. This allows you to also use development datacubes if required. product : str A product name (or list of product names) to load CRSs from. query : dict A datacube query including x, y and time range to assess for the most common CRS Returns ------- epsg_string : str An EPSG string giving the most common CRS from all datasets returned by the query above """ # Find list of datasets matching query for either product or # list of products if isinstance(product, list): matching_datasets = [] for i in product: matching_datasets.extend(dc.find_datasets(product=i, **query)) else: matching_datasets = dc.find_datasets(product=product, **query) # Extract all CRSs crs_list = [str(i.crs) for i in matching_datasets] # If CRSs are returned if len(crs_list) > 0: # Identify most common CRS crs_counts = Counter(crs_list) crs_mostcommon = crs_counts.most_common(1)[0][0] # Warn user if multiple CRSs are encountered if len(crs_counts.keys()) > 1: warnings.warn( f"Multiple UTM zones {list(crs_counts.keys())} " f"were returned for this query. Defaulting to " f"the most common zone: {crs_mostcommon}", UserWarning, ) return crs_mostcommon else: raise ValueError( f"No CRS was returned as no data was found for " f"the supplied product ({product}) and query. " f"Please ensure that data is available for " f"{product} for the spatial extents and time " f"period specified in the query (e.g. by using " f"the Data Cube Explorer for this datacube " f"instance)." )
[docs] def download_unzip(url, output_dir=None, remove_zip=True): """ Downloads and unzips a .zip file from an external URL to a local directory. Parameters ---------- url : str A string giving a URL path to the zip file you wish to download and unzip output_dir : str, optional An optional string giving the directory to unzip files into. Defaults to None, which will unzip files in the current working directory remove_zip : bool, optional An optional boolean indicating whether to remove the downloaded .zip file after files are unzipped. Defaults to True, which will delete the .zip file. """ # Get basename for zip file zip_name = os.path.basename(url) # Raise exception if the file is not of type .zip if not zip_name.endswith(".zip"): raise ValueError( f"The URL provided does not point to a .zip " f"file (e.g. {zip_name}). Please specify a " f"URL path to a valid .zip file" ) # Download zip file print(f"Downloading {zip_name}") r = requests.get(url) with open(zip_name, "wb") as f: f.write(r.content) # Extract into output_dir with zipfile.ZipFile(zip_name, "r") as zip_ref: zip_ref.extractall(output_dir) print( f"Unzipping output files to: " f"{output_dir if output_dir else os.getcwd()}" ) # Optionally cleanup if remove_zip: os.remove(zip_name)
[docs] def wofs_fuser(dest, src): """ Fuse two WOfS water measurements represented as ``ndarray``s. Note: this is a copy of the function located here: https://github.com/GeoscienceAustralia/digitalearthau/blob/develop/digitalearthau/utils.py """ empty = (dest & 1).astype(bool) both = ~empty & ~((src & 1).astype(bool)) dest[empty] = src[empty] dest[both] |= src[both]
[docs] def dilate(array, dilation=10, invert=True): """ Dilate a binary array by a specified nummber of pixels using a disk-like radial dilation. By default, invalid (e.g. False or 0) values are dilated. This is suitable for applications such as cloud masking (e.g. creating a buffer around cloudy or shadowed pixels). This functionality can be reversed by specifying `invert=False`. Parameters ---------- array : array The binary array to dilate. dilation : int, optional An optional integer specifying the number of pixels to dilate by. Defaults to 10, which will dilate `array` by 10 pixels. invert : bool, optional An optional boolean specifying whether to invert the binary array prior to dilation. The default is True, which dilates the invalid values in the array (e.g. False or 0 values). Returns ------- An array of the same shape as `array`, with valid data pixels dilated by the number of pixels specified by `dilation`. """ y, x = np.ogrid[ -dilation : (dilation + 1), -dilation : (dilation + 1), ] # disk-like radial dilation kernel = (x * x) + (y * y) <= (dilation + 0.5) ** 2 # If invert=True, invert True values to False etc if invert: array = ~array return ~binary_dilation( array.astype(bool), structure=kernel.reshape((1,) + kernel.shape) )
[docs] def paths_to_datetimeindex(paths, string_slice=(0, 10)): """ Helper function to generate a Pandas datetimeindex object from dates contained in a file path string. Parameters ---------- paths : list of strings A list of file path strings that will be used to extract times string_slice : tuple An optional tuple giving the start and stop position that contains the time information in the provided paths. These are applied to the basename (i.e. file name) in each path, not the path itself. Defaults to (0, 10). Returns ------- datetime : pandas.DatetimeIndex A pandas.DatetimeIndex object containing a 'datetime64[ns]' derived from the file paths provided by `paths`. """ date_strings = [os.path.basename(i)[slice(*string_slice)] for i in paths] return pd.to_datetime(date_strings)
def _select_along_axis(values, idx, axis): other_ind = np.ix_(*[np.arange(s) for s in idx.shape]) sl = other_ind[:axis] + (idx,) + other_ind[axis:] return values[sl]
[docs] def first(array: xr.DataArray, dim: str, index_name: str = None) -> xr.DataArray: """ Finds the first occuring non-null value along the given dimension. Parameters ---------- array : xr.DataArray The array to search. dim : str The name of the dimension to reduce by finding the first non-null value. Returns ------- reduced : xr.DataArray An array of the first non-null values. The `dim` dimension will be removed, and replaced with a coord of the same name, containing the value of that dimension where the last value was found. """ axis = array.get_axis_num(dim) idx_first = np.argmax(~pd.isnull(array), axis=axis) reduced = array.reduce(_select_along_axis, idx=idx_first, axis=axis) reduced[dim] = array[dim].isel({dim: xr.DataArray(idx_first, dims=reduced.dims)}) if index_name is not None: reduced[index_name] = xr.DataArray(idx_first, dims=reduced.dims) return reduced
[docs] def last(array: xr.DataArray, dim: str, index_name: str = None) -> xr.DataArray: """ Finds the last occuring non-null value along the given dimension. Parameters ---------- array : xr.DataArray The array to search. dim : str The name of the dimension to reduce by finding the last non-null value. index_name : str, optional If given, the name of a coordinate to be added containing the index of where on the dimension the nearest value was found. Returns ------- reduced : xr.DataArray An array of the last non-null values. The `dim` dimension will be removed, and replaced with a coord of the same name, containing the value of that dimension where the last value was found. """ axis = array.get_axis_num(dim) rev = (slice(None),) * axis + (slice(None, None, -1),) idx_last = -1 - np.argmax(~pd.isnull(array)[rev], axis=axis) reduced = array.reduce(_select_along_axis, idx=idx_last, axis=axis) reduced[dim] = array[dim].isel({dim: xr.DataArray(idx_last, dims=reduced.dims)}) if index_name is not None: reduced[index_name] = xr.DataArray(idx_last, dims=reduced.dims) return reduced
[docs] def nearest( array: xr.DataArray, dim: str, target, index_name: str = None ) -> xr.DataArray: """ Finds the nearest values to a target label along the given dimension, for all other dimensions. E.g. For a DataArray with dimensions ('time', 'x', 'y') nearest_array = nearest(array, 'time', '2017-03-12') will return an array with the dimensions ('x', 'y'), with non-null values found closest for each (x, y) pixel to that location along the time dimension. The returned array will include the 'time' coordinate for each x,y pixel that the nearest value was found. Parameters ---------- array : xr.DataArray The array to search. dim : str The name of the dimension to look for the target label. target : same type as array[dim] The value to look up along the given dimension. index_name : str, optional If given, the name of a coordinate to be added containing the index of where on the dimension the nearest value was found. Returns ------- nearest_array : xr.DataArray An array of the nearest non-null values to the target label. The `dim` dimension will be removed, and replaced with a coord of the same name, containing the value of that dimension closest to the given target label. """ before_target = slice(None, target) after_target = slice(target, None) da_before = array.sel({dim: before_target}) da_after = array.sel({dim: after_target}) da_before = last(da_before, dim, index_name) if da_before[dim].shape[0] else None da_after = first(da_after, dim, index_name) if da_after[dim].shape[0] else None if da_before is None and da_after is not None: return da_after if da_after is None and da_before is not None: return da_before target = array[dim].dtype.type(target) is_before_closer = abs(target - da_before[dim]) < abs(target - da_after[dim]) nearest_array = xr.where(is_before_closer, da_before, da_after, keep_attrs=True) nearest_array[dim] = xr.where( is_before_closer, da_before[dim], da_after[dim], keep_attrs=True ) if index_name is not None: nearest_array[index_name] = xr.where( is_before_closer, da_before[index_name], da_after[index_name], keep_attrs=True, ) return nearest_array
[docs] def parallel_apply(ds, dim, func, use_threads=False, *args, **kwargs): """ Applies a custom function in parallel along the dimension of an xarray.Dataset or xarray.DataArray. The function can be any function that can be applied to an individual xarray.Dataset or xarray.DataArray (e.g. data for a single timestep). The function should also return data in xarray.Dataset or xarray.DataArray format. This function is useful as a simple method for parallising code that cannot easily be parallised using Dask. Parameters ---------- ds : xarray.Dataset or xarray.DataArray xarray data with a dimension `dim` to apply the custom function along. dim : string The dimension along which the custom function will be applied. func : function The function that will be applied in parallel to each array along dimension `dim`. The first argument passed to this function should be the array along `dim`. use_threads : bool, optional Whether to use threads instead of processes for parallelisation. Defaults to False, which means it'll use multi-processing. In brief, the difference between threads and processes is that threads share memory, while processes have separate memory. *args : Any number of arguments that will be passed to `func`. **kwargs : Any number of keyword arguments that will be passed to `func`. Returns ------- xarray.Dataset A concatenated dataset containing an output for each array along the input `dim` dimension. """ from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor from functools import partial from itertools import repeat from tqdm import tqdm # Use threads or processes if use_threads: Executor = ThreadPoolExecutor else: Executor = ProcessPoolExecutor with Executor as executor: # Update func to add kwargs func = partial(func, **kwargs) # Apply func in parallel groups = [group for (i, group) in ds.groupby(dim)] to_iterate = (groups, *(repeat(i, len(groups)) for i in args)) out_list = list(tqdm(executor.map(func, *to_iterate), total=len(groups))) # Combine to match the original dataset return xr.concat(out_list, dim=ds[dim])
def _apply_weights(da, band_weights): """ Apply weights from a dictionary to the bands of a multispectral xarray.DataArray. Raises a ValueError if any bands in `da` are not present in the `band_weights` dictionary. Parameters ---------- da : xarray.DataArray object DataArray containing multispectral data. The dataarray should contain a "variable" dimension that corresponds to the different bands of the data. band_weights : dict Mapping of band names to weights to be applied. The keys of the dictionary should be the names of the bands in the "variable" dimension of `da`, and the values should be the weights to be applied to each band. Returns ------- xarray.DataArray object DataArray with weights applied to the bands. """ # Identify any bands without weights, and raise an # error if they exist bands_without_weights = set(da["variable"].values) - set(band_weights.keys()) if len(bands_without_weights) > 0: raise ValueError( f"The following multispectral bands are missing from the " f"`band_weights` dictionary: {bands_without_weights}.\n" f"Ensure that weights are supplied for all multispectral " f"bands in `ds`, or set `band_weights=None`." ) # Create xr.DataArray with weights for each variable # along the "variable" dimension weights_da = xr.DataArray( data=list(band_weights.values()), coords={"variable": list(band_weights.keys())}, dims="variable", ) # Apply weights return da.weighted(weights_da) def _brovey_pansharpen(ds, pan_band, band_weights=None): """ Perform pansharpening on multiple timesteps of a multispectral dataset using the Brovey transform (with optional per-band weights). Source: https://pro.arcgis.com/en/pro-app/latest/help/analysis/ raster-functions/fundamentals-of-pan-sharpening-pro.htm Parameters ---------- ds : xarray.Dataset Dataset containing multispectral and panchromatic bands. pan_band : str Name of the panchromatic band in the dataset. band_weights : dict, optional Mapping of band names to weights to be applied to each band when calculating the sum of all multispectral bands. The keys of the dictionary should be the names of the bands, and the values should be the weights to apply to each band, e.g.: ``{"nbart_red": 0.4, "nbart_green": 0.4, "nbart_blue": 0.2}``. The default accounts for Landsat 8 and 9's pan band only partially overlapping with the blue band; this may not be suitable for all applications. Setting `band_weights=None` will use a simple unweighted sum. Returns ------- ds_pansharpened : xarray.Dataset Pansharpened dataset with the same dimensions as the input dataset. """ # Create new dataarrays with and without pan band da_nopan = ds.drop(pan_band).to_array() da_pan = ds[pan_band] # Calculate weighted sum if band_weights is not None: da_total = _apply_weights(da_nopan, band_weights).sum(dim="variable") else: da_total = da_nopan.sum(dim="variable") # Perform Brovey Transform in form of: band / total * panchromatic da_pansharpened = da_nopan / da_total * da_pan ds_pansharpened = da_pansharpened.to_dataset("variable") return ds_pansharpened def _esri_pansharpen(ds, pan_band, band_weights=None): """ Perform pansharpening on multiple timesteps of a multispectral dataset using the ESRI transform (with optional per-band weights). Source: https://pro.arcgis.com/en/pro-app/latest/help/analysis/ raster-functions/fundamentals-of-pan-sharpening-pro.htm Parameters ---------- ds : xarray.Dataset Dataset containing multispectral and panchromatic bands. pan_band : str Name of the panchromatic band in the dataset. band_weights : dict, optional Mapping of band names to weights to be applied to each band when calculating the mean of all multispectral bands. The keys of the dictionary should be the names of the bands, and the values should be the weights to apply to each band, e.g.: ``{"nbart_red": 0.4, "nbart_green": 0.4, "nbart_blue": 0.2}``. The default accounts for Landsat 8 and 9's pan band only partially overlapping with the blue band; this may not be suitable for all applications. Setting `band_weights=None` will use a simple unweighted mean. Returns ------- ds_pansharpened : xarray.Dataset Pansharpened dataset with the same dimensions as the input dataset. """ # Create new dataarrays with and without pan band da_nopan = ds.drop(pan_band).to_array() da_pan = ds[pan_band] # Calculate weighted sum if band_weights is not None: da_mean = _apply_weights(da_nopan, band_weights).mean(dim="variable") else: da_mean = da_nopan.mean(dim="variable") # Calculate adjustment and apply to multispectral bands adj = da_pan - da_mean da_pansharpened = da_nopan + adj ds_pansharpened = da_pansharpened.to_dataset("variable") return ds_pansharpened def _simple_mean_pansharpen(ds, pan_band): """ Perform pansharpening on multiple timesteps of a multispectral dataset using the Simple Mean transform. Source: https://pro.arcgis.com/en/pro-app/latest/help/analysis/ raster-functions/fundamentals-of-pan-sharpening-pro.htm Parameters ---------- ds : xarray.Dataset Dataset containing multispectral and panchromatic bands. pan_band : str Name of the panchromatic band in the dataset. Returns ------- ds_pansharpened : xarray.Dataset Pansharpened dataset with the same dimensions as the input dataset. """ # Create new dataarrays with and without pan band ds_nopan = ds.drop(pan_band) da_pan = ds[pan_band] # Take mean of pan band and RGBs ds_pansharpened = (ds_nopan + da_pan) / 2.0 return ds_pansharpened def _hsv_timestep_pansharpen(ds_i, pan_band): """ Perform pansharpening on a single timestep of a multispectral dataset using the Hue Saturation Value (HSV) transform. Parameters ---------- ds : xarray.Dataset Dataset containing multispectral and panchromatic bands. pan_band : str Name of the panchromatic band in the dataset. Returns ------- ds_pansharpened : xarray.Dataset Pansharpened dataset with the same dimensions as the input dataset. """ # Convert to an xr.DataArray and move "variable" to end da_i = ds_i.to_array().transpose(..., "variable") # Create new dataarrays with and without pan band da_i_nopan = da_i.drop(pan_band, dim="variable") da_i_pan = da_i.sel(variable=pan_band) # Convert to HSV colour space hsv = rgb2hsv(da_i_nopan) # Replace value (lightness) channel with pan band data hsv[:, :, 2] = da_i_pan.values # Convert back to RGB colour space pansharped_array = hsv2rgb(hsv) # Add back into original array, reshape and return dataframe da_i_nopan[:] = pansharped_array ds_pansharpened = da_i_nopan.to_dataset("variable") return ds_pansharpened def _pca_timestep_pansharpen(ds_i, pan_band, pca_rescaling="histogram"): """ Perform pansharpening on a single timestep of a multispectral dataset using the principal component analysis (PCA) transform. Parameters ---------- ds : xarray.Dataset Dataset containing multispectral and panchromatic bands. pan_band : str Name of the panchromatic band in the dataset. pca_rescaling : str, optional Method to use for rescaling pan band to more closely match the distribution of values in the first PCA component. "simple" scales the pan band values to more closely match the first PCA component by subtracting the mean of the pan band values from each value, scaling the resulting values by the ratio of the standard deviations of the first PCA component and the pan band, and adding back the mean of the first PCA component. "histogram" uses a histogram matching technique to adjust the pan band values so that the resulting histogram more closely matches the histogram of the first PCA component. Returns ------- ds_pansharpened : xarray.Dataset Pansharpened dataset with the same dimensions as the input dataset. """ # Reshape to 2D by stacking x and y dimensions to prepare it # as an input to PCA. Drop NA rows as these are not supported # by `pca.fit_transform`. da_2d = ( ds_i.to_array() .stack(pixel=("y", "x")) .transpose("pixel", "variable") .dropna(dim="pixel") ) # Create new dataarrays with and without pan band da_2d_nopan = da_2d.drop(pan_band, dim="variable") da_2d_pan = da_2d.sel(variable=pan_band) # Apply PCA transformation pca = sklearn.decomposition.PCA() pca_array = pca.fit_transform(da_2d_nopan) # Rescale pan band to more closely match the first PCA component if pca_rescaling == "simple": pca_array[:, 0] = (da_2d_pan.values - da_2d_pan.values.mean()) * ( pca_array[:, 0].std() / da_2d_pan.values.std() ) + pca_array[:, 0].mean() elif pca_rescaling == "histogram": pca_array[:, 0] = match_histograms(da_2d_pan.values, pca_array[:, 0]) # Apply reverse PCA transform to restore multispectral array pansharped_array = pca.inverse_transform(pca_array) # Add back into original array, reshape and return dataframe da_2d_nopan[:] = pansharped_array ds_pansharpened = da_2d_nopan.unstack("pixel").to_dataset("variable") return ds_pansharpened
[docs] def xr_pansharpen( ds, transform, pan_band="nbart_panchromatic", return_pan=False, output_dtype=None, parallelise=False, band_weights={"nbart_red": 0.4, "nbart_green": 0.4, "nbart_blue": 0.2}, pca_rescaling="histogram", ): """ Apply pan-sharpening to multispectral satellite data with one or more timesteps. The following pansharpening transforms are currently supported: - Brovey ("brovey"), with optional band weighting - ESRI ("esri"), with optional band weighting - Simple mean ("simple mean") - PCA ("pca") - HSV ("hsv"), similar to IHS Note: Pan-sharpening transforms do not necessarily maintain the spectral integrity of the input satellite data, and may be more suitable for visualisation than quantitative work. Parameters ---------- ds : xarray.Dataset An xarrray dataset containing the three input multispectral bands, and a panchromatic band. This dataset should have already been resampled to the spatial resolution of the panchromatic band (15 m for Landsat). Due to differences in the electromagnetic spectrum covered by the panchromatic band, Landsat 8 and 9 data should be supplied with 'blue', 'green', and 'red' multispectral bands, while Landsat 7 should be supplied with 'green', 'red' and 'NIR'. transform : string The pansharpening transform to apply to the data. Valid options include "brovey", "esri", "simple mean", "pca", "hsv". pan_band : string, optional The name of the panchromatic band that will be used to pansharpen the multispectral data. return_pan : bool, optional Whether to return the panchromatic band in the output dataset. Defaults to False. output_dtype : string or numpy.dtype, optional The dtype used for the output values. Defaults to the input dtype of the multispectral bands in `ds`. parallelise: bool, optional Whether to parallelise transformations across multiple cores. Used for PCA and HSV transforms that are applied to each timestep in `ds` individually; defaults to False. band_weights : dict, optional Used for the Brovey and ESRI transforms. Mapping of band names to weights to be applied to each band when calculating the sum (Brovey) or mean (ESRI) of all multispectral bands. The keys of the dictionary should be the names of the bands, and the values should be the weights to apply to each band, e.g.: ``{"nbart_red": 0.4, "nbart_green": 0.4, "nbart_blue": 0.2}``. The default accounts for Landsat 8 and 9's pan band only partially overlapping with the blue band; this may not be suitable for all applications. Setting `band_weights=None` will use a simple unweighted sum (for the Brovey transform) or unweighted mean (for the ESRI transform). pca_rescaling : str, optional Used for the PCA transform. The method to use for rescaling pan band to more closely match the distribution of values in the first PCA component. "simple" scales the pan band values to more closely match the first PCA component by subtracting the mean of the pan band values from each value, scaling the resulting values by the ratio of the standard deviations of the first PCA component and the pan band, and adding back the mean of the first PCA component. "histogram" uses a histogram matching technique to adjust the pan band values so that the resulting histogram more closely matches the histogram of the first PCA component. Returns ------- ds_pansharpened : xarray.Dataset An xarrray dataset containing the three pansharpened input multispectral bands and optionally the panchromatic band (if `return_pan=True`). """ # Assert whether pan band exists in the dataset if pan_band not in ds.data_vars: raise ValueError( f"The specified panchromatic band '{pan_band}' cannot be found in `ds`. " f"Specify a panchromatic band name that exists in the dataset using `pan_band=...`." ) # Assert whether exactly three multispectral bands are included in `ds` n_multi = len(ds.drop(pan_band).data_vars) if n_multi != 3: raise ValueError( f"`ds` should contain exactly three multispectral bands (not " f"including the panchromatic band). However, {n_multi} " f"multispectral bands were found: {list(ds.drop(pan_band).data_vars)}. " ) # Define dict linking functions to each transform transform_dict = { "brovey": _brovey_pansharpen, "esri": _esri_pansharpen, "simple mean": _simple_mean_pansharpen, "pca": _pca_timestep_pansharpen, "hsv": _hsv_timestep_pansharpen, } # If Brovey, ESRI or Simple Mean pansharpening is specified, apply to # entire `xr.Dataset` in one go (with optional weights for Brovey, ESRI) if transform in ("brovey", "esri", "simple mean"): print(f"Applying {transform.capitalize()} pansharpening") extra_params = ( {"band_weights": band_weights} if transform in ("brovey", "esri") else {} ) ds_pansharpened = transform_dict[transform]( ds, pan_band=pan_band, **extra_params, ) # Otherwise, apply PCA or HSV pansharpening to each # timestep in the `xr.Dataset` using `.apply` elif transform in ("pca", "hsv"): extra_params = {"pca_rescaling": pca_rescaling} if transform == "pca" else {} # Apply pansharpening to all timesteps in data in parallel if ("time" in ds.dims) and parallelise: print(f"Applying {transform.upper()} pansharpening in parallel") ds_pansharpened = parallel_apply( ds, "time", transform_dict[transform], pan_band, *extra_params.values(), # TODO: Update once `parallel_apply` supports kwargs ) # Apply pansharpening to all timesteps in data sequentially elif ("time" in ds.dims) and not parallelise: print(f"Applying {transform.upper()} pansharpening") ds_pansharpened = ds.groupby("time").apply( transform_dict[transform], pan_band=pan_band, **extra_params, ) # Otherwise, apply func directly if only one timestep else: print(f"Applying {transform.upper()} pansharpening") ds_pansharpened = transform_dict[transform]( ds, pan_band=pan_band, **extra_params ) else: raise ValueError( f"Unsupported value '{transform}' passed to `method`. Please " f"provide one of {list(transform_dict.keys())}." ) # Optionally insert pan band back into dataset if return_pan: ds_pansharpened[pan_band] = ds[pan_band] # Return data in original or requested dtype return ds_pansharpened.astype( ds.to_array().dtype if output_dtype is None else output_dtype )
[docs] def load_reproject( path, how, resolution="auto", tight=False, resampling="nearest", chunks={"x": 2048, "y": 2048}, bands=None, masked=True, reproject_kwds=None, **kwargs, ): """ Load and reproject part of a raster dataset into a given GeoBox or custom CRS/resolution. Parameters ---------- path : str Path to the raster dataset to be loaded and reprojected. how : GeoBox, str or int How to reproject the raster. Can be a GeoBox or a CRS (e.g. "ESPG:XXXX" string or integer). resolution : str or int, optional The resolution to reproject the raster dataset into if `how` is a CRS, by default "auto". Supports: - "same" use exactly the same resolution as the input raster - "fit" use center pixel to determine required scale change - "auto" uses the same resolution on the output if CRS units are the same between source and destination; otherwise "fit" - Else, a specific resolution in the units of the output crs tight : bool, optional By default output pixel grid is adjusted to align pixel edges to X/Y axis, suppling tight=True produces an unaligned geobox. resampling : str, optional Resampling method to use when reprojecting data, by default "nearest", supports all standard GDAL options ("average", "bilinear", "min", "max", "cubic" etc). chunks : dict, optional The size of the Dask chunks to load the data with, by default {"x": 2048, "y": 2048}. bands : str or list, optional Bands to optionally filter to when loading data. masked : bool, optional Whether to mask the data by its nodata value, by default True. reproject_kwds : dict, optional Additional keyword arguments to pass to the `.odc.reproject()` method, by default None. **kwargs : dict Additional keyword arguments to be passed to the `rioxarray.open_rasterio` function. Returns ------- xarray.Dataset The reprojected raster dataset. """ # Use empty kwds if not provided reproject_kwds = {} if reproject_kwds is None else reproject_kwds # Load data with rasterio da = rioxarray.open_rasterio( filename=path, masked=masked, chunks=chunks, **kwargs, ) # Optionally filter to bands if bands is not None: da = da.sel(band=bands) # Reproject into GeoBox da = da.odc.reproject( how=how, resolution=resolution, tight=tight, resampling=resampling, dst_nodata=np.NaN if masked else None, **reproject_kwds, ) # Squeeze if only one band da = da.squeeze() return da