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

# enpt_enmapboxapp, A QGIS EnMAPBox plugin providing a GUI for the EnMAP processing tools (EnPT)
#
Daniel Scheffler's avatar
Daniel Scheffler committed
5
# Copyright (C) 2018-2021 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#
# 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 Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# 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/>.

"""This module provides the base class for EnPTAlgorithm and ExternalEnPTAlgorithm."""

import os
from os.path import expanduser
29
import psutil
30
from pkgutil import get_loader
31
32
from datetime import date
from multiprocessing import cpu_count
33
34
35
from threading import Thread
from queue import Queue
from subprocess import Popen, PIPE
36
from glob import glob
37

38
39
40
41
42
43
44
45
from qgis.core import (
    QgsProcessingAlgorithm,
    QgsProcessingParameterFile,
    QgsProcessingParameterNumber,
    QgsProcessingParameterFolderDestination,
    QgsProcessingParameterBoolean,
    QgsProcessingParameterDefinition,
    QgsProcessingParameterRasterLayer,
46
47
    QgsProcessingParameterEnum,
    NULL
48
)
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

from .version import __version__


class _EnPTBaseAlgorithm(QgsProcessingAlgorithm):
    # NOTE: The parameter assignments made here follow the parameter names in enpt/options/options_schema.py

    # Input parameters
    P_json_config = 'json_config'
    P_CPUs = 'CPUs'
    P_path_l1b_enmap_image = 'path_l1b_enmap_image'
    P_path_l1b_enmap_image_gapfill = 'path_l1b_enmap_image_gapfill'
    P_path_dem = 'path_dem'
    P_average_elevation = 'average_elevation'
    P_output_dir = 'output_dir'
    P_working_dir = 'working_dir'
    P_n_lines_to_append = 'n_lines_to_append'
66
    P_drop_bad_bands = 'drop_bad_bands'
67
    P_disable_progress_bars = 'disable_progress_bars'
68
    P_output_format = 'output_format'
69
70
71
72
73
74
75
    P_path_earthSunDist = 'path_earthSunDist'
    P_path_solar_irr = 'path_solar_irr'
    P_scale_factor_toa_ref = 'scale_factor_toa_ref'
    P_enable_keystone_correction = 'enable_keystone_correction'
    P_enable_vnir_swir_coreg = 'enable_vnir_swir_coreg'
    P_path_reference_image = 'path_reference_image'
    P_enable_ac = 'enable_ac'
76
77
78
79
80
    P_mode_ac = 'mode_ac'
    P_polymer_root = 'polymer_root'
    P_threads = 'threads'
    P_blocksize = 'blocksize'
    P_auto_download_ecmwf = 'auto_download_ecmwf'
81
82
83
84
85
86
87
88
    P_scale_factor_boa_ref = 'scale_factor_boa_ref'
    P_run_smile_P = 'run_smile_P'
    P_run_deadpix_P = 'run_deadpix_P'
    P_deadpix_P_algorithm = 'deadpix_P_algorithm'
    P_deadpix_P_interp_spectral = 'deadpix_P_interp_spectral'
    P_deadpix_P_interp_spatial = 'deadpix_P_interp_spatial'
    P_ortho_resampAlg = 'ortho_resampAlg'
    P_vswir_overlap_algorithm = 'vswir_overlap_algorithm'
89
90
    P_target_projection_type = 'target_projection_type'
    P_target_epsg = 'target_epsg'
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

    # # Output parameters
    P_OUTPUT_RASTER = 'outraster'
    # P_OUTPUT_VECTOR = 'outvector'
    # P_OUTPUT_FILE = 'outfile'
    P_OUTPUT_FOLDER = 'outfolder'

    def group(self):
        return 'Pre-Processing'

    def groupId(self):
        return 'PreProcessing'

    def name(self):
        return 'EnPTAlgorithm'

    def displayName(self):
108
        return f'EnPT - EnMAP Processing Tool (v{__version__})'
109
110
111
112

    def createInstance(self, *args, **kwargs):
        return type(self)()

113
114
115
116
117
    @staticmethod
    def _get_default_polymer_root():
        try:
            path_polymer = os.path.abspath(
                os.path.join(os.path.dirname(get_loader("polymer").path), os.pardir))
Daniel Scheffler's avatar
Daniel Scheffler committed
118
        except AttributeError:
119
120
121
            path_polymer = ''
        return path_polymer

122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
    @staticmethod
    def _get_default_output_dir():
        userhomedir = expanduser('~')

        default_enpt_dir = \
            os.path.join(userhomedir, 'Documents', 'EnPT', 'Output') if os.name == 'nt' else\
            os.path.join(userhomedir, 'EnPT', 'Output')

        outdir_nocounter = os.path.join(default_enpt_dir, date.today().strftime('%Y%m%d'))

        counter = 1
        while os.path.isdir('%s__%s' % (outdir_nocounter, counter)):
            counter += 1

        return '%s__%s' % (outdir_nocounter, counter)

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    def addParameter(self, param, *args, advanced=False, **kwargs):
        """Add a parameter to the QgsProcessingAlgorithm.

        This overrides the parent method to make it accept an 'advanced' parameter.

        :param param:       the parameter to be added
        :param args:        arguments to be passed to the parent method
        :param advanced:    whether the parameter should be flagged as 'advanced'
        :param kwargs:      keyword arguments to be passed to the parent method
        """
        if advanced:
            param.setFlags(param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)

        super(_EnPTBaseAlgorithm, self).addParameter(param, *args, **kwargs)

153
    def initAlgorithm(self, configuration=None):
154
155
156
157
158
159
160
        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_json_config,
                description='Configuration JSON template file',
                behavior=QgsProcessingParameterFile.File,
                extension='json',
                defaultValue=None,
Daniel Scheffler's avatar
Daniel Scheffler committed
161
                optional=True))
162

163
164
165
166
167
        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_CPUs,
                description='Number of CPU cores to be used for processing',
                type=QgsProcessingParameterNumber.Integer,
168
                defaultValue=cpu_count(), minValue=0, maxValue=cpu_count()),
169
            advanced=True)
170

171
172
173
174
        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_path_l1b_enmap_image,
                description='EnMAP Level-1B image (zip-archive or root directory)'))
175

176
177
178
179
180
        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_path_l1b_enmap_image_gapfill,
                description='Adjacent EnMAP Level-1B image to be used for gap-filling (zip-archive or root directory)',
                optional=True),
181
            advanced=True)
182

183
184
185
186
187
188
189
190
191
192
193
194
195
196
        self.addParameter(
            QgsProcessingParameterRasterLayer(
                name=self.P_path_dem,
                description='Input path of digital elevation model in map or sensor geometry; GDAL compatible file '
                            'format \n(must cover the EnMAP L1B data completely if given in map geometry or must have '
                            'the same \npixel dimensions like the EnMAP L1B data if given in sensor geometry)',
                optional=True))

        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_average_elevation,
                description='Average elevation in meters above sea level \n'
                            '(may be provided if no DEM is available and ignored if DEM is given)',
                type=QgsProcessingParameterNumber.Integer,
197
                defaultValue=0),
198
            advanced=True)
199

200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
        self.addParameter(
            QgsProcessingParameterFolderDestination(
                name=self.P_output_dir,
                description='Output directory where processed data and log files are saved',
                defaultValue=self._get_default_output_dir(),
                optional=True))

        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_working_dir,
                description='Directory to be used for temporary files',
                behavior=QgsProcessingParameterFile.Folder,
                defaultValue=None,
                optional=True))

        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_n_lines_to_append,
                description='Number of lines to be added to the main image [if not given, use the whole imgap]',
                type=QgsProcessingParameterNumber.Integer,
                defaultValue=None,
                optional=True),
222
            advanced=True)
223

224
225
226
227
228
229
        self.addParameter(
            QgsProcessingParameterBoolean(
                name=self.P_drop_bad_bands,
                description='Do not include bad bands (water absorption bands 1358-1453 nm / 1814-1961 nm) '
                            'in the L2A product',
                defaultValue=True),
230
            advanced=True)
231

232
233
234
235
236
        self.addParameter(
            QgsProcessingParameterBoolean(
                name=self.P_disable_progress_bars,
                description='Disable all progress bars during processing',
                defaultValue=True),
237
            advanced=True)
238

239
240
241
242
243
244
        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_output_format,
                description="Output format (file format of all raster output files).",
                options=['GTiff', 'ENVI'],
                defaultValue=0),
245
            advanced=True)
246

247
        # output_interleave?
248

249
250
251
252
253
254
        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_path_earthSunDist,
                description='Input path of the earth sun distance model',
                defaultValue=None,
                optional=True),
255
            advanced=True)
256

257
258
259
260
261
262
        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_path_solar_irr,
                description='Input path of the solar irradiance model',
                defaultValue=None,
                optional=True),
263
            advanced=True)
264

265
266
267
268
269
        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_scale_factor_toa_ref,
                description='Scale factor to be applied to TOA reflectance result',
                type=QgsProcessingParameterNumber.Integer,
270
                defaultValue=10000, minValue=0),
271
            advanced=True)
272

273
274
275
276
277
278
279
280
281
282
283
        # self.addParameter(
        #     QgsProcessingParameterBoolean(
        #         name=self.P_enable_keystone_correction,
        #         description='Keystone correction',
        #         defaultValue=False))

        # self.addParameter(
        #     QgsProcessingParameterBoolean(
        #         name=self.P_enable_vnir_swir_coreg,
        #         description='VNIR/SWIR co-registration',
        #         defaultValue=False))
284
285
286
287

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                name=self.P_path_reference_image,
288
                description='Reference image for absolute co-registration.',
289
290
291
292
293
294
                defaultValue=None,
                optional=True))

        self.addParameter(
            QgsProcessingParameterBoolean(
                name=self.P_enable_ac,
295
                description='Enable atmospheric correction',
296
297
                defaultValue=True))

298
299
300
301
        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_mode_ac,
                description="Atmospheric correction mode",
302
                options=['land - SICOR is applied to land AND water',
303
                         'water - POLYMER is applied to water only; land is cleared ',
304
                         'combined - SICOR is applied to land and POLYMER to water'],
305
                defaultValue=2))
306
307
308
309
310
311

        self.addParameter(
            QgsProcessingParameterFile(
                name=self.P_polymer_root,
                description='Polymer root directory (that contains the subdirectory for ancillary data)',
                behavior=QgsProcessingParameterFile.Folder,
312
313
                defaultValue=self._get_default_polymer_root(),
                optional=True),
314
315
316
317
318
            advanced=True)

        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_threads,
319
                description='Number of threads for multiprocessing when running ACwater/Polymer \n'
320
321
                            "('0: no threads', '-1: automatic', '>0: number of threads')",
                type=QgsProcessingParameterNumber.Integer,
322
                defaultValue=-1, minValue=-1, maxValue=cpu_count()),
323
324
325
326
327
            advanced=True)

        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_blocksize,
328
                description='Block size for multiprocessing when running ACwater/Polymer',
329
                type=QgsProcessingParameterNumber.Integer,
330
                defaultValue=100, minValue=1),
331
332
333
334
335
            advanced=True)

        self.addParameter(
            QgsProcessingParameterBoolean(
                name=self.P_auto_download_ecmwf,
336
337
                description='Automatically download ECMWF data for atmospheric correction '
                            'of water surfaces in ACwater/Polymer',
338
339
340
341
342
343
344
345
                defaultValue=False),
            advanced=True)

        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_scale_factor_boa_ref,
                description='Scale factor to be applied to BOA reflectance result',
                type=QgsProcessingParameterNumber.Integer,
346
                defaultValue=10000, minValue=0),
347
            advanced=True)
348

349
350
351
352
353
        # self.addParameter(
        #     QgsProcessingParameterBoolean(
        #         name=self.P_run_smile_P,
        #         description='Smile detection and correction (provider smile coefficients are ignored)',
        #         defaultValue=False))
354
355
356
357
358
359
360
361
362
363
364
365
366
367

        self.addParameter(
            QgsProcessingParameterBoolean(
                name=self.P_run_deadpix_P,
                description='Dead pixel correction',
                defaultValue=True))

        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_deadpix_P_algorithm,
                description="Algorithm for dead pixel correction",
                options=['spectral', 'spatial'],
                defaultValue=0),
            advanced=True)
368

369
370
371
372
373
374
375
        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_deadpix_P_interp_spectral,
                description="Spectral interpolation algorithm to be used during dead pixel correction ",
                options=['linear', 'bilinear', 'cubic', 'spline'],
                defaultValue=0),
            advanced=True)
376

377
378
379
380
381
382
        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_deadpix_P_interp_spatial,
                description="Spatial interpolation algorithm to be used during dead pixel correction",
                options=['linear', 'bilinear', 'cubic', 'spline'],
                defaultValue=0),
383
            advanced=True)
384

385
386
387
388
389
390
        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_ortho_resampAlg,
                description="Ortho-rectification resampling algorithm",
                options=['nearest', 'bilinear', 'gauss'],
                defaultValue=2),
391
            advanced=True)
392

393
394
395
396
397
398
399
400
        self.addParameter(
            QgsProcessingParameterEnum(
                name=self.P_vswir_overlap_algorithm,
                description="Algorithm specifying how to deal with the spectral bands "
                            "in the VNIR/SWIR spectral overlap region",
                options=['VNIR and SWIR bands, order by wavelength', 'average VNIR and SWIR bands',
                         'VNIR bands only', 'SWIR bands only'],
                defaultValue=3),
401
            advanced=True)
402

403
404
405
406
407
408
        self.addParameter(
            QgsProcessingParameterEnum(
                self.P_target_projection_type,
                description='Projection type of the raster output files',
                options=['UTM', 'Geographic'],
                defaultValue=0),
409
            advanced=True)
410

411
412
413
414
415
416
417
        self.addParameter(
            QgsProcessingParameterNumber(
                name=self.P_target_epsg,
                description='Custom EPSG code of the target projection (overrides target_projection_type)',
                type=QgsProcessingParameterNumber.Integer,
                defaultValue=None,
                optional=True),
418
            advanced=True)
419

420
421
422
423
    # TODO:
    #   "target_coord_grid": "None"  /*custom target coordinate grid to which the output is resampled
    #                        ([x0, x1, y0, y1], e.g., [0, 30, 0, 30])*/

424
425
426
427
428
429
430
431
432
433
434
435
436
437
    @staticmethod
    def shortHelpString(*args, **kwargs):
        """Example:

        '<p>Here comes the HTML documentation.</p>' \
        '<h3>With Headers...</h3>' \
        '<p>and Hyperlinks: <a href="www.google.de">Google</a></p>'

        :param args:
        :param kwargs:
        """

        text = \
            '<p>General information about this EnMAP box app can be found ' \
438
439
440
441
442
443
444
            '<a href="https://enmap.git-pages.gfz-potsdam.de/GFZ_Tools_EnMAP_BOX/enpt_enmapboxapp/doc/">here</a>. ' \
            'For details, e.g., about all the algorithms implemented in EnPT, take a look at the ' \
            '<a href="https://enmap.git-pages.gfz-potsdam.de/GFZ_Tools_EnMAP_BOX/EnPT/doc/index.html">EnPT backend ' \
            'documentation</a>.</p>' \
            '<p>Type <i>enpt -h</i> into a shell to get further information about individual parameters or check out ' \
            'the <a href="https://enmap.git-pages.gfz-potsdam.de/GFZ_Tools_EnMAP_BOX/EnPT/doc/usage.html#' \
            'command-line-utilities">documentation</a>.</p>'
445
446
447
448
449
450
451
452
453

        return text

    def helpString(self):
        return self.shortHelpString()

    @staticmethod
    def helpUrl(*args, **kwargs):
        return 'https://enmap.git-pages.gfz-potsdam.de/GFZ_Tools_EnMAP_BOX/enpt_enmapboxapp/doc/'
454

455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
    @staticmethod
    def _get_preprocessed_parameters(parameters):
        # replace Enum parameters with corresponding strings (not needed in case of unittest)
        for n, opts in [
            ('output_format', {0: 'GTiff', 1: 'ENVI'}),
            ('mode_ac', {0: 'land', 1: 'water', 2: 'combined'}),
            ('deadpix_P_algorithm', {0: 'spectral', 1: 'spatial'}),
            ('deadpix_P_interp_spectral', {0: 'linear', 1: 'bilinear', 2: 'cubic', 3: 'spline'}),
            ('deadpix_P_interp_spatial', {0: 'linear', 1: 'bilinear', 2: 'cubic', 3: 'spline'}),
            ('ortho_resampAlg', {0: 'nearest', 1: 'bilinear', 2: 'gauss'}),
            ('vswir_overlap_algorithm', {0: 'order_by_wvl', 1: 'average', 2: 'vnir_only', 3: 'swir_only'}),
            ('target_projection_type', {0: 'UTM', 1: 'Geographic'}),
        ]:
            if isinstance(parameters[n], int):
                parameters[n] = opts[parameters[n]]

        # remove all parameters not to be forwarded to the EnPT CLI
        parameters = {k: v for k, v in parameters.items()
                      if k not in ['anaconda_root']
                      and v not in [None, NULL, 'NULL', '']}

        return parameters

478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
    @staticmethod
    def _run_cmd(cmd, qgis_feedback=None, **kwargs):
        """Execute external command and get its stdout, exitcode and stderr.

        Code based on: https://stackoverflow.com/a/31867499

        :param cmd: a normal shell command including parameters
        """
        def reader(pipe, queue):
            try:
                with pipe:
                    for line in iter(pipe.readline, b''):
                        queue.put((pipe, line))
            finally:
                queue.put(None)

        process = Popen(cmd, stdout=PIPE, stderr=PIPE, shell=True, **kwargs)
        q = Queue()
        Thread(target=reader, args=[process.stdout, q]).start()
        Thread(target=reader, args=[process.stderr, q]).start()

499
500
501
        stdout_qname = None
        stderr_qname = None

502
503
504
505
506
507
508
509
510
511
512
513
514
515
        # for _ in range(2):
        for source, line in iter(q.get, None):
            if qgis_feedback.isCanceled():
                # qgis_feedback.reportError('CANCELED')

                proc2kill = psutil.Process(process.pid)
                for proc in proc2kill.children(recursive=True):
                    proc.kill()
                proc2kill.kill()

                raise KeyboardInterrupt

            linestr = line.decode('latin-1').rstrip()
            # print("%s: %s" % (source, linestr))
516
517
518
519

            # source name seems to be platfor/environment specific, so grab it from dummy STDOUT/STDERR messages.
            if linestr == 'Connecting to EnPT STDOUT stream.':
                stdout_qname = source.name
520
                continue
521
522
            if linestr == 'Connecting to EnPT STDERR stream.':
                stderr_qname = source.name
523
                continue
524
525

            if source.name == stdout_qname:
526
                qgis_feedback.pushInfo(linestr)
527
528
529
            elif source.name == stderr_qname:
                qgis_feedback.reportError(linestr)
            else:
530
531
532
533
534
                qgis_feedback.reportError(linestr)

        exitcode = process.poll()

        return exitcode
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549

    def _handle_results(self, parameters: dict, feedback, exitcode: int) -> dict:
        success = False

        if exitcode:
            feedback.reportError("\n" +
                                 "=" * 60 +
                                 "\n" +
                                 "An exception occurred. Processing failed.")

        # list output dir
        if 'output_dir' in parameters:
            outdir = parameters['output_dir']
            outraster_matches = \
                glob(os.path.join(outdir, '*', '*SPECTRAL_IMAGE.TIF')) or \
550
551
552
                glob(os.path.join(outdir, '*', '*SPECTRAL_IMAGE.bsq')) or \
                glob(os.path.join(outdir, '*', '*SPECTRAL_IMAGE.bil')) or \
                glob(os.path.join(outdir, '*', '*SPECTRAL_IMAGE.bip'))
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
            outraster = outraster_matches[0] if len(outraster_matches) > 0 else None

            if os.path.isdir(outdir):
                if os.listdir(outdir):
                    feedback.pushInfo("The output folder '%s' contains:\n" % outdir)
                    feedback.pushCommandInfo('\n'.join([os.path.basename(f) for f in os.listdir(outdir)]) + '\n')

                    if outraster:
                        subdir = os.path.dirname(outraster_matches[0])
                        feedback.pushInfo(subdir)
                        feedback.pushInfo("...where the folder '%s' contains:\n" % os.path.split(subdir)[-1])
                        feedback.pushCommandInfo('\n'.join(sorted([os.path.basename(f)
                                                                   for f in os.listdir(subdir)])) + '\n')
                        success = True
                    else:
                        feedback.reportError("No output raster was written.")

                else:
                    feedback.reportError("The output folder is empty.")

            else:
                feedback.reportError("No output folder created.")

            # return outputs
            if success:
                return {
                    'success': True,
                    self.P_OUTPUT_RASTER: outraster,
                    # self.P_OUTPUT_VECTOR: parameters[self.P_OUTPUT_RASTER],
                    # self.P_OUTPUT_FILE: parameters[self.P_OUTPUT_RASTER],
                    self.P_OUTPUT_FOLDER: outdir
                }
            else:
                return {'success': False}

        else:
            feedback.pushInfo('The output was skipped according to user setting.')
            return {'success': True}