Source code for enpt.model.images.image_baseclasses

# -*- coding: utf-8 -*-

# EnPT, EnMAP Processing Tool - A Python package for pre-processing of EnMAP Level-1B data
#
# Copyright (C) 2018-2024 Karl Segl (GFZ Potsdam, segl@gfz-potsdam.de), Daniel Scheffler
# (GFZ Potsdam, danschef@gfz-potsdam.de), Niklas Bohn (GFZ Potsdam, nbohn@gfz-potsdam.de),
# Stéphane Guillaso (GFZ Potsdam, stephane.guillaso@gfz-potsdam.de)
#
# This software was developed within the context of the EnMAP project supported
# by the DLR Space Administration with funds of the German Federal Ministry of
# Economic Affairs and Energy (on the basis of a decision by the German Bundestag:
# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG.
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version. Please note the following exception: `EnPT` depends on tqdm, which
# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files
# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore".
# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <https://www.gnu.org/licenses/>.

"""EnPT EnMAP object base classes."""

from types import SimpleNamespace
from typing import List, Optional  # noqa: F401
import numpy as np

# noinspection PyPackageRequirements
from skimage import exposure  # contained in package requirements as scikit-image

from geoarray import GeoArray

from ...model.metadata import EnMAP_Metadata_L2A_MapGeo  # noqa: F401  # only used for type hint

__author__ = ['Daniel Scheffler', 'Stéphane Guillaso', 'André Hollstein']


[docs] class _EnMAP_Image(object): """EnPT base class for all kinds of EnMAP images. NOTE: - Basic functionality that all EnMAP image objects have in common is to be implemented here. - All EnMAP image classes should (directly or indirectly) inherit from _EnMAP_image. Attributes: - to be listed here. Check help(_EnMAP_image) in the meanwhile! """ def __init__(self): """Load and hold data needed for processing EnMAP Level-1B data to Level-2A. Intended usage: - only to be used as a base class to be inherited from! - Example: class EnMAP_Image_Subclass(_EnMAP_Image): def __init__(self): super(EnMAP_Image_Subclass, self).__init__() """ # add EnMAP specific attributes self.paths = SimpleNamespace() # private attributes self._data = None self._mask_landwater = None self._mask_clouds = None self._mask_cloudshadow = None self._mask_haze = None self._mask_snow = None self._mask_cirrus = None self._dem = None self._deadpixelmap = None self._polymer_logchl = None self._polymer_logfb = None self._polymer_rgli = None self._polymer_rnir = None self._polymer_bitmask = None self._subset = None # FIXME how is _subset to be set? # defaults self.entity_ID = '' self.basename = '' @property def data(self) -> GeoArray: """Return the actual EnMAP image data. Bundled with all the corresponding metadata. Attributes and functions (most important; for a full list check help(self.data)!): - ALL attributes of numpy.ndarray! - is_inmem(bool): True if the image data are completely loaded into memory; False if GeoArray only holds a link to a file on disk. - arr: np.ndarray holding the pixel values (if is_mem is True) - rows(int) - cols(int) - bands(int) - shape(tuple) - gt(list): GDAL geotransform: contains the geocoding - prj(str): WKT projection string - show(): plot the image - show_map(): plot a map of the image (based on cartopy library) - reproject_to_new_grid() Usage (there will soon be detailed instructions on usage at https://git.gfz-potsdam.de/danschef/geoarray): - Use self.data like a normal numpy.ndarray! - NOTE: Operators like *, /, + , - will soon be implemented. In the meanwhile use: result = self.data[:] *10 - How to set self.data? - Link an image file to self.data -> all raster data is read into memory ON DEMAND: self.data = '/path/to/image.tif' # sets self.data to GeoArray('/path/to/image.tif') - Link a numpy.ndarray instance with self.data (remaining attributes like geocoding, projection, etc. are copied from the previous self.data attribute). self.data = numpy.array([[1,2,3],[4,5,6]]) - Set self.data to an existing instance of GeoArray (allows to set specific attributes of GeoArray by yourself) self.data = GeoArray('/path/to/image.tif', geotransform=[...], projection='WKTString') - Delete self.data: del self.data or 'self.data = None' :return: instance of geoarray.GeoArray """ if self._subset is None: return self._data return GeoArray(self._data[self._subset[0]:self._subset[1], self._subset[2]:self._subset[3], :], geotransform=self._data.gt, projection=self._data.prj) @data.setter def data(self, *geoArr_initArgs): if geoArr_initArgs[0] is not None: # TODO this must be able to handle subset inputs in tiled processing if isinstance(geoArr_initArgs[0], np.ndarray): self._data = GeoArray(geoArr_initArgs[0], geotransform=self._data.gt, projection=self._data.prj) else: self._data = GeoArray(*geoArr_initArgs) else: del self.data @data.deleter def data(self): self._data = None @property def mask_landwater(self) -> GeoArray: """Return the land/water mask. pixel values: - 0: background within scene dimensions, e.g. due to missing values/errors - 1: no water - 2: water - 3: background outside the scene dimensions (artifact from resampling between map and sensor geometry) :return: geoarray.GeoArray """ return self._mask_landwater @mask_landwater.setter def mask_landwater(self, *geoArr_initArgs): self._mask_landwater = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_landwater', nodataVal=3) @mask_landwater.deleter def mask_landwater(self): self._mask_landwater = None @property def mask_clouds(self) -> GeoArray: """Return the cloud mask. :return: geoarray.GeoArray """ return self._mask_clouds @mask_clouds.setter def mask_clouds(self, *geoArr_initArgs): self._mask_clouds = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_clouds', nodataVal=0) @mask_clouds.deleter def mask_clouds(self): self._mask_clouds = None @property def mask_cloudshadow(self) -> GeoArray: """Return the cloud shadow mask (0=no cloud shadow, 1=cloud shadow). :return: geoarray.GeoArray """ return self._mask_cloudshadow @mask_cloudshadow.setter def mask_cloudshadow(self, *geoArr_initArgs): self._mask_cloudshadow = \ self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cloudshadow', nodataVal=0) @mask_cloudshadow.deleter def mask_cloudshadow(self): self._mask_cloudshadow = None @property def mask_haze(self) -> GeoArray: """Return the haze mask (0=no haze, 1=haze).. :return: geoarray.GeoArray """ return self._mask_haze @mask_haze.setter def mask_haze(self, *geoArr_initArgs): self._mask_haze = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_haze', nodataVal=0) @mask_haze.deleter def mask_haze(self): self._mask_haze = None @property def mask_snow(self) -> GeoArray: """Return the snow mask (0=no snow, 1=snow). :return: geoarray.GeoArray """ return self._mask_snow @mask_snow.setter def mask_snow(self, *geoArr_initArgs): self._mask_snow = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_snow', nodataVal=0) @mask_snow.deleter def mask_snow(self): self._mask_snow = None @property def mask_cirrus(self) -> GeoArray: """Return the cirrus mask (0=none, 1=thin, 2=medium, 3=thick). :return: geoarray.GeoArray """ return self._mask_cirrus @mask_cirrus.setter def mask_cirrus(self, *geoArr_initArgs): self._mask_cirrus = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cirrus', nodataVal=0) @mask_cirrus.deleter def mask_cirrus(self): self._mask_cirrus = None @property def dem(self) -> GeoArray: """Return a DEM in the exact dimension and pixel grid of self.data. :return: geoarray.GeoArray """ if self._dem is None: raise NotImplementedError('DEM is missing. An automatic DEM getter is currently not implemented.') return self._dem @dem.setter def dem(self, *geoArr_initArgs): self._dem = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'dem', nodataVal=0) # FIXME 0? @dem.deleter def dem(self): self._dem = None @property def deadpixelmap(self) -> GeoArray: """Return the dead pixel map. :return: geoarray.GeoArray """ if self._deadpixelmap is not None: self._deadpixelmap.arr = self._deadpixelmap[:].astype(bool) # ensure boolean map return self._deadpixelmap @deadpixelmap.setter def deadpixelmap(self, *geoArr_initArgs): if geoArr_initArgs[0] is not None: dpm = GeoArray(*geoArr_initArgs) if dpm.ndim == 3 and dpm.shape != self.data.shape: raise ValueError("The 'deadpixelmap' GeoArray can only be instanced with a 3D array with the same size " "like %s.data, i.e.: %s Received %s." % (self.__class__.__name__, str(self.data.shape), str(dpm.shape))) elif dpm.ndim == 2 and dpm.shape != (self.data.bands, self.data.cols): raise ValueError("The 'deadpixelmap' GeoArray can only be instanced with an array with the size " "'bands x columns' of the GeoArray %s.data. Received %s. Expected %s" % (self.__class__.__name__, str(dpm.shape), str((self.data.bands, self.data.cols)))) self._deadpixelmap = dpm else: del self.deadpixelmap @deadpixelmap.deleter def deadpixelmap(self): self._deadpixelmap = None @property def polymer_logchl(self) -> GeoArray: """Chlorophyll concentration, in mg/m3, in 10-based logarithm (generated by Polymer AC for water areas). Data range: -3 to 2.3 :return: geoarray.GeoArray """ return self._polymer_logchl @polymer_logchl.setter def polymer_logchl(self, *geoArr_initArgs): self._polymer_logchl = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_logchl', nodataVal=-9999) @polymer_logchl.deleter def polymer_logchl(self): self._polymer_logchl = None @property def polymer_logfb(self) -> GeoArray: """Particle scattering factor fb in Park & Ruddick (2005) (in 10-based logarithm; generated by Polymer AC). Data range: positive and negative values, maximum range not defined :return: geoarray.GeoArray """ return self._polymer_logfb @polymer_logfb.setter def polymer_logfb(self, *geoArr_initArgs): self._polymer_logfb = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_logfb', nodataVal=-9999) @polymer_logfb.deleter def polymer_logfb(self): self._polymer_logfb = None @property def polymer_rgli(self) -> GeoArray: """Reflectance of the sun glint predicted from ECMWF wind speed (generated by Polymer AC for water areas). Typical data range: 0-0.2 :return: geoarray.GeoArray """ return self._polymer_rgli @polymer_rgli.setter def polymer_rgli(self, *geoArr_initArgs): self._polymer_rgli = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_rgli', nodataVal=-9999) @polymer_rgli.deleter def polymer_rgli(self): self._polymer_rgli = None @property def polymer_rnir(self) -> GeoArray: """TOA reflectance at 865 nm corrected for Rayleigh scattering (generated by Polymer AC for water areas). Typical data range: 0-0.2 :return: geoarray.GeoArray """ return self._polymer_rnir @polymer_rnir.setter def polymer_rnir(self, *geoArr_initArgs): self._polymer_rnir = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_rnir', nodataVal=-9999) @polymer_rnir.deleter def polymer_rnir(self): self._polymer_rnir = None @property def polymer_bitmask(self) -> GeoArray: """Quality flags as generated by Polymer AC (bit-encoded mask). Value 0 represents water, all fine, no flags. Other pixel values are explained below: +--------------------+-------------+--------------------------------------------+ | Flag name | Flag value | Description | +====================+=============+============================================+ | LAND | 1 | Land mask | +--------------------+-------------+--------------------------------------------+ | CLOUD_BASE | 2 | Polymer's basic cloud mask | +--------------------+-------------+--------------------------------------------+ | L1_INVALID | 4 | Invalid level1 pixel | +--------------------+-------------+--------------------------------------------+ | NEGATIVE_BB | 8 | (deprecated flag) | +--------------------+-------------+--------------------------------------------+ | OUT_OF_BOUNDS | 16 | Retrieved marine parameters are outside | | | | valid bounds | +--------------------+-------------+--------------------------------------------+ | EXCEPTION | 32 | A processing error was encountered | +--------------------+-------------+--------------------------------------------+ | THICK_AEROSOL | 64 | Thick aerosol flag | +--------------------+-------------+--------------------------------------------+ | HIGH_AIR_MASS | 128 | Air mass exceeds 5 | +--------------------+-------------+--------------------------------------------+ | EXTERNAL_MASK | 512 | Pixel was masked using external mask | +--------------------+-------------+--------------------------------------------+ | CASE2 | 1024 | Pixel was processed in "case2" mode | +--------------------+-------------+--------------------------------------------+ | INCONSISTENCY | 2048 | Inconsistent result was detected | | | | (atmospheric reflectance out of bounds) | +--------------------+-------------+--------------------------------------------+ | ANOMALY_RWMOD_BLUE | 4096 | Excessive difference was found at 412nm | | | | between Rw and Rwmod | +--------------------+-------------+--------------------------------------------+ :return: geoarray.GeoArray """ return self._polymer_bitmask @polymer_bitmask.setter def polymer_bitmask(self, *geoArr_initArgs): self._polymer_bitmask = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_bitmask', nodataVal=-9999) @polymer_bitmask.deleter def polymer_bitmask(self): self._polymer_bitmask = None
[docs] def _get_geoarray_with_datalike_geometry(self, geoArr_initArgs: tuple, attrName: str, nodataVal: int = None, specialclass=None) -> Optional[GeoArray]: if geoArr_initArgs[0] is None: return None else: GeoArrayOrSubclass = GeoArray if not specialclass else specialclass gA = GeoArrayOrSubclass(*geoArr_initArgs) if gA.shape[:2] != self.data.shape[:2]: raise ValueError("The '%s' GeoArray can only be instanced with an array with " "the same X/Y dimensions like %s.data %s. Got %s." % (attrName, self.__class__.__name__, str(self.data.shape[:2]), str(gA.shape[:2]))) # noinspection PyProtectedMember if gA._nodata is None and nodataVal is not None: gA.nodata = nodataVal gA.gt = self.data.gt gA.prj = self.data.prj return gA
[docs] def generate_quicklook(self, bands2use: List[int]) -> GeoArray: """ Generate image quicklook and save into a file. :param bands2use: band indices of self.data to be used as (red, green, blue) for quicklook image, e.g., [3, 2, 1] :return: GeoArray """ red, green, blue = [self.data[:, :, bandidx] for bandidx in bands2use] def rescale(bandarray): pixvals = np.ma.masked_equal(bandarray, self.data.nodata).compressed() \ if self.data.nodata is not None else bandarray p2, p98 = np.nanpercentile(pixvals, 2), np.nanpercentile(pixvals, 98) return exposure.rescale_intensity(bandarray, in_range=(p2, p98), out_range=(0, 255)) pix = np.dstack((rescale(red), rescale(green), rescale(blue))).astype(np.uint8) return GeoArray(pix)