|
| 1 | +#!/usr/bin/env python |
| 2 | +# -*- coding: utf-8 -*- |
| 3 | +# Copyright (c) 2023-2024 HyperGas developers |
| 4 | +# |
| 5 | +# This file is part of hypergas. |
| 6 | +# |
| 7 | +# hypergas is a library to retrieve trace gases from hyperspectral satellite data |
| 8 | +"""Create a 2D quality mask for hyperspectral satellite data.""" |
| 9 | + |
| 10 | +import logging |
| 11 | +import numpy as np |
| 12 | +import xarray as xr |
| 13 | + |
| 14 | +from .unit_spectrum import Unit_spec |
| 15 | + |
| 16 | +LOG = logging.getLogger(__name__) |
| 17 | + |
| 18 | + |
| 19 | +class QualityMask(): |
| 20 | + """Quality mask for hyperspectral satellite scenes. |
| 21 | +
|
| 22 | + Computes per-pixel boolean flags for water, cloud, and cirrus |
| 23 | + contamination based on top-of-atmosphere (TOA) reflectance thresholds, |
| 24 | + and combines them into a single ``qmask`` DataArray. |
| 25 | +
|
| 26 | + Parameters |
| 27 | + ---------- |
| 28 | + scn : :class:`~satpy.Scene` |
| 29 | + Satpy Scene that must contain the following datasets: |
| 30 | +
|
| 31 | + * ``'radiance'`` – spectral radiance in W m-2 sr-1 um-1, |
| 32 | + with dimensions ``(bands, y, x)``. |
| 33 | + * ``'sza'`` – solar zenith angle in degrees, shape ``(y, x)``. |
| 34 | +
|
| 35 | + Attributes |
| 36 | + ---------- |
| 37 | + scn : :class:`~satpy.Scene` |
| 38 | + The input scene (stored for reference). |
| 39 | + rad : :class:`~xarray.DataArray` |
| 40 | + Radiance converted to uW cm-2 sr-1 nm-1 (divided by 10). |
| 41 | + sza : float |
| 42 | + Scene-mean solar zenith angle in radians. |
| 43 | + rho : :class:`~xarray.DataArray` |
| 44 | + TOA apparent reflectance, shape ``(bands, y, x)``. |
| 45 | + water : :class:`~xarray.DataArray` |
| 46 | + Boolean water mask, shape ``(y, x)``. |
| 47 | + cloud : :class:`~xarray.DataArray` |
| 48 | + Boolean cloud mask, shape ``(y, x)``. |
| 49 | + cirrus : :class:`~xarray.DataArray` |
| 50 | + Boolean cirrus mask, shape ``(y, x)``. |
| 51 | + qmask : :class:`~xarray.DataArray` |
| 52 | + Combined quality mask with a ``quality_flag`` dimension whose |
| 53 | + labels are ``['water', 'cloud', 'cirrus', 'invalid']``. |
| 54 | + """ |
| 55 | + |
| 56 | + def __init__(self, scn): |
| 57 | + # Store scene and convert radiance units: |
| 58 | + # W m-2 sr-1 um-1 -> uW cm-2 sr-1 nm-1 (factor = 1/10) |
| 59 | + self.scn = scn |
| 60 | + self.rad = scn['radiance'] / 10 |
| 61 | + |
| 62 | + # Scene-mean solar zenith angle (degrees -> radians) |
| 63 | + self.sza = np.deg2rad(float(scn['sza'].mean())) |
| 64 | + |
| 65 | + # Pre-compute TOA reflectance used by all mask methods |
| 66 | + self._toa() |
| 67 | + |
| 68 | + def _toa(self): |
| 69 | + """Compute top-of-atmosphere (apparent) reflectance. |
| 70 | +
|
| 71 | + The TOA reflectance is defined as: |
| 72 | +
|
| 73 | + .. math:: |
| 74 | +
|
| 75 | + \\rho(\\lambda) = |
| 76 | + \\frac{\\pi \\cdot L(\\lambda)}{E_0(\\lambda) \\cdot \\cos(\\theta_s)} |
| 77 | +
|
| 78 | + where :math:`L` is the at-sensor radiance (uW cm-2 sr-1 nm-1), |
| 79 | + :math:`E_0` is the extraterrestrial solar irradiance convolved to |
| 80 | + the sensor's spectral response function (uW cm-2 nm-1), and |
| 81 | + :math:`\\theta_s` is the solar zenith angle. |
| 82 | +
|
| 83 | + The result is stored in ``self.rho`` with shape ``(bands, y, x)``. |
| 84 | + """ |
| 85 | + # Build a Unit_spec object to access the solar irradiance spectrum |
| 86 | + unit = Unit_spec( |
| 87 | + self.scn['radiance'], |
| 88 | + self.scn['radiance'].coords['bands'], |
| 89 | + self.scn['radiance'].coords['bands'].min(), |
| 90 | + self.scn['radiance'].coords['bands'].max(), |
| 91 | + ) |
| 92 | + irr = unit.solar_irradiance |
| 93 | + |
| 94 | + # convert to uW cm-2 nm-1 |
| 95 | + irr /= 10 |
| 96 | + irr.attrs['units'] = 'uW cm-2 nm-1' |
| 97 | + |
| 98 | + # Convolve the high-resolution solar spectrum with the sensor SRF |
| 99 | + irr_resampled = unit._convolve( |
| 100 | + unit.wvl_sensor, |
| 101 | + unit.fwhm_sensor, |
| 102 | + irr.coords['wavelength'].values, |
| 103 | + irr.values, |
| 104 | + ) |
| 105 | + |
| 106 | + # Compute TOA reflectance: ρ = (π · L) / (E₀ · cos θ_s) |
| 107 | + # Transpose operations keep xarray dimension alignment correct |
| 108 | + # Keep rho as a lazy dask array — don't load the full cube here |
| 109 | + irr_resampled = xr.DataArray( |
| 110 | + irr_resampled, |
| 111 | + dims=["bands"], |
| 112 | + coords={"bands": self.rad["bands"]} |
| 113 | + ) |
| 114 | + |
| 115 | + self.rho = (np.pi * self.rad) / (irr_resampled * np.cos(self.sza)) |
| 116 | + |
| 117 | + def mask(self): |
| 118 | + """Compute all individual masks and combine into ``self.qmask``. |
| 119 | +
|
| 120 | + Calls :meth:`water_mask`, :meth:`cloud_mask`, and |
| 121 | + :meth:`cirrus_mask` in sequence, then concatenates the results |
| 122 | + along a new ``quality_flag`` dimension. An additional |
| 123 | + ``'invalid'`` flag is appended that is ``True`` wherever *any* |
| 124 | + of the three masks is ``True``. |
| 125 | +
|
| 126 | + After calling this method the combined mask is available as |
| 127 | + ``self.qmask`` with shape ``(quality_flag, y, x)`` and |
| 128 | + ``quality_flag`` labels ``['water', 'cloud', 'cirrus', 'invalid']``. |
| 129 | + """ |
| 130 | + LOG.info('Generating quality masks using TOA ...') |
| 131 | + |
| 132 | + # Load all required bands in ONE dask compute instead of 5 separate ones |
| 133 | + bands_needed = [450, 1000, 1250, 1380, 1650] |
| 134 | + rho_subset = self.rho.sel(bands=bands_needed, method='nearest').load() |
| 135 | + |
| 136 | + self.water_mask(rho_subset) |
| 137 | + self.cloud_mask(rho_subset) |
| 138 | + self.cirrus_mask(rho_subset) |
| 139 | + |
| 140 | + # Drop any scalar/spectral coordinates inherited from .sel(bands=...) |
| 141 | + # (e.g. 'bands', 'fwhm', 'wavelength') so that xr.concat finds a |
| 142 | + # consistent set of coordinates across all three masks. |
| 143 | + spatial_masks = [ |
| 144 | + m.drop_vars([c for c in m.coords if c not in m.dims]) |
| 145 | + for m in (self.water, self.cloud, self.cirrus) |
| 146 | + ] |
| 147 | + |
| 148 | + # Stack the three boolean masks along a new 'quality_flag' dimension |
| 149 | + qmask = xr.concat( |
| 150 | + spatial_masks, |
| 151 | + dim=xr.DataArray( |
| 152 | + ['water', 'cloud', 'cirrus'], |
| 153 | + dims='quality_flag', |
| 154 | + name='quality_flag', |
| 155 | + ), |
| 156 | + ) |
| 157 | + |
| 158 | + # Derive a single 'invalid' flag: True if any individual flag is True |
| 159 | + any_mask = qmask.any(dim='quality_flag') |
| 160 | + any_mask = any_mask.expand_dims(quality_flag=['invalid']) |
| 161 | + |
| 162 | + # Append 'invalid' to produce the final 4-flag mask |
| 163 | + qmask = xr.concat([qmask, any_mask], dim='quality_flag').astype(float) |
| 164 | + self.qmask = qmask.rename('quality_mask') |
| 165 | + |
| 166 | + def water_mask(self, rho=None): |
| 167 | + """Identify water pixels using TOA reflectance at 1000 nm. |
| 168 | +
|
| 169 | + A pixel is flagged as water when its near-infrared reflectance |
| 170 | + falls below 0.05, exploiting the strong absorption of liquid |
| 171 | + water beyond 900 nm. |
| 172 | +
|
| 173 | + Sets ``self.water`` to a boolean :class:`~xarray.DataArray` |
| 174 | + of shape ``(y, x)``. |
| 175 | + """ |
| 176 | + rho = rho if rho is not None else self.rho.load() |
| 177 | + rho_1000 = rho.sel(bands=1000, method='nearest') |
| 178 | + self.water = rho_1000 < 0.05 |
| 179 | + |
| 180 | + def cloud_mask(self, rho=None): |
| 181 | + """Identify cloud pixels using multi-band TOA reflectance thresholds. |
| 182 | +
|
| 183 | + Three independent reflectance tests are applied: |
| 184 | +
|
| 185 | + * 450 nm > 0.28 (high visible reflectance) |
| 186 | + * 1250 nm > 0.46 (high short-wave infrared reflectance) |
| 187 | + * 1650 nm > 0.22 (high short-wave infrared reflectance) |
| 188 | +
|
| 189 | + A pixel is flagged as cloudy only when **all three** conditions |
| 190 | + are satisfied (majority vote ≥ 3 out of 3), reducing false |
| 191 | + positives over bright land surfaces. |
| 192 | +
|
| 193 | + Sets ``self.cloud`` to a boolean :class:`~xarray.DataArray` |
| 194 | + of shape ``(y, x)``. |
| 195 | +
|
| 196 | + References |
| 197 | + ---------- |
| 198 | + Sandford et al., *AMT*, 13, 7047–7057, 2020. |
| 199 | + https://doi.org/10.5194/amt-13-7047-2020 |
| 200 | + """ |
| 201 | + rho = rho if rho is not None else self.rho.load() |
| 202 | + rho_450 = rho.sel(bands=450, method='nearest') |
| 203 | + rho_1250 = rho.sel(bands=1250, method='nearest') |
| 204 | + rho_1650 = rho.sel(bands=1650, method='nearest') |
| 205 | + |
| 206 | + # Cast each threshold test to int so they can be summed |
| 207 | + self.cloud = ( |
| 208 | + (rho_450 > 0.28).astype(int) |
| 209 | + + (rho_1250 > 0.46).astype(int) |
| 210 | + + (rho_1650 > 0.22).astype(int) |
| 211 | + ) >= 3 |
| 212 | + |
| 213 | + def cirrus_mask(self, rho=None): |
| 214 | + """Identify cirrus cloud pixels using TOA reflectance at 1380 nm. |
| 215 | +
|
| 216 | + The 1380 nm water-vapour absorption band is used as a cirrus |
| 217 | + proxy: surface-leaving radiance is almost entirely absorbed by |
| 218 | + atmospheric water vapour at this wavelength, so any residual |
| 219 | + reflectance above the threshold is attributed to high-altitude |
| 220 | + cirrus ice clouds. |
| 221 | +
|
| 222 | + Sets ``self.cirrus`` to a boolean :class:`~xarray.DataArray` |
| 223 | + of shape ``(y, x)``. |
| 224 | +
|
| 225 | + References |
| 226 | + ---------- |
| 227 | + Gao & Goetz, *GRL*, 20(4), 301–304, 1993. |
| 228 | + https://doi.org/10.1029/93GL00106 |
| 229 | + """ |
| 230 | + rho = rho if rho is not None else self.rho.load() |
| 231 | + rho_1380 = rho.sel(bands=1380, method='nearest') |
| 232 | + self.cirrus = rho_1380 > 0.1 |
0 commit comments