Merge branch 'feature/spectral_homogenization'

examples/notebooks/GMS beta usecase.ipynb

@@ -2089,7 +2089,7 @@
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 3.0
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",

gms_preprocessing/algorithms/L2B_P.py

@@ -12,9 +12,11 @@ import matplotlib.pyplot as plt
 from sklearn.cluster import KMeans
 from pandas import DataFrame
 from typing import Union  # noqa F401  # flake8 issue
+from tqdm import tqdm
 
 from sklearn.cluster import k_means_  # noqa F401  # flake8 issue
 from geoarray import GeoArray  # noqa F401  # flake8 issue
+from py_tools_ds.numeric.array import get_array_tilebounds
 
 from ..config import GMS_config as CFG
 from ..io.input_reader import SRF  # noqa F401  # flake8 issue
@@ -83,6 +85,10 @@ class SpectralResampler(object):
         if wvl_unit not in ['micrometers', 'nanometers']:
             raise ValueError('Unknown wavelength unit %s.' % wvl_unit)
 
+        # privates
+        self._wvl_1nm = None
+        self._srf_1nm = {}
+
         wvl = np.array(wvl_src, dtype=np.float).flatten()
         if srf_tgt.wvl_unit != 'nanometers':
             srf_tgt.convert_wvl_unit()
@@ -92,6 +98,27 @@ class SpectralResampler(object):
         self.wvl_unit = wvl_unit
         self.logger = logger or logging.getLogger(__name__)
 
+    @property
+    def wvl_1nm(self):
+        # spectral resampling of input image to 1 nm resolution
+        if self._wvl_1nm is None:
+            self._wvl_1nm = np.arange(np.ceil(self.wvl_src_nm.min()), np.floor(self.wvl_src_nm.max()), 1)
+        return self._wvl_1nm
+
+    @property
+    def srf_1nm(self):
+        if not self._srf_1nm:
+            for band in self.srf_tgt.bands:
+                # resample srf to 1 nm
+                self._srf_1nm[band] = \
+                    sp.interpolate.interp1d(self.srf_tgt.srfs_wvl, self.srf_tgt.srfs[band],
+                                            bounds_error=False, fill_value=0, kind='linear')(self.wvl_1nm)
+
+                # validate
+                assert len(self._srf_1nm[band]) == len(self.wvl_1nm)
+
+        return self._srf_1nm
+
     def resample_signature(self, spectrum, scale_factor=10000, v=False):
         # type: (np.ndarray, int, bool) -> np.ndarray
         """Resample the given spectrum according to the spectral response functions of the target instument.
@@ -108,75 +135,82 @@ class SpectralResampler(object):
         assert spectrum.size == self.wvl_src_nm.size
 
         # resample input spectrum and wavelength to 1nm
-        wvl_1nm = np.arange(np.ceil(self.wvl_src_nm.min()), np.floor(self.wvl_src_nm.max()), 1)
         spectrum_1nm = interp1d(self.wvl_src_nm, spectrum,
-                                bounds_error=False, fill_value=0, kind='linear')(wvl_1nm)
+                                bounds_error=False, fill_value=0, kind='linear')(self.wvl_1nm)
 
         if v:
             plt.figure()
-            plt.plot(wvl_1nm, spectrum_1nm/scale_factor, '.')
+            plt.plot(self.wvl_1nm, spectrum_1nm/scale_factor, '.')
 
         spectrum_rsp = []
 
         for band, wvl_center in zip(self.srf_tgt.bands, self.srf_tgt.wvl):
-            # resample srf to 1 nm
-            srf_1nm = interp1d(self.srf_tgt.srfs_wvl, self.srf_tgt.srfs[band],
-                               bounds_error=False, fill_value=0, kind='linear')(wvl_1nm)
-
             # compute the resampled spectral value (np.average computes the weighted mean value)
-            specval_rsp = np.average(spectrum_1nm, weights=srf_1nm)
+            specval_rsp = np.average(spectrum_1nm, weights=self.srf_1nm[band])
 
             if v:
-                plt.plot(wvl_1nm, srf_1nm/max(srf_1nm))
+                plt.plot(self.wvl_1nm, self.srf_1nm/max(self.srf_1nm))
                 plt.plot(wvl_center, specval_rsp/scale_factor, 'x', color='r')
 
             spectrum_rsp.append(specval_rsp)
 
         return np.array(spectrum_rsp)
 
-    def resample_image(self, image_cube):
-        # type: (np.ndarray, int, bool) -> np.ndarray
-        """Resample the given spectral image cube according to the spectral response functions of the target instument.
+    def resample_image(self, image_cube, tiledims=(20, 20)):
+        # type: (Union[GeoArray, np.ndarray], tuple) -> np.ndarray
+        """Resample the given spectral image cube according to the spectral response functions of the target instrument.
 
         :param image_cube:      image (3D array) containing the spectral information in the third dimension
+        :param tiledims:        dimension of tiles to be used during computation
         :return:    resampled spectral image cube
         """
         # input validation
-        if not isinstance(image_cube, np.ndarray):
+        if not isinstance(image_cube, (GeoArray, np.ndarray)):
             raise TypeError(image_cube)
         if not image_cube.ndim == 3:
             raise ValueError("The given image cube must be 3-dimensional. Received a %s-dimensional array."
                              % image_cube.ndim)
         assert image_cube.shape[2] == self.wvl_src_nm.size
 
-        # resample input image and wavelength to 1nm
-        wvl_1nm = np.arange(np.ceil(self.wvl_src_nm.min()), np.floor(self.wvl_src_nm.max()), 1)
-        image_1nm = interp1d(self.wvl_src_nm, image_cube,
-                             axis=2, bounds_error=False, fill_value=0, kind='linear')(wvl_1nm)
+        tilebounds = get_array_tilebounds(array_shape=image_cube.shape, tile_shape=tiledims)
 
         (R, C), B = image_cube.shape[:2], len(self.srf_tgt.bands)
         image_rsp = np.zeros((R, C, B), dtype=image_cube.dtype)
+        for bounds in tilebounds:
+            (rs, re), (cs, ce) = bounds
+            tile = image_cube[rs: re+1, cs: ce+1, :] if image_cube.ndim == 3 else image_cube[rs: re+1, cs: ce+1]
 
-        for band_idx, (band, wvl_center) in enumerate(zip(self.srf_tgt.bands, self.srf_tgt.wvl)):
-            self.logger.info('Applying spectral resampling to band %s...' % band_idx)
+            tile_rsp = self._specresample(tile)
 
-            # resample srf to 1 nm
-            srf_1nm = sp.interpolate.interp1d(self.srf_tgt.srfs_wvl, self.srf_tgt.srfs[band],
-                                              bounds_error=False, fill_value=0, kind='linear')(wvl_1nm)
+            if image_cube.ndim == 3:
+                image_rsp[rs: re + 1, cs: ce + 1, :] = tile_rsp
+            else:
+                image_rsp[rs: re + 1, cs: ce + 1] = tile_rsp
 
+        return image_rsp
+
+    def _specresample(self, tile):
+        # spectral resampling of input image to 1 nm resolution
+        tile_1nm = interp1d(self.wvl_src_nm, tile,
+                            axis=2, bounds_error=False, fill_value=0, kind='linear')(self.wvl_1nm)
+
+        tile_rsp = np.zeros((*tile_1nm.shape[:2], len(self.srf_tgt.bands)), dtype=tile.dtype)
+        for band_idx, band in enumerate(self.srf_tgt.bands):
             # compute the resampled image cube (np.average computes the weighted mean value)
-            image_rsp[:, :, band_idx] = np.average(image_1nm, weights=srf_1nm, axis=2)
+            tile_rsp[:, :, band_idx] = np.average(tile_1nm, weights=self.srf_1nm[band], axis=2)
 
-        return image_rsp
+        return tile_rsp
 
 
 class KMeansRSImage(object):
-    _clusters = None
-    _im_clust = None
-    _spectra = None
-
     def __init__(self, im, n_clusters):
         # type: (GeoArray, int) -> None
+
+        # privates
+        self._clusters = None
+        self._im_clust = None
+        self._spectra = None
+
         self.im = im
         self.n_clusters = n_clusters
 
@@ -334,11 +368,17 @@ class SpecHomo_Classifier(object):
             self.logger.info('Reading the input image %s...' % im_name)
             im_gA = GeoArray(im_ref)
             im_gA.cwl = np.array(im_gA.meta.loc['wavelength'], dtype=np.float).flatten()
+            # im_gA.to_mem()
+            wvl_unit = 'nanometers' if max(im_gA.cwl) > 15 else 'micrometers'
 
             # perform spectral resampling of input image to match spectral properties of target sensor
             self.logger.info('Performing spectral resampling to match spectral properties of target sensor...')
-            SR = SpectralResampler(im_gA.cwl, tgt_srf)
-            im_tgt = GeoArray(SR.resample_image(im_gA[:]), geotransform=im_gA.gt, projection=im_gA.prj)
+            SR = SpectralResampler(im_gA.cwl, tgt_srf, wvl_unit=wvl_unit)
+
+            im_tgt = np.empty((*im_gA.shape[:2], len(tgt_srf.bands)))
+            for ((rS, rE), (cS, cE)), tiledata in tqdm(im_gA.tiles((1000, 1000)), total=900):
+                im_tgt[rS: rE + 1, cS: cE + 1, :] = SR.resample_image(tiledata)
+            im_tgt = GeoArray(im_tgt, im_gA.gt, im_gA.prj)
 
             # compute KMeans clusters for the spectrally resampled image
             self.logger.info('Computing %s KMeans clusters...' % n_clusters)

gms_preprocessing/io/input_reader.py

@@ -220,9 +220,9 @@ class SRF(object):
         if wvl_unit not in ['micrometers', 'nanometers']:
             raise ValueError('Unknown wavelength unit %s.' % wvl_unit)
 
-        self.srfs_wvl = []
-        self.srfs = {}
-        self.srfs_norm01 = {}
+        self.srfs_wvl = []  # wavelength positions with 1 nm precision
+        self.srfs = {}  # srf values with 1 nm precision
+        self.srfs_norm01 = {}  # srf values with 1 nm precision
         self.bands = []
         self.wvl = []
         self.wvl_unit = wvl_unit

gms_preprocessing/misc/helper_functions.py

@@ -134,101 +134,6 @@ def convert_absPathArchive_to_GDALvsiPath(path_archive):
     return os.path.join(gdal_prefix_dict[file_suffix], os.path.basename(path_archive))
 
 
-def get_image_tileborders(target_tileShape, target_tileSize, path_GMS_file=None, shape_fullArr=None):
-    """Calculates row/col bounds for image tiles according to the given parameters.
-
-    :param target_tileShape:    <str> 'cube','row','col','band','block','pixel' or 'custom'
-    :param target_tileSize:     None or <tuple>. Only needed if target_tileShape=='block': (rows,cols)
-    :param path_GMS_file:       <str> path to the *.gms file. Only needed if shape_fullArr is not given
-    :param shape_fullArr:       <tuple> (rows,columns,bands)
-    """
-
-    if shape_fullArr is None:
-        shape_fullArr = json.load(open(path_GMS_file))['shape_fullArr']
-    rows, cols = shape_fullArr[:2]
-
-    if target_tileShape == 'row':
-        return [i for i in range(rows)]
-
-    elif target_tileShape == 'col':
-        return [i for i in range(cols)]
-
-    elif target_tileShape == 'pixel':
-        return [[r, c] for r in range(rows) for c in range(cols)]
-
-    elif target_tileShape == 'block':
-        row_bounds = [0]
-        while row_bounds[-1] + target_tileSize[0] < rows:
-            row_bounds.append(row_bounds[-1] + target_tileSize[0] - 1)
-            row_bounds.append(row_bounds[-2] + target_tileSize[0])
-        else:
-            row_bounds.append(rows - 1)
-
-        col_bounds = [0]
-        while col_bounds[-1] + target_tileSize[1] < cols:
-            col_bounds.append(col_bounds[-1] + target_tileSize[1] - 1)
-            col_bounds.append(col_bounds[-2] + target_tileSize[1])
-        else:
-            col_bounds.append(cols - 1)
-        return [[tuple([row_bounds[r], row_bounds[r + 1]]), tuple([col_bounds[c], col_bounds[c + 1]])]
-                for r in range(0, len(row_bounds), 2) for c in range(0, len(col_bounds), 2)]
-
-
-def cut_GMS_obj_into_blocks(tuple__In_obj__blocksize_RowsCols):
-    # type: (tuple) -> list
-    """Cut a GMS object into tiles with respect to raster attributes as well as scene wide attributes.
-
-    :param tuple__In_obj__blocksize_RowsCols: a tuple with GMS_obj as first and [rows,cols] as second element"""
-
-    In_obj, blocksize_RowsCols = tuple__In_obj__blocksize_RowsCols
-    assert type(blocksize_RowsCols) in [list, tuple] and len(blocksize_RowsCols) == 2, \
-        "The argument 'blocksize_RowsCols' must represent a list of size 2."
-    # tilepos = get_image_tileborders('block', blocksize_RowsCols, shape_fullArr=In_obj.shape_fullArr)
-
-    return list(In_obj.to_tiles(blocksize=blocksize_RowsCols))
-
-    # NOTE:  it's better to call get_subset_GMS_obj (also takes care of tile map infos)
-    # for tp in tilepos:
-    #     (rS,rE),(cS,cE) = tp
-    #     tile = parentObjDict[In_obj.proc_level](*initArgsDict[In_obj.proc_level])
-    #     [setattr(tile, i, getattr(In_obj,i)) for i in In_obj.__dict__ \
-    #             if not callable(getattr(In_obj,i)) and not isinstance(getattr(In_obj,i),(GeoArray, np.ndarray))]
-    #     [setattr(tile, i, getattr(In_obj,i)[rS:rE+1, cS:cE+1]) \
-    #             for i in In_obj.__dict__ \
-    #             if not callable(getattr(In_obj,i)) and isinstance(getattr(In_obj,i),(GeoArray, np.ndarray))]
-    #     tile.arr_shape = 'block'
-    #     tile.arr_pos   = tp
-    #
-    #     yield tile
-
-
-def merge_GMS_tiles_to_GMS_obj(list_GMS_tiles):
-    # type: (list) -> L1A_object
-    """Merge separate GMS objects belonging to the same scene-ID to ONE GMS object
-
-    :param list_GMS_tiles: <list> of GMS objects that have been created by cut_GMS_obj_into_blocks()
-    """
-
-    warnings.warn("'merge_GMS_tiles_to_GMS_obj' is deprecated. Use GMS_object.from_tiles() instead!",
-                  DeprecationWarning)
-
-    if 'IMapUnorderedIterator' in str(type(list_GMS_tiles)):
-        list_GMS_tiles = list(list_GMS_tiles)
-    procLvl = list_GMS_tiles[0].proc_level
-    GMS_obj = parentObjDict[procLvl](*initArgsDict[procLvl])
-    [setattr(GMS_obj, i, getattr(list_GMS_tiles[0], i)) for i in list_GMS_tiles[0].__dict__
-     if not callable(getattr(list_GMS_tiles[0], )) and not isinstance(getattr(list_GMS_tiles[0], i), np.ndarray)]
-    GMS_obj.arr_shape = 'cube'
-    GMS_obj.arr_pos = None
-
-    list_arraynames = [i for i in list_GMS_tiles[0].__dict__ if not callable(getattr(list_GMS_tiles[0], i)) and
-                       isinstance(getattr(list_GMS_tiles[0], i), np.ndarray)]
-
-    if list_arraynames:
-        GMS_obj = _numba_array_merger(GMS_obj, list_arraynames, list_GMS_tiles)
-    return GMS_obj
-
-
 @jit
 def _numba_array_merger(GMS_obj, list_arraynames, list_GMS_tiles):
     """private function, e.g. called by merge_GMS_tiles_to_GMS_obj() in order to fasten array merging"""

gms_preprocessing/model/dataset.py

@@ -12,6 +12,7 @@ from py_tools_ds.geo.projection import WKT2EPSG, get_UTMzone
 from py_tools_ds.geo.coord_calc import calc_FullDataset_corner_positions
 from py_tools_ds.geo.coord_trafo import pixelToLatLon, pixelToMapYX, imXY2mapXY
 from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
+from py_tools_ds.numeric.array import get_array_tilebounds
 
 from ..misc.logging import GMS_logger as DatasetLogger
 from ..model.metadata import METADATA, get_LayerBandsAssignment
@@ -432,8 +433,7 @@ class Dataset(object):
 
     def get_tilepos(self, target_tileshape, target_tilesize):
         self.tile_pos = [[target_tileshape, tb]
-                         for tb in HLP_F.get_image_tileborders(target_tileshape, target_tilesize,
-                                                               shape_fullArr=self.shape_fullArr)]
+                         for tb in get_array_tilebounds(array_shape=self.shape_fullArr, tile_shape=target_tilesize)]
 
     @staticmethod
     def rescale_array(inArray, outScaleFactor, inScaleFactor=1):
@@ -578,6 +578,8 @@ class Dataset(object):
         # type: (tuple) -> self
         """Returns a generator object where items represent tiles of the given block size for the GMS object.
 
+        # NOTE:  it's better to call get_subset_obj (also takes care of tile map infos)
+
         :param blocksize:   target dimensions of the generated block tile (rows, columns)
         :return:            <list> of GMS_object tiles
         """
@@ -585,7 +587,7 @@ class Dataset(object):
         assert type(blocksize) in [list, tuple] and len(blocksize) == 2, \
             "The argument 'blocksize_RowsCols' must represent a tuple of size 2."
 
-        tilepos = HLP_F.get_image_tileborders('block', blocksize, shape_fullArr=self.shape_fullArr)
+        tilepos = get_array_tilebounds(array_shape=self.shape_fullArr, tile_shape=blocksize)
 
         for tp in tilepos:
             (xmin, xmax), (ymin, ymax) = tp  # e.g. [(0, 1999), (0, 999)] at a blocksize of 2000*1000 (rowsxcols)

gms_preprocessing/processing/process_controller.py

@@ -28,6 +28,8 @@ from ..config import set_config, GMS_config
 from .multiproc import MAP
 from ..misc.definition_dicts import proc_chain, db_jobs_statistics_def
 
+from py_tools_ds.numeric.array import get_array_tilebounds
+
 if TYPE_CHECKING:
     from collections import OrderedDict  # noqa F401  # flake8 issue
 
@@ -327,8 +329,8 @@ class process_controller(object):
             else:
                 # define tile positions and size
                 def get_tilepos_list(GMSfile):
-                    return HLP_F.get_image_tileborders('block', blocksize,
-                                                       shape_fullArr=INP_R.GMSfile2dict(GMSfile)['shape_fullArr'])
+                    return get_array_tilebounds(array_shape=INP_R.GMSfile2dict(GMSfile)['shape_fullArr'],
+                                                tile_shape=blocksize)
 
                 # get input parameters for creating GMS objects as blocks
                 work = [[GMSfile, ['block', tp]] for GMSfile in GMSfile_list_prevLvl_inDB