input_reader.py 25.1 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
3
4

# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data
#
5
# Copyright (C) 2020  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
6
7
8
9
10
11
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
12
# the terms of the GNU General Public License as published by the Free Software
13
14
# Foundation, either version 3 of the License, or (at your option) any later version.
# Please note the following exception: `gms_preprocessing` depends on tqdm, which
15
16
17
# 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.
18
19
20
21
22
23
24
25
26
#
# 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 <http://www.gnu.org/licenses/>.
# -*- coding: utf-8 -*-
27
28
"""GeoMultiSens Input Reader:  Universal reader for all kinds of GeoMultiSens intermediate results."""

29
import collections
30
import fnmatch
31
32
import json
import os
33
import re
34
35
import tarfile
import zipfile
36
from tempfile import NamedTemporaryFile as tempFile
Daniel Scheffler's avatar
Daniel Scheffler committed
37
from logging import Logger
38
from typing import List, Tuple, TYPE_CHECKING  # noqa F401  # flake8 issue
39
from datetime import datetime
40
import logging
41
42
43
44
45

import dill
import numpy as np
import spectral
from spectral.io import envi as envi
46
from pyproj import CRS
47
48
49

from geoarray import GeoArray
from py_tools_ds.geo.coord_calc import corner_coord_to_minmax
50
from py_tools_ds.geo.vector.geometry import boxObj
51
from py_tools_ds.geo.coord_trafo import transform_any_prj
52
from py_tools_ds.geo.projection import isProjectedOrGeographic
53
from pyrsr.rsr import RSR_reader
54

55
from ..options.config import GMS_config as CFG
56
57
from ..misc import path_generator as PG
from ..misc import helper_functions as HLP_F
58
from ..misc.logging import GMS_logger, close_logger
59
from ..misc.database_tools import get_overlapping_scenes_from_postgreSQLdb
60
from ..misc.path_generator import path_generator
61
from ..misc.spatial_index_mediator import SpatialIndexMediator
62
from ..misc.locks import DatabaseLock
63

64
65
66
if TYPE_CHECKING:
    from ..model.gms_object import GMS_identifier  # noqa F401  # flake8 issue

67

68
def read_ENVIfile(path, arr_shape, arr_pos, logger=None, return_meta=False, q=0):
69
    hdr_path = os.path.splitext(path)[0] + '.hdr' if not os.path.splitext(path)[1] == '.hdr' else path
70
    return read_ENVI_image_data_as_array(hdr_path, arr_shape, arr_pos, logger=logger, return_meta=return_meta, q=q)
71

72

73
def read_ENVIhdr_to_dict(hdr_path, logger=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
74
    # type: (str, Logger) -> dict
75
    if not os.path.isfile(hdr_path):
76
77
78
79
80
        if logger:
            logger.critical('read_ENVIfile: Input data not found at %s.' % hdr_path)
        else:
            print('read_ENVIfile: Input data not found at %s.' % hdr_path)
        raise FileNotFoundError(hdr_path)
81
82
83
84
    else:
        SpyFileheader = envi.open(hdr_path)
        return SpyFileheader.metadata

85

86
def read_ENVI_image_data_as_array(path, arr_shape, arr_pos, logger=None, return_meta=False, q=0):
87
    """Read ENVI image data as array using a specified read pattern.
88

89
    :param path:        <str> Path of the ENVI image or header file
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    :param arr_shape:   <str> 'cube','row','col','band','block','pixel' or 'custom'
    :param arr_pos:     None, <int> or <list>. The content of the list depends on the chosen arr_shape as follows:
                        - 'cube':   No array position neccessary. Set arr_pos to None.
                        - 'row':    single int. -> Row with this index will be read.
                        - 'col':    single int. -> Column with this index will be read.
                        - 'band':   single int. -> Band with this index will be read.
                        - 'block':  row_bounds (2-tuple of ints): (a, b) -> Rows a through b-1 will be read.
                                    col_bounds (2-tuple of ints): (a, b) -> Columns a through b-1 will be read.
                        - 'pixel':  row, col (int): Indices of the row & column for the pixel to be read.
                        - 'custom': row_bounds (2-tuple of ints): (a, b) -> Rows a through b-1 will be read.
                                    col_bounds (2-tuple of ints): (a, b) -> Columns a through b-1 will be read.
                                    bands (list of ints):  Optional list of bands to read.
                                                            If not specified, all bands are read.
    :param logger:      <instance> of logging.logger (optional)
    :param return_meta: <bool> whether to return not only raster data but also meta data (optional, default=False)
    :param q:           <bool> quiet mode (supresses all console or logging output) (optional, default=False)
106
    """
107
    hdr_path = os.path.splitext(path)[0] + '.hdr' if not os.path.splitext(path)[1] == '.hdr' else path
Daniel Scheffler's avatar
Daniel Scheffler committed
108
    if not os.path.isfile(hdr_path):
109
        if not q:
110
111
112
113
            if logger:
                logger.critical('read_ENVIfile: Input data not found at %s.' % hdr_path)
            else:
                print('read_ENVIfile: Input data not found at %s.' % hdr_path)
Daniel Scheffler's avatar
Daniel Scheffler committed
114
    else:
115
116
117
118
119
120
        assert arr_shape in ['cube', 'row', 'col', 'band', 'block', 'pixel', 'custom'], \
            "Array shape '%s' is not known. Known array shapes are cube, row, col, band, block, pixel, custom." \
            % arr_shape
        if logger and not q:
            logger.info('Reading %s ...' % (os.path.basename(hdr_path)))
        File_obj = spectral.open_image(hdr_path)
Daniel Scheffler's avatar
Daniel Scheffler committed
121
        SpyFileheader = envi.open(hdr_path)
122
123
124
125
126
127
128
129
130
        image_data = File_obj.read_bands(range(SpyFileheader.nbands)) if arr_shape == 'cube' else \
            File_obj.read_subimage([arr_pos], range(SpyFileheader.ncols)) if arr_shape == 'row' else \
            File_obj.read_subimage(range(SpyFileheader.nrows), [arr_pos]) if arr_shape == 'col' else \
            File_obj.read_band(arr_pos) if arr_shape == 'band' else \
            File_obj.read_subregion((arr_pos[0][0], arr_pos[0][1] + 1),
                                    (arr_pos[1][0], arr_pos[1][1] + 1)) if arr_shape == 'block' else \
            File_obj.read_pixel(arr_pos[0], arr_pos[1]) if arr_shape == 'pixel' else \
            File_obj.read_subregion((arr_pos[0][0], arr_pos[0][1] + 1),
                                    (arr_pos[1][0], arr_pos[1][1] + 1), arr_pos[2])  # custom
131
        return (image_data, SpyFileheader.metadata) if return_meta else image_data
Daniel Scheffler's avatar
Daniel Scheffler committed
132

133

134
135
136
def read_mask_subset(path_masks, bandname, logger, subset=None):
    subset = subset if subset else ['cube', None]
    assert subset[0] in ['cube', 'block', 'MGRS_tile'], \
137
        "INP_R.read_mask_subset(): '%s' subset is not supported." % subset[0]
138
139
140
141
142
    path_masks_hdr = \
        os.path.splitext(path_masks)[0] + '.hdr' if not os.path.splitext(path_masks)[1] == '.hdr' else path_masks
    hdrDict = read_ENVIhdr_to_dict(path_masks_hdr, logger)
    (rS, rE), (cS, cE) = ((0, hdrDict['lines']), (0, hdrDict['samples'])) if subset[0] == 'cube' else subset[1]
    band_idx = hdrDict['band names'].index(bandname) if bandname in hdrDict['band names'] else None
143
    if band_idx is None:
144
        logger.warning("No band called '%s' in %s. Attribute set to <None>." % (bandname, path_masks))
145
        mask_sub = None
146
    else:
147
148
149
150
151
        if subset is None or subset[0] == 'cube':
            mask_sub = read_ENVI_image_data_as_array(path_masks_hdr, 'band', band_idx, logger=logger, q=1)
        else:
            mask_sub = read_ENVI_image_data_as_array(
                path_masks_hdr, 'custom', ((rS, rE), (cS, cE), [band_idx]), logger=logger, q=1)
152
        mask_sub = mask_sub[:, :, 0] if len(mask_sub.shape) == 3 and mask_sub.shape[2] == 1 else mask_sub
153
154
    return mask_sub

155

Daniel Scheffler's avatar
Daniel Scheffler committed
156
def GMSfile2dict(path_GMSfile):
157
    """ Converts a JSON file (like the GMS file) to a Python dictionary with keys and values.
158

159
160
    :param path_GMSfile:    absolute path on disk
    :return:                the corresponding Python dictionary
161
    """
162
163
164
    with open(path_GMSfile) as inF:
        GMSdict = json.load(inF)
    return GMSdict
Daniel Scheffler's avatar
Daniel Scheffler committed
165

166

167
168
def unify_envi_header_keys(header_dict):
    """Ensures the compatibility of ENVI header keys written by Spectral-Python the code internal attribute names.
169
    (ENVI header keys are always lowercase in contrast to the attribute names used in code).
170

171
172
    :param header_dict:
    """
173
    refkeys = ['AcqDate', 'AcqTime', 'Additional', 'FieldOfView', 'IncidenceAngle', 'Metafile', 'PhysUnit',
174
               'ProcLCode', 'Quality', 'Satellite', 'Sensor', 'SunAzimuth', 'SunElevation', 'ViewingAngle']
Daniel Scheffler's avatar
Daniel Scheffler committed
175
176
177
    unified_header_dict = header_dict
    for key in header_dict.keys():
        for refkey in refkeys:
178
            if re.match(key, refkey, re.I) and key != refkey:
Daniel Scheffler's avatar
Daniel Scheffler committed
179
180
181
182
                unified_header_dict[refkey] = header_dict[key]
                del unified_header_dict[key]
    return unified_header_dict

183

184
185
def get_list_GMSfiles(dataset_list, target):
    """Returns a list of absolute paths pointing to gms-files of truely written datasets that fullfill certain criteria.
186

187
188
189
    :param dataset_list: [dataset1_dictionary, dataset2_dictionary]
    :param target:       target GMS processing level
    :return              [/path/to/gms_file1.gms, /path/to/gms_file1.gms]
190
    """
191
    dataset_list = [dataset_list] if not isinstance(dataset_list, list) else dataset_list
192

193
194
    def get_gmsP(ds, tgt):
        return PG.path_generator(ds, proc_level=tgt).get_path_gmsfile()
195

196
    GMS_list = [p for p in [get_gmsP(ds, target) for ds in dataset_list] if os.path.exists(p)]
197

Daniel Scheffler's avatar
Daniel Scheffler committed
198
    return GMS_list
199
200


201
def SRF_Reader(GMS_id, no_thermal=None, no_pan=None, v=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
202
    # type: (GMS_identifier, bool, bool, bool) -> collections.OrderedDict
203
    """Read SRF for any sensor and return a dictionary containing band names as keys and SRF numpy arrays as values.
204

205
    :param GMS_id:
206
207
208
209
210
    :param no_thermal:      whether to exclude thermal bands from the returned bands list
                            (default: CFG.skip_thermal)
    :param no_pan:          whether to exclude panchromatic bands from the returned bands list
                            (default: CFG.skip_pan)
    :param v:               verbose mode
211
    """
212
213
214
215
    # set defaults
    # NOTE: these values cannot be set in function signature because CFG is not present at library import time
    no_thermal = no_thermal if no_thermal is not None else CFG.skip_thermal
    no_pan = no_pan if no_pan is not None else CFG.skip_pan
216

217
    logger = GMS_id.logger or Logger(__name__)
Daniel Scheffler's avatar
Daniel Scheffler committed
218

219
220
221
222
223
224
225
226
227
228
    SRF_dict = RSR_reader(satellite=GMS_id.satellite,
                          sensor=GMS_id.sensor,
                          subsystem=GMS_id.subsystem,
                          no_thermal=no_thermal,
                          no_pan=no_pan,
                          after_ac=False,
                          sort_by_cwl=True,
                          tolerate_missing=True,
                          logger=logger,
                          v=v)
229

230
    return SRF_dict
Daniel Scheffler's avatar
Daniel Scheffler committed
231

232

233
def pickle_SRF_DB(L1A_Instances, dir_out):
Daniel Scheffler's avatar
Daniel Scheffler committed
234
235
    list_GMS_identifiers = [i.GMS_identifier for i in L1A_Instances]
    out_dict = collections.OrderedDict()
236
    logger = GMS_logger('log__SRF2PKL', path_logfile=os.path.join(dir_out, 'log__SRF2PKL.log'),
237
                        log_level=CFG.log_level, append=False)
238
    for Id, Inst in zip(list_GMS_identifiers, L1A_Instances):
Daniel Scheffler's avatar
Daniel Scheffler committed
239
        Id['logger'] = logger
240
241
        out_dict[
            Inst.satellite + '_' + Inst.sensor + (('_' + Inst.subsystem) if Inst.subsystem not in ['', None] else '')] \
242
            = SRF_Reader(Id)
243
    print(list(out_dict.keys()))
244
    outFilename = os.path.join(dir_out, 'SRF_DB.pkl')
Daniel Scheffler's avatar
Daniel Scheffler committed
245
    with open(outFilename, 'wb') as outFile:
246
247
        dill.dump(out_dict, outFile)
        print('Saved SRF dictionary to %s' % outFilename)
248
249
250
    # with open(outFilename, 'rb') as inFile:
    #     readFile = pickle.load(inFile)
    # [print(i) for i in readFile.keys()]
251
    logger.close()
Daniel Scheffler's avatar
Daniel Scheffler committed
252

253

254
def Solar_Irradiance_reader(resol_nm=None, wvl_min_nm=None, wvl_max_nm=None):
255
256
257
258
259
260
261
    """Read the solar irradiance file and return an array of irradiances.

    :param resol_nm:    spectral resolution for returned irradiances [nanometers]
    :param wvl_min_nm:  minimum wavelength of returned irradiances [nanometers]
    :param wvl_max_nm:  maximum wavelength of returned irradiances [nanometers]
    :return:
    """
262
263
    from scipy.interpolate import interp1d

264
    sol_irr = np.loadtxt(CFG.path_solar_irr, skiprows=1)
265
266
267
268
    if resol_nm is not None and isinstance(resol_nm, (float, int)):
        wvl_min = (np.min(sol_irr[:, 0]) if wvl_min_nm is None else wvl_min_nm)
        wvl_max = (np.max(sol_irr[:, 0]) if wvl_max_nm is None else wvl_max_nm)
        wvl_rsp = np.arange(wvl_min, wvl_max, resol_nm)
269
        sol_irr = interp1d(sol_irr[:, 0], sol_irr[:, 1], kind='linear')(wvl_rsp)
270
    return sol_irr
271
272


273
274
275
276
277
278
279
def open_specific_file_within_archive(path_archive, matching_expression, read_mode='r'):
    # type: (str, str, str) -> (str, str)
    """Finds a specific file within an archive using a given matching expression and returns its content as string.

    :param path_archive:        the file path of the archive
    :param matching_expression: the matching expession to find the file within the archive
    :param read_mode:           the read mode used to open the archive (default: 'r')
Daniel Scheffler's avatar
Daniel Scheffler committed
280
    :return: tuple(content_file, filename_file)
281
282
    """

283
284
    file_suffix = os.path.splitext(path_archive)[1][1:]
    file_suffix = 'tar.gz' if path_archive.endswith('tar.gz') else file_suffix
285
    assert file_suffix in ['zip', 'tar', 'gz', 'tgz', 'tar.gz'], '*.%s files are not supported.' % file_suffix
286

287
    if file_suffix == 'zip':
288
289
290
        archive = zipfile.ZipFile(path_archive, 'r')
        # [print(i) for i in archive.namelist()]
        matching_files = fnmatch.filter(archive.namelist(), matching_expression)
291
292
293
294
295
296
297

        # NOTE: flink deactivates assertions via python -O flag. So, a usual 'assert matching_files' does NOT work here.
        if not matching_files:
            raise RuntimeError('Matching expression matches no file. Please revise your expression!')
        if len(matching_files) > 1:
            raise RuntimeError('Matching expression matches more than 1 file. Please revise your expression!')

298
299
        content_file = archive.read(matching_files[0])
        filename_file = os.path.join(path_archive, matching_files[0])
300

301
302
    else:  # 'tar','gz','tgz', 'tar.gz'
        archive = tarfile.open(path_archive, 'r|gz')  # open in stream mode is much faster than normal mode
303
304
        count_matching_files = 0
        for F in archive:
305
306
307
308
            if fnmatch.fnmatch(F.name, matching_expression):
                content_file = archive.extractfile(F)
                content_file = content_file.read()
                filename_file = os.path.join(path_archive, F.name)
309
                count_matching_files += 1
310
311
312
313

        # NOTE: flink deactivates assertions via python -O flag. So, a usual 'assert matching_files' does NOT work here.
        if count_matching_files == 0:
            raise RuntimeError('Matching expression matches no file. Please revise your expression!')
Daniel Scheffler's avatar
Daniel Scheffler committed
314
        if count_matching_files > 1:
315
            raise RuntimeError('Matching expression matches more than 1 file. Please revise your expression!')
316

317
    archive.close()
318
319
    content_file = content_file.decode('latin-1') \
        if isinstance(content_file, bytes) and read_mode == 'r' else content_file  # Python3
320

321
322
323
    return content_file, filename_file


324
class DEM_Creator(object):
325
    def __init__(self, dem_sensor='SRTM', db_conn='', logger=None):
326
327
328
329
330
        """Creator class for digital elevation models based on ASTER or SRTM.

        :param dem_sensor:  'SRTM' or 'ASTER'
        :param db_conn:     database connection string
        """
331
        if dem_sensor not in ['SRTM', 'ASTER']:
Daniel Scheffler's avatar
Daniel Scheffler committed
332
333
            raise ValueError('%s is not a supported DEM sensor. Choose between SRTM and ASTER (both 30m native GSD).'
                             % dem_sensor)
334
335

        self.dem_sensor = dem_sensor
336
        self.db_conn = db_conn if db_conn else CFG.conn_database
337
        self.logger = logger or logging.getLogger('DEM_Creator')
338
339
340
341
342
343
344

        self.project_dir = os.path.abspath(os.path.curdir)
        self.rootpath_DEMs = ''
        self.imNames = []
        self.dsID_dic = dict(ASTER=2, SRTM=225)
        self.DEM = None

345
346
347
348
349
350
351
352
353
354
355
356
    def __getstate__(self):
        """Defines how the attributes of DEM_Creator are pickled."""

        if self.logger not in [None, 'not set']:
            close_logger(self.logger)
            self.logger = None
        return self.__dict__

    def __del__(self):
        close_logger(self.logger)
        self.logger = None

357
358
359
360
361
    @staticmethod
    def _get_corner_coords_lonlat(cornerCoords_tgt, prj):
        # transform to Longitude/Latitude coordinates
        tgt_corner_coord_lonlat = [transform_any_prj(prj, 4326, X, Y) for X, Y in cornerCoords_tgt]

362
        # handle coordinates crossing the 180 degrees meridian
363
364
365
366
367
368
        xvals = [x for x, y in tgt_corner_coord_lonlat]
        if max(xvals) - min(xvals) > 180:
            tgt_corner_coord_lonlat = [(x, y) if x > 0 else (x + 360, y) for x, y in tgt_corner_coord_lonlat]

        return tgt_corner_coord_lonlat

369
370
371
    def get_overlapping_DEM_tiles(self, tgt_corner_coord_lonlat, timeout_sec=30, use_index_mediator=True):
        # type: (List[Tuple], int, bool) -> list
        """Get the overlapping DEM tiles for the given extent.
372
373

        :param tgt_corner_coord_lonlat: list of longitude/latitude target coordinates [[X,Y], [X,Y], ...]]
374
375
376
        :param timeout_sec:             database query timeout (seconds)
        :param use_index_mediator:      whether to use or not to use the SpatialIndexMediator (default: True)
        :return: list of matching DEM tile scene IDs
377
        """
378
        if use_index_mediator:
379
            SpIM = SpatialIndexMediator(host=CFG.spatial_index_server_host, port=CFG.spatial_index_server_port,
380
                                        timeout=timeout_sec, retries=10)
381
382
383
384
385
386
            with DatabaseLock(allowed_slots=1, logger=self.logger):
                scenes = SpIM.getFullSceneDataForDataset(envelope=corner_coord_to_minmax(tgt_corner_coord_lonlat),
                                                         timeStart=datetime(1970, 1, 1, 0, 0, 0),
                                                         timeEnd=datetime(2100, 12, 31, 0, 0, 0),
                                                         minCloudCover=0, maxCloudCover=100,
                                                         datasetid=self.dsID_dic[self.dem_sensor])
387
388
            sceneIDs_srtm = [scene.sceneid for scene in scenes]

389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
        else:
            sceneIDs_srtm = get_overlapping_scenes_from_postgreSQLdb(self.db_conn,
                                                                     table='scenes',
                                                                     tgt_corners_lonlat=tgt_corner_coord_lonlat,
                                                                     conditions=['datasetid=%s'
                                                                                 % self.dsID_dic[self.dem_sensor]],
                                                                     timeout=timeout_sec*1000)  # milliseconds
            sceneIDs_srtm = [i[0] for i in sceneIDs_srtm]

        return sceneIDs_srtm

    def _get_DEMPathes_to_include(self, tgt_corner_coord_lonlat, timeout_sec=30):
        # type: (List[Tuple], int) -> list
        """Return the paths of DEM files to merge in order to generate a DEM covering the given area of interest.

        :param tgt_corner_coord_lonlat:     list of longitude/latitude target coordinates [(X,Y), (X,Y), ...]]
        :param timeout_sec:                 database query timeout (seconds)
        :return: list of GDAL readable paths
        """
        # get overlapping SRTM scene IDs from GMS database
        try:
            # try to use the SpatialIndexMediator
411
            # noinspection PyBroadException
412
413
414
415
            try:
                sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, timeout_sec)
            except ConnectionRefusedError:
                # fallback to plain pgSQL
416
                self.logger.warning('SpatialIndexMediator refused connection. Falling back to plain postgreSQL query.')
417
                sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, use_index_mediator=False)
418
419
420
421
422
423
            except Exception as err:
                # fallback to plain pgSQL
                self.logger.warning('Error while running SpatialIndexMediator database query. '
                                    'Falling back to plain postgreSQL query. '
                                    'Error message was: %s' % str(repr(err)))
                sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, use_index_mediator=False)
424

425
            if not sceneIDs_srtm:
426
                # fallback to plain pgSQL
427
428
                self.logger.warning('SpatialIndexMediator did not return matching DEM tiles. '
                                    'Trying plain postgreSQL query..')
429
                sceneIDs_srtm = self.get_overlapping_DEM_tiles(tgt_corner_coord_lonlat, use_index_mediator=False)
430
431

        except TimeoutError:
432
433
            raise TimeoutError('Spatial database query for %s DEM generation timed out after %s seconds.'
                               % (self.dem_sensor, timeout_sec))
434

435
        if not sceneIDs_srtm:
436
            raise RuntimeError('No matching %s scene IDs for DEM generation found.' % self.dem_sensor)
437
438
439
440

        # get corresponding filenames on disk and generate GDALvsiPathes pointing to raster files within archives
        archivePaths = [path_generator(scene_ID=ID).get_local_archive_path_baseN() for ID in sceneIDs_srtm]
        self.rootpath_DEMs = os.path.dirname(archivePaths[0])
441
        imNameMatchExp = '*.bil' if self.dem_sensor == 'SRTM' else '*dem.tif'
442
443
        self.imNames = [fnmatch.filter(HLP_F.get_zipfile_namelist(aP), imNameMatchExp)[0] for aP in archivePaths]
        gdalImPaths = [os.path.join(HLP_F.convert_absPathArchive_to_GDALvsiPath(aP), bN)
444
                       for aP, bN in zip(archivePaths, self.imNames)]
445
446
447

        return gdalImPaths

448
    def _run_cmd(self, cmd):
449
450
        output, exitcode, err = HLP_F.subcall_with_output(cmd)
        if exitcode:
451
452
453
            self.logger.error('\nAn error occurred during DEM creation.')
            self.logger.warning("Print traceback in case you care:")
            self.logger.warning(err.decode('latin-1'))
454
455
456
457
458
459
        if output:
            return output.decode('UTF-8')

    def from_extent(self, cornerCoords_tgt, prj, tgt_xgsd, tgt_ygsd):
        """Returns a GeoArray of a DEM according to the given target coordinates

460
        :param cornerCoords_tgt:    list of target coordinates [[X,Y], [X,Y], ...]] (at least 2 coordinates)
461
462
463
464
465
466
        :param prj:                 WKT string of the projection belonging cornerCoords_tgt
        :param tgt_xgsd:            output X GSD
        :param tgt_ygsd:            output Y GSD
        :return: DEM GeoArray
        """

467
468
469
470
471
        # generate at least 4 coordinates in case less coords have been given in order to avoid nodata triangles in DEM
        if len(cornerCoords_tgt) < 4 and isProjectedOrGeographic(prj) == 'projected':
            co_yx = [(y, x) for x, y in cornerCoords_tgt]
            cornerCoords_tgt = boxObj(boxMapYX=co_yx).boxMapXY

472
473
474
475
476
477
478
        # handle coordinate infos
        tgt_corner_coord_lonlat = self._get_corner_coords_lonlat(cornerCoords_tgt, prj)

        # get GDAL readable pathes
        pathes = self._get_DEMPathes_to_include(tgt_corner_coord_lonlat)

        # Build GDAL VRT from pathes and create output DEM
479
480
481
482
        if not os.path.exists(CFG.path_tempdir):
            os.makedirs(CFG.path_tempdir)  # directory where tempfile is created must exist (CentOS)
        with tempFile(dir=CFG.path_tempdir, prefix='GeoMultiSens_', suffix='_dem_merged.vrt') as tFm, \
                tempFile(dir=CFG.path_tempdir, prefix='GeoMultiSens_', suffix='_dem_out.vrt') as tFo:
483
484
485
486
487
488
489
490
491
492
493
494

            try:
                os.chdir(self.rootpath_DEMs)

                # create a merged VRT of all input DEMs
                t_xmin, t_xmax, t_ymin, t_ymax = corner_coord_to_minmax(tgt_corner_coord_lonlat)
                self._run_cmd('gdalbuildvrt %s %s -te %s %s %s %s -vrtnodata 0'
                              % (tFm.name, ' '.join(pathes), t_xmin, t_ymin, t_xmax, t_ymax))

                # run gdalwarp to match target grid and extent
                merged_prj = GeoArray(tFm.name).prj
                t_xmin, t_xmax, t_ymin, t_ymax = corner_coord_to_minmax(cornerCoords_tgt)
495
                self._run_cmd('gdalwarp -r average -of VRT -srcnodata 0 -dstnodata 0 '
496
497
                              '-tr %s %s -s_srs EPSG:%s -t_srs EPSG:%s -te %s %s %s %s %s %s'
                              % (tgt_xgsd, tgt_ygsd, CRS(merged_prj).to_epsg(), CRS(prj).to_epsg(),
498
                                 t_xmin, t_ymin, t_xmax, t_ymax, tFm.name, tFo.name))
499
500
501
502
503
504
505
506
                assert os.path.exists(tFo.name)

                self.DEM = GeoArray(tFo.name).to_mem()

            finally:
                os.chdir(self.project_dir)

        return self.DEM