Commit 49941090 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'feature/spectral_homogenization'


Former-commit-id: f75e9ef7
Former-commit-id: 42822934
parents 711986fd 2166856c
...@@ -2089,7 +2089,7 @@ ...@@ -2089,7 +2089,7 @@
"language_info": { "language_info": {
"codemirror_mode": { "codemirror_mode": {
"name": "ipython", "name": "ipython",
"version": 3.0 "version": 3
}, },
"file_extension": ".py", "file_extension": ".py",
"mimetype": "text/x-python", "mimetype": "text/x-python",
......
...@@ -12,9 +12,11 @@ import matplotlib.pyplot as plt ...@@ -12,9 +12,11 @@ import matplotlib.pyplot as plt
from sklearn.cluster import KMeans from sklearn.cluster import KMeans
from pandas import DataFrame from pandas import DataFrame
from typing import Union # noqa F401 # flake8 issue from typing import Union # noqa F401 # flake8 issue
from tqdm import tqdm
from sklearn.cluster import k_means_ # noqa F401 # flake8 issue from sklearn.cluster import k_means_ # noqa F401 # flake8 issue
from geoarray import GeoArray # 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 ..config import GMS_config as CFG
from ..io.input_reader import SRF # noqa F401 # flake8 issue from ..io.input_reader import SRF # noqa F401 # flake8 issue
...@@ -83,6 +85,10 @@ class SpectralResampler(object): ...@@ -83,6 +85,10 @@ class SpectralResampler(object):
if wvl_unit not in ['micrometers', 'nanometers']: if wvl_unit not in ['micrometers', 'nanometers']:
raise ValueError('Unknown wavelength unit %s.' % wvl_unit) 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() wvl = np.array(wvl_src, dtype=np.float).flatten()
if srf_tgt.wvl_unit != 'nanometers': if srf_tgt.wvl_unit != 'nanometers':
srf_tgt.convert_wvl_unit() srf_tgt.convert_wvl_unit()
...@@ -92,6 +98,27 @@ class SpectralResampler(object): ...@@ -92,6 +98,27 @@ class SpectralResampler(object):
self.wvl_unit = wvl_unit self.wvl_unit = wvl_unit
self.logger = logger or logging.getLogger(__name__) 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): def resample_signature(self, spectrum, scale_factor=10000, v=False):
# type: (np.ndarray, int, bool) -> np.ndarray # type: (np.ndarray, int, bool) -> np.ndarray
"""Resample the given spectrum according to the spectral response functions of the target instument. """Resample the given spectrum according to the spectral response functions of the target instument.
...@@ -108,75 +135,82 @@ class SpectralResampler(object): ...@@ -108,75 +135,82 @@ class SpectralResampler(object):
assert spectrum.size == self.wvl_src_nm.size assert spectrum.size == self.wvl_src_nm.size
# resample input spectrum and wavelength to 1nm # 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, 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: if v:
plt.figure() plt.figure()
plt.plot(wvl_1nm, spectrum_1nm/scale_factor, '.') plt.plot(self.wvl_1nm, spectrum_1nm/scale_factor, '.')
spectrum_rsp = [] spectrum_rsp = []
for band, wvl_center in zip(self.srf_tgt.bands, self.srf_tgt.wvl): 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) # 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: 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') plt.plot(wvl_center, specval_rsp/scale_factor, 'x', color='r')
spectrum_rsp.append(specval_rsp) spectrum_rsp.append(specval_rsp)
return np.array(spectrum_rsp) return np.array(spectrum_rsp)
def resample_image(self, image_cube): def resample_image(self, image_cube, tiledims=(20, 20)):
# type: (np.ndarray, int, bool) -> np.ndarray # type: (Union[GeoArray, np.ndarray], tuple) -> np.ndarray
"""Resample the given spectral image cube according to the spectral response functions of the target instument. """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 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 :return: resampled spectral image cube
""" """
# input validation # input validation
if not isinstance(image_cube, np.ndarray): if not isinstance(image_cube, (GeoArray, np.ndarray)):
raise TypeError(image_cube) raise TypeError(image_cube)
if not image_cube.ndim == 3: if not image_cube.ndim == 3:
raise ValueError("The given image cube must be 3-dimensional. Received a %s-dimensional array." raise ValueError("The given image cube must be 3-dimensional. Received a %s-dimensional array."
% image_cube.ndim) % image_cube.ndim)
assert image_cube.shape[2] == self.wvl_src_nm.size assert image_cube.shape[2] == self.wvl_src_nm.size
# resample input image and wavelength to 1nm tilebounds = get_array_tilebounds(array_shape=image_cube.shape, tile_shape=tiledims)
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)
(R, C), B = image_cube.shape[:2], len(self.srf_tgt.bands) (R, C), B = image_cube.shape[:2], len(self.srf_tgt.bands)
image_rsp = np.zeros((R, C, B), dtype=image_cube.dtype) 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)): tile_rsp = self._specresample(tile)
self.logger.info('Applying spectral resampling to band %s...' % band_idx)
# resample srf to 1 nm if image_cube.ndim == 3:
srf_1nm = sp.interpolate.interp1d(self.srf_tgt.srfs_wvl, self.srf_tgt.srfs[band], image_rsp[rs: re + 1, cs: ce + 1, :] = tile_rsp
bounds_error=False, fill_value=0, kind='linear')(wvl_1nm) 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) # 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): class KMeansRSImage(object):
_clusters = None
_im_clust = None
_spectra = None
def __init__(self, im, n_clusters): def __init__(self, im, n_clusters):
# type: (GeoArray, int) -> None # type: (GeoArray, int) -> None
# privates
self._clusters = None
self._im_clust = None
self._spectra = None
self.im = im self.im = im
self.n_clusters = n_clusters self.n_clusters = n_clusters
...@@ -334,11 +368,17 @@ class SpecHomo_Classifier(object): ...@@ -334,11 +368,17 @@ class SpecHomo_Classifier(object):
self.logger.info('Reading the input image %s...' % im_name) self.logger.info('Reading the input image %s...' % im_name)
im_gA = GeoArray(im_ref) im_gA = GeoArray(im_ref)
im_gA.cwl = np.array(im_gA.meta.loc['wavelength'], dtype=np.float).flatten() 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 # 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...') self.logger.info('Performing spectral resampling to match spectral properties of target sensor...')
SR = SpectralResampler(im_gA.cwl, tgt_srf) SR = SpectralResampler(im_gA.cwl, tgt_srf, wvl_unit=wvl_unit)
im_tgt = GeoArray(SR.resample_image(im_gA[:]), geotransform=im_gA.gt, projection=im_gA.prj)
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 # compute KMeans clusters for the spectrally resampled image
self.logger.info('Computing %s KMeans clusters...' % n_clusters) self.logger.info('Computing %s KMeans clusters...' % n_clusters)
......
...@@ -220,9 +220,9 @@ class SRF(object): ...@@ -220,9 +220,9 @@ class SRF(object):
if wvl_unit not in ['micrometers', 'nanometers']: if wvl_unit not in ['micrometers', 'nanometers']:
raise ValueError('Unknown wavelength unit %s.' % wvl_unit) raise ValueError('Unknown wavelength unit %s.' % wvl_unit)
self.srfs_wvl = [] self.srfs_wvl = [] # wavelength positions with 1 nm precision
self.srfs = {} self.srfs = {} # srf values with 1 nm precision
self.srfs_norm01 = {} self.srfs_norm01 = {} # srf values with 1 nm precision
self.bands = [] self.bands = []
self.wvl = [] self.wvl = []
self.wvl_unit = wvl_unit self.wvl_unit = wvl_unit
......
...@@ -134,101 +134,6 @@ def convert_absPathArchive_to_GDALvsiPath(path_archive): ...@@ -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)) 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 @jit
def _numba_array_merger(GMS_obj, list_arraynames, list_GMS_tiles): 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""" """private function, e.g. called by merge_GMS_tiles_to_GMS_obj() in order to fasten array merging"""
......
...@@ -12,6 +12,7 @@ from py_tools_ds.geo.projection import WKT2EPSG, get_UTMzone ...@@ -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_calc import calc_FullDataset_corner_positions
from py_tools_ds.geo.coord_trafo import pixelToLatLon, pixelToMapYX, imXY2mapXY 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.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 ..misc.logging import GMS_logger as DatasetLogger
from ..model.metadata import METADATA, get_LayerBandsAssignment from ..model.metadata import METADATA, get_LayerBandsAssignment
...@@ -432,8 +433,7 @@ class Dataset(object): ...@@ -432,8 +433,7 @@ class Dataset(object):
def get_tilepos(self, target_tileshape, target_tilesize): def get_tilepos(self, target_tileshape, target_tilesize):
self.tile_pos = [[target_tileshape, tb] self.tile_pos = [[target_tileshape, tb]
for tb in HLP_F.get_image_tileborders(target_tileshape, target_tilesize, for tb in get_array_tilebounds(array_shape=self.shape_fullArr, tile_shape=target_tilesize)]
shape_fullArr=self.shape_fullArr)]
@staticmethod @staticmethod
def rescale_array(inArray, outScaleFactor, inScaleFactor=1): def rescale_array(inArray, outScaleFactor, inScaleFactor=1):
...@@ -578,6 +578,8 @@ class Dataset(object): ...@@ -578,6 +578,8 @@ class Dataset(object):
# type: (tuple) -> self # type: (tuple) -> self
"""Returns a generator object where items represent tiles of the given block size for the GMS object. """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) :param blocksize: target dimensions of the generated block tile (rows, columns)
:return: <list> of GMS_object tiles :return: <list> of GMS_object tiles
""" """
...@@ -585,7 +587,7 @@ class Dataset(object): ...@@ -585,7 +587,7 @@ class Dataset(object):
assert type(blocksize) in [list, tuple] and len(blocksize) == 2, \ assert type(blocksize) in [list, tuple] and len(blocksize) == 2, \
"The argument 'blocksize_RowsCols' must represent a tuple of size 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: for tp in tilepos:
(xmin, xmax), (ymin, ymax) = tp # e.g. [(0, 1999), (0, 999)] at a blocksize of 2000*1000 (rowsxcols) (xmin, xmax), (ymin, ymax) = tp # e.g. [(0, 1999), (0, 999)] at a blocksize of 2000*1000 (rowsxcols)
......
...@@ -28,6 +28,8 @@ from ..config import set_config, GMS_config ...@@ -28,6 +28,8 @@ from ..config import set_config, GMS_config
from .multiproc import MAP from .multiproc import MAP
from ..misc.definition_dicts import proc_chain, db_jobs_statistics_def from ..misc.definition_dicts import proc_chain, db_jobs_statistics_def
from py_tools_ds.numeric.array import get_array_tilebounds
if TYPE_CHECKING: if TYPE_CHECKING:
from collections import OrderedDict # noqa F401 # flake8 issue from collections import OrderedDict # noqa F401 # flake8 issue
...@@ -327,8 +329,8 @@ class process_controller(object): ...@@ -327,8 +329,8 @@ class process_controller(object):
else: else:
# define tile positions and size # define tile positions and size
def get_tilepos_list(GMSfile): def get_tilepos_list(GMSfile):
return HLP_F.get_image_tileborders('block', blocksize, return get_array_tilebounds(array_shape=INP_R.GMSfile2dict(GMSfile)['shape_fullArr'],
shape_fullArr=INP_R.GMSfile2dict(GMSfile)['shape_fullArr']) tile_shape=blocksize)
# get input parameters for creating GMS objects as blocks # get input parameters for creating GMS objects as blocks
work = [[GMSfile, ['block', tp]] for GMSfile in GMSfile_list_prevLvl_inDB work = [[GMSfile, ['block', tp]] for GMSfile in GMSfile_list_prevLvl_inDB
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment