Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / pycama / transform / Filters.py @ 838:750fb4349dae

History | View | Annotate | Download (40.1 KB)

1
#!/usr/bin/env python3
2
# -*- coding: utf-8 -*-
3

    
4
# Copyright 2016-2017 Maarten Sneep, KNMI
5
#
6
# Redistribution and use in source and binary forms, with or without
7
# modification, are permitted provided that the following conditions are met:
8
#
9
# 1. Redistributions of source code must retain the above copyright notice,
10
#    this list of conditions and the following disclaimer.
11
#
12
# 2. Redistributions in binary form must reproduce the above copyright notice,
13
#    this list of conditions and the following disclaimer in the documentation
14
#    and/or other materials provided with the distribution.
15
#
16
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
20
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23
# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26

    
27
## \file Filters.py
28
#  Filter functions for PyCAMA.
29
#
30
#  These subclasses of AbstractTransformer are filters available to PyCAMA by default.
31
#  User filters can be placed in a direcory pointed to by the `PYCAMA_PLUGIN_PATH`
32
#  environment variable (can be a list of directories). These can be accessed
33
#  through the same name space as the standard filters here, i.e. by specifying
34
#  <tt>transform._Classname_</tt> in the configuration.
35
# @author Maarten Sneep
36

    
37
import numpy as np
38
import numexpr
39

    
40
from .Transformer import AbstractTransformer
41
from ..utilities import CAMAException
42

    
43
## Unity transformation function for PyCAMA
44
#
45
# No-op filter for PyCAMA.
46
class Unity(AbstractTransformer):
47
    ## The actual callable interface
48
    #
49
    # This function simply returns the input array. The arguments are as defined in the superclass.
50
    def __call__(self, primary,
51
                 secondary,
52
                 latitude,
53
                 longitude,
54
                 solar_zenith_angle,
55
                 viewing_zenith_angle,
56
                 solar_azimuth_angle,
57
                 viewing_azimuth_angle,
58
                 processing_quality_flags,
59
                 geolocation_flags,
60
                 surface_classification):
61
        return primary
62

    
63

    
64
## Select only the central part of the swath in PyCAMA
65
#
66
# Cut off the edges of the swath.
67
class SwathPartSelection(AbstractTransformer):
68
    ## The actual callable interface
69
    #
70
    # This function simply returns the input array. The arguments are as defined in the superclass.
71
    def __init__(self, **kwargs):
72
        super().__init__()
73
        
74
        if 'count' in kwargs:
75
            self.count = kwargs['count']
76
        else:
77
            self.count = 200
78
        
79
    def __call__(self, primary,
80
                 secondary,
81
                 latitude,
82
                 longitude,
83
                 solar_zenith_angle,
84
                 viewing_zenith_angle,
85
                 solar_azimuth_angle,
86
                 viewing_azimuth_angle,
87
                 processing_quality_flags,
88
                 geolocation_flags,
89
                 surface_classification):
90
        
91
        n_pixels = primary.shape[0]
92
        offset = (n_pixels - self.count)//2
93
        try:
94
            mask = primary.mask
95
            dummy = mask[0]
96
        except (AttributeError, IndexError):
97
            primary = np.ma.asarray(primary)
98
            mask = np.zeros(primary.shape, dtype=np.bool)
99
        
100
        mask[0:offset, ...] = True
101
        mask[-offset:, ...] = True
102
        primary.mask = mask
103
        
104
        return primary
105

    
106
## %Unity transformation function for PyCAMA with extra debugging output.
107
#
108
# This no-op filter that prints out extra statistics on the data passing through it.
109
class DebugFilter(AbstractTransformer):
110
    ## The constructor
111
    #
112
    # @param[in] name A string label for printing, to identify the output if multiple variables use this filter.
113
    # @param[in] min  Count the number of pixels with a value smaller than this minimum.
114
    # @param[in] max  Count the number of pixels with a value larger than this maximum.
115
    def __init__(self, **kwargs):
116
        super().__init__()
117

    
118
        if 'name' in kwargs:
119
            ## The label for printing output.
120
            self.name = kwargs['name']
121
        else:
122
            self.name = ''
123

    
124
        if 'min' in kwargs:
125
            ## The minimum threshold. Report the number of observations that are below this threshold.
126
            self.minimum = kwargs['min']
127
        else:
128
            self.minimum = None
129

    
130
        if 'max' in kwargs:
131
            ## The maximum threshold. Report the number of observations that are above this threshold.
132
            self.maximum = kwargs['max']
133
        else:
134
            self.maximum = None
135

    
136
    ## The actual callable interface
137
    #
138
    # This function simply returns the input array after printing extra statistics
139
    # on the data passing through it. Note that standard out is used directly as the
140
    # logging system isn't available from within the filter sub-system.
141
    #
142
    # The arguments are as defined in the superclass.
143
    #
144
    def __call__(self, primary,
145
                 secondary,
146
                 latitude,
147
                 longitude,
148
                 solar_zenith_angle,
149
                 viewing_zenith_angle,
150
                 solar_azimuth_angle,
151
                 viewing_azimuth_angle,
152
                 processing_quality_flags,
153
                 geolocation_flags,
154
                 surface_classification):
155

    
156
        flattened = primary.flatten()
157
        if not hasattr(flattened, 'mask'):
158
            flattened = np.ma.asarray(flattened)
159
            flattened.mask = np.zeros(flattened.shape, dtype=np.bool)
160

    
161
        if 'min' in kwargs:
162
            select = flattened[np.logical_and(np.logical_not(flattened.mask), flattened < self.minimum)]
163
            if len(select) > 0:
164
                print("DebugFilter({3}): {0}/{2} values smaller than minimum {1}".format(len(select),
165
                    self.minimum, len(flattened), self.name))
166
        print("DebugFilter({1}: mimimum value {0}".format(np.min(flattened), self.name))
167

    
168
        if 'max' in kwargs:
169
            select = flattened[np.logical_and(np.logical_not(flattened.mask), flattened > self.maximum)]
170
            if len(select) > 0:
171
                print("DebugFilter({3}): {0}/{2} values larger than maximum {1}".format(len(select),
172
                    self.maximum, len(flattened), self.name))
173
        print("DebugFilter({1}): maximum value {0}".format(np.max(flattened), self.name))
174

    
175
        print("DebugFilter({0}): fill value count {1}/{2}".format(self.name, sum(flattened.mask), len(flattened)))
176

    
177
        return primary
178

    
179
## a surface filter
180
#
181
#  The filter removes pixels that do not match the set citeria.
182
class SurfaceFilter(AbstractTransformer):
183
    ## The constructor
184
    #
185
    #  @param select what to select: 'land' or 'water'
186
    #  @param source Allow other land/sea masks to be used. Currently accepts 'OMI', 'TROPOMI'/'S5P' (default)
187
    def __init__(self, **kwargs):
188
        super().__init__()
189

    
190
        if 'select' in kwargs:
191
            self.select = kwargs['select'].lower()
192
        else:
193
            self.select = 'land'
194

    
195
        if self.select not in ('land', 'water', 'both'):
196
            self.logger.warning("The select argument for SurfaceFilter the contains an illegal value")
197
            self.select = 'land'
198

    
199
        if 'source' in kwargs:
200
            self.source = kwargs['source'].upper()
201
            if self.source == 'TROPOMI':
202
                self.source = 'S5P'
203
        else:
204
            self.source = 'S5P'
205

    
206
        if self.source not in ('OMI', 'S5P'):
207
            self.logger.warning("The source argument for SurfaceFilter the contains an illegal value")
208
            self.source = 'S5P'
209

    
210
    ## land/sea filter for PyCAMA
211
    #
212
    # The arguments are as defined in the superclass.
213
    #
214
    def __call__(self,
215
                 primary,
216
                 secondary,
217
                 latitude,
218
                 longitude,
219
                 solar_zenith_angle,
220
                 viewing_zenith_angle,
221
                 solar_azimuth_angle,
222
                 viewing_azimuth_angle,
223
                 processing_quality_flags,
224
                 geolocation_flags,
225
                 surface_classification):
226

    
227
        if self.select == "both":
228
            return primary
229
        try:
230
            mask = primary.mask
231
            dummy = mask[0]
232
        except (AttributeError, IndexError):
233
            primary = np.ma.asarray(primary)
234
            mask = np.zeros(primary.shape, dtype=np.bool)
235

    
236
        if self.source == "S5P":
237
            if self.select == "land":
238
                mask = np.logical_or(mask, (surface_classification & 1) == 1)
239
            else:
240
                mask = np.logical_or(mask, (surface_classification & 1) == 0)
241

    
242
        elif self.source == "OMI":
243
            if self.select == "land":
244
                mask = np.logical_or(mask, (surface_classification & 0x000f) != 1)
245
            else:
246
                mask = np.logical_or(mask, (surface_classification & 0x000F) == 1)
247

    
248
        primary.mask = mask
249

    
250
        return primary
251

    
252
## A threshold filter for PyCAMA
253
#
254
#  The filter removes values on the 'wrong' side of the threshold (single sided).
255
class ThresholdFilter(AbstractTransformer):
256
    ## The constructor
257
    #
258
    #  @param threshold The cut off value for the filter
259
    #  @param comparison Binary comparison, resolving to 'True' for elements
260
    #    that must be removed, i.e. use '<' to remove values smaller than
261
    #    the threshold. This is a string which is used to build the comparison operation via `numexpr`.
262
    def __init__(self, **kwargs):
263
        super().__init__()
264

    
265
        if 'threshold' in kwargs:
266
            ## Threshold to apply in the filter.
267
            self.threshold = kwargs['threshold']
268
        else:
269
            self.threshold = 0.2
270

    
271
        if 'comparison' in kwargs:
272
            ## Comparison symbol for the threshold.
273
            self.comparison = kwargs['comparison']
274
        else:
275
            self.comparison = "<"
276

    
277
    ## Threshold filter for PyCAMA
278
    #
279
    # The arguments are as defined in the superclass.
280
    #
281
    def __call__(self,
282
                 primary,
283
                 secondary,
284
                 latitude,
285
                 longitude,
286
                 solar_zenith_angle,
287
                 viewing_zenith_angle,
288
                 solar_azimuth_angle,
289
                 viewing_azimuth_angle,
290
                 processing_quality_flags,
291
                 geolocation_flags,
292
                 surface_classification):
293

    
294
        # import into current namespace
295
        threshold = self.threshold
296

    
297
        if secondary is not None:
298
            mask = numexpr.evaluate("secondary {0} threshold".format(self.comparison))
299
        else:
300
            mask = numexpr.evaluate("primary {0} threshold".format(self.comparison))
301
        n_mask = np.sum(mask)
302
        p_mask = np.sum(primary.mask)
303

    
304
        self.logger.debug("{0}: {1} values masked ({2} before)".format(self.__class__.__name__, n_mask, p_mask))
305

    
306
        primary.mask = np.logical_or(mask, primary.mask)
307
        return primary
308

    
309
## A band-pass filter for PyCAMA
310
#
311
#  The filter removes values outside of the specified band.
312
class BandFilter(AbstractTransformer):
313
    ## The constructor
314
    #
315
    # @param lower The lower threshold value for the filter. The default is \f$-\infty\f$.
316
    # @param upper The upper threshold value for the filter. The default is \f$+\infty\f$.
317
    # @param parameter parameter on which to filter (string, one of 'primary', 'secondary',
318
    #            'latitude', 'longitude', 'solar_zenith_angle', 'viewing_zenith_angle',
319
    #             'solar_azimuth_angle', 'viewing_azimuth_angle', 'processing_quality_flags',
320
    #             'geolocation_flags', 'surface_classification'.
321
    #             Defaults to 'secondary if given, or primary otherwise).
322
    #
323
    # Keep in mind that the filter functions are never called on the base parameters
324
    # (basically the list of the parameter keyword argument).
325
    def __init__(self, **kwargs):
326
        super().__init__()
327

    
328
        if 'lower' in kwargs and kwargs['lower'] is not None:
329
            ## Lower limit of the band filter. Values below or equal to this threshold are removed.
330
            self.lower = kwargs['lower']
331
        else:
332
            self.lower = -np.inf
333

    
334
        if 'upper' in kwargs and kwargs['upper'] is not None:
335
            ## Upper limit of the band filter. Values above this threshold are removed.
336
            self.upper = kwargs['upper']
337
        else:
338
            self.upper = np.inf
339

    
340
        if 'parameter' in kwargs and kwargs['parameter'] is not None:
341
            ## The parameter to apply the filter to.
342
            self.parameter = kwargs['parameter']
343
        else:
344
            self.parameter = None
345

    
346
    ## A band-pass filter for PyCAMA
347
    #
348
    # The arguments are as defined in the superclass.
349
    def __call__(self,
350
                 primary,
351
                 secondary,
352
                 latitude,
353
                 longitude,
354
                 solar_zenith_angle,
355
                 viewing_zenith_angle,
356
                 solar_azimuth_angle,
357
                 viewing_azimuth_angle,
358
                 processing_quality_flags,
359
                 geolocation_flags,
360
                 surface_classification):
361

    
362
        lower = self.lower
363
        upper = self.upper
364

    
365
        if self.parameter is not None:
366
            parameter = self.parameter
367
        else:
368
            if secondary is not None:
369
                parameter = 'secondary'
370
            else:
371
                parameter = 'primary'
372

    
373
        mask = np.logical_not(np.logical_and(numexpr.evaluate("lower <= {0}".format(parameter)),
374
                                             numexpr.evaluate("{0} < upper".format(parameter))))
375

    
376
        if hasattr(primary.mask, 'shape'):
377
            primary.mask = np.logical_or(mask, primary.mask)
378
        else:
379
            primary = np.ma.asarray(primary)
380
            primary.mask = mask
381

    
382
        return primary
383

    
384
## Filter on (warning) flags.
385
#
386
#  Error flags lead to a fill value anyway, so those are irrelevant.
387
class ProcessingQualityFlagsFilter(AbstractTransformer):
388
    ## The constructor
389
    #
390
    #  @param[in] mask  The warning(s) that should be treated as errors.
391
    #         A bitwise AND of the values you want to remove from the `processing_quality_flags`
392
    #         variable in the S5P level 2 files. The default is to remove _all_ pixels with a warning.
393
    #
394
    #  Note that errors are automagically filtered out, as the primary
395
    #  already contains a fill value for those pixels.
396
    #
397
    def __init__(self, **kwargs):
398
        super().__init__()
399

    
400
        if 'mask' in kwargs and kwargs['mask'] is not None:
401
            ## mask, if the bitwise AND of this `mask` and the `processing_quality_flags` is not zero, the ground pixel is removed.
402
            self.mask = kwargs['mask']
403
        else:
404
            self.mask = 0xffffff00
405

    
406
    ## Processing quality flags filter for PyCAMA
407
    #
408
    #  The arguments are as defined in the superclass.
409
    #
410
    def __call__(self,
411
                 primary,
412
                 secondary,
413
                 latitude,
414
                 longitude,
415
                 solar_zenith_angle,
416
                 viewing_zenith_angle,
417
                 solar_azimuth_angle,
418
                 viewing_azimuth_angle,
419
                 processing_quality_flags,
420
                 geolocation_flags,
421
                 surface_classification):
422

    
423
        if processing_quality_flags is not None:
424
            outmask = np.asarray((processing_quality_flags & self.mask) != 0, dtype=np.bool)
425
            primary.mask = np.logical_or(outmask, primary.mask)
426

    
427
        return primary
428

    
429
## Apply a scaling to the variable.
430
#
431
#  Scale a variable by multiplication or division with a constant.
432
#  The `numexpr` module is used for the operation.
433
class Multiplier(AbstractTransformer):
434
    ## The constructor
435
    #
436
    #  @param[in] scalefactor The constant for the multiplication or division. Defaults to 100.0
437
    #  @param[in] operator    The mathematical operator to use, must be a string. Defaults to '*'.
438
    def __init__(self, **kwargs):
439
        super().__init__()
440

    
441
        if 'scalefactor' in kwargs:
442
            ## The constant for the multiplication or division.
443
            self.scalefactor = kwargs['scalefactor']
444
        else:
445
            self.scalefactor = 100.0
446

    
447
        if 'operator' in kwargs:
448
            ## The mathematical operator to use, a string
449
            self.operator = kwargs['operator']
450
        else:
451
            self.operator = "*"
452

    
453
    ## Multiplication filter.
454
    #
455
    #  The arguments are as defined in the superclass.
456
    #
457
    def __call__(self,
458
                 primary,
459
                 secondary,
460
                 latitude,
461
                 longitude,
462
                 solar_zenith_angle,
463
                 viewing_zenith_angle,
464
                 solar_azimuth_angle,
465
                 viewing_azimuth_angle,
466
                 processing_quality_flags,
467
                 geolocation_flags,
468
                 surface_classification):
469

    
470
        scalefactor = self.scalefactor
471
        mask = primary.mask
472
        primary = np.where(primary.mask, np.nan, primary)
473
        primary = numexpr.evaluate("primary {0} scalefactor".format(self.operator))
474
        primary = np.ma.asarray(primary)
475
        primary.mask = mask
476
        return primary
477

    
478

    
479
## Reduce the dimensionality of a data set by selecting a slice from the full array. This is a PyCAMA filter.
480
#
481
#  After the filtering stage, all data sets are supposed to have the same dimensionality,
482
#  this is one of the methods to acchieve that.
483
class Select(AbstractTransformer):
484
    ## Reduction of dimensionality filter for PyCAMA
485
    #
486
    # @param dimension The dimension (index) to slice. Default = -1 (last dimension).
487
    # @param index     Which index in the slice dimension to keep. Default = 0 (first item).
488
    # @param collapse  Instead of selection use mean along given dimension. This is a boolean flag.
489
    #
490
    def __init__(self, **kwargs):
491
        super().__init__()
492

    
493
        if 'dimension' in kwargs:
494
            ## The dimension in which to extract
495
            self.dimension = kwargs['dimension']
496
        else:
497
            self.dimension = None
498

    
499
        if 'index' in kwargs:
500
            ## the index in the dimension
501
            self.index = kwargs['index']
502
        else:
503
            self.index = None
504

    
505
        if 'collapse' in kwargs:
506
            ## Flag to indicate mean instead of slice
507
            self.collapse = kwargs['collapse']
508
        else:
509
            self.collapse = None
510

    
511
        if 'sum' in kwargs:
512
            ## Flag to indicate mean instead of slice
513
            self.sum = kwargs['sum']
514
        else:
515
            self.sum = None
516

    
517

    
518
    ## Select filter.
519
    #
520
    #  The arguments are as defined in the superclass.
521
    #
522
    def __call__(self,
523
                 primary,
524
                 secondary,
525
                 latitude,
526
                 longitude,
527
                 solar_zenith_angle,
528
                 viewing_zenith_angle,
529
                 solar_azimuth_angle,
530
                 viewing_azimuth_angle,
531
                 processing_quality_flags,
532
                 geolocation_flags,
533
                 surface_classification):
534

    
535
        if self.dimension is not None:
536
            dimension = self.dimension
537
        else:
538
            dimension = -1
539

    
540
        if self.index is not None:
541
            index = self.index
542
        else:
543
            index = 0
544

    
545
        if self.collapse is not None:
546
            collapse = self.collapse
547
        else:
548
            collapse = False
549

    
550
        if self.sum is not None:
551
            do_sum = self.collapse
552
        else:
553
            do_sum = False
554

    
555
        self.logger.debug("Select filter on dimension {0}, index {1}, {2}".format(dimension, index, "collapse" if collapse else "select"))
556

    
557
        if collapse:
558
            primary[np.logical_not(primary.mask)] = np.nan
559
            primary = np.nanmean(primary, axis=dimension)
560
            return primary
561

    
562
        if do_sum:
563
            primary[np.logical_not(primary.mask)] = np.nan
564
            primary = np.nansum(primary, axis=dimension)
565
            return primary
566

    
567
        if dimension == -1:
568
            if len(primary.shape) > 1:
569
                primary = primary[..., index]
570
            else:
571
                primary = primary[index]
572
        elif dimension == -2:
573
            if len(primary.shape) > 2:
574
                primary = primary[..., index, :]
575
            else:
576
                primary = primary[index, :]
577
        elif dimension == -3:
578
            if len(primary.shape) > 3:
579
                primary = primary[..., index, :, :]
580
            else:
581
                primary = primary[index, :, :]
582
        elif dimension == -4:
583
            if len(primary.shape) > 4:
584
                primary = primary[..., index, :, :, :]
585
            else:
586
                primary = primary[index, :, :, :]
587
        elif dimension == 1:
588
            if len(primary.shape) > 1:
589
                primary = primary[index, ...]
590
            else:
591
                primary = primary[index]
592
        elif dimension == 2:
593
            if len(primary.shape) > 2:
594
                primary = primary[:, index, ...]
595
            else:
596
                primary = primary[:, index]
597
        elif dimension == 3:
598
            if len(primary.shape) > 3:
599
                primary = primary[:, :, index, ...]
600
            else:
601
                primary = primary[:, :, index]
602
        elif dimension == 4:
603
            if len(primary.shape) > 4:
604
                primary = primary[:, :, :, index, ...]
605
            else:
606
                primary = primary[:, :, :, index]
607
        elif dimension == 5:
608
            if len(primary.shape) > 5:
609
                primary = primary[:, :, :, :, index, ...]
610
            else:
611
                primary = primary[:, :, :, :, index]
612
        elif dimension == 6:
613
            if len(primary.shape) > 6:
614
                primary = primary[:, :, :, :, :, index, ...]
615
            else:
616
                primary = primary[:, :, :, :, :, index]
617
        elif dimension == 7:
618
            if len(primary.shape) > 7:
619
                primary = primary[:, :, :, :, :, :, index, ...]
620
            else:
621
                primary = primary[:, :, :, :, :, :, index]
622
        elif dimension == 8:
623
            if len(primary.shape) > 8:
624
                primary = primary[:, :, :, :, :, :, :, index, ...]
625
            else:
626
                primary = primary[:, :, :, :, :, :, :, index]
627
        else:
628
            raise RuntimeError("Dimension index '{0}' not supported".format(dimension))
629

    
630
        return primary
631

    
632
## Reduce dimensionality by integrating over a vertical column.
633
#
634
#  This filter produces an integrated vertical column.
635
#
636
#  \f[
637
#  C = \sum_i \frac{10 N_{A} c_i p_i}{g*M}
638
#  \f]
639
#
640
#  with \f$c_i\f$ the vertical profile, \f$p_i\f$ the pressure, both at level \f$i\f$,
641
#  \f$N_A\f$ Avogadro's constant (6.0221367e23), \f$g\f$ the gravitational accelleration
642
#  (9.80665), and \f$M\f$ the molar mass of air (28.9644).
643
#
644
#  The pressure levels \f$p_i\f$ use the surface pressure and two constant arrays
645
#  \f$A_i\f$ and \f$B_i\f$:
646
#
647
#  \f[
648
#  p_i = \frac{A_i + p_s B_i}{100}
649
#  \f]
650
#
651
#  The result of this filter is \f$C\f$.
652
#
653
class IntegratedColumn(AbstractTransformer):
654
    ## The constructor
655
    #
656
    # @param[in] dimension the dimension (index) over which to integrate. Default = -1 (last dimension).
657
    # @param[in] coefficients_a The constant array \f$A_i\f$.
658
    def __init__(self, *args, **kwargs):
659
        super().__init__()
660

    
661
        if 'dimension' in kwargs:
662
            ## The dimension over which to integrate
663
            self.dimension = kwargs['dimension']
664
        else:
665
            self.dimension = None
666

    
667
        if 'coefficients_a' in kwargs:
668
            ## The A constant (array) for the pressure levels in the profile. After initialization this is a string, a call to the read_additional_variables() method reads the actual data.
669
            self.tm5_constant_a = kwargs['coefficients_a']
670
        else:
671
            raise CAMAException("Required argument 'coefficients_a' for the IntegratedColumn filter is missing.")
672

    
673
        if 'coefficients_b' in kwargs:
674
            ## The B constant (array) for the pressure levels in the profile. After initialization this is a string, a call to the read_additional_variables() method reads the actual data.
675
            self.tm5_constant_b = kwargs['coefficients_b']
676
        else:
677
            raise CAMAException("Required argument 'coefficients_b' for the IntegratedColumn filter is missing.")
678

    
679
    ## Read additional parameters from an input file
680
    #
681
    #  Some parameters are only known at runtime, or must be read from file
682
    #  (other than a normal variable).
683
    #
684
    #  Overridden here because we need access to the file to read `tm5_constant_a` and `tm5_constant_b`.
685
    #
686
    # @param f The pycama.File.File object.
687
    #
688
    def read_additional_variables(self, f):
689
        """
690
        Read TM5 a and b constants from file, and calculate the values at the interface.
691
        """
692
        if hasattr(self.tm5_constant_a, 'startswith') and self.tm5_constant_a.startswith('read_from_file("'):
693
            varname = self.tm5_constant_a.replace('read_from_file("', '').replace('")', '')
694
            self.logger.debug("Loading additional variable '{0}'".format(varname))
695
            var = f.find_variable(varname)
696
            if var is None:
697
                raise CAMAException("Skipping '%s' because the variable could not be found in %s.", self.tm5_constant_a[16:-2], f.product)
698
            self.tm5_constant_a = var[0, :]
699
            self.tm5_constant_ai = [0.0]
700
            for fa in self.tm5_constant_a:
701
                self.tm5_constant_ai.append(2*fa - self.tm5_constant_ai[-1])
702
            self.tm5_constant_ai[-1] = 0.0
703
            self.tm5_constant_ai = np.asarray(self.tm5_constant_ai)
704
            self.logger.debug("hyai = [{0}]".format(", ".join([str(a) for a in self.tm5_constant_ai])))
705

    
706
        if hasattr(self.tm5_constant_b, 'startswith') and self.tm5_constant_b.startswith('read_from_file("'):
707
            varname = self.tm5_constant_b.replace('read_from_file("', '').replace('")', '')
708
            self.logger.debug("Loading additional variable '{0}'".format(varname))
709
            var = f.find_variable(varname)
710
            if var is None:
711
                raise CAMAException("Skipping '%s' because the variable could not be found in %s.", self.tm5_constant_b[16:-2], f.product)
712
            self.tm5_constant_b = var[0, :]
713
            self.tm5_constant_bi = [1.0]
714
            for fb in self.tm5_constant_b:
715
                self.tm5_constant_bi.append(2*fb - self.tm5_constant_bi[-1])
716
                if self.tm5_constant_bi[-1] < 1e-7: self.tm5_constant_bi[-1] = 0.0
717
            self.tm5_constant_bi = np.asarray(self.tm5_constant_bi)
718
            self.logger.debug("hybi = [{0}]".format(", ".join([str(b) for b in self.tm5_constant_bi])))
719

    
720
    ## The callable interface for the integration filter.
721
    #
722
    # Note that the surface pressure _must_ be supplied in the `secondary` parameter.
723
    # Otherwise the arguments are as defined in the superclass.
724
    #
725
    # @returns the integrated vertical column in mol/2.
726
    #
727
    def __call__(self,
728
                 primary,
729
                 secondary,
730
                 latitude,
731
                 longitude,
732
                 solar_zenith_angle,
733
                 viewing_zenith_angle,
734
                 solar_azimuth_angle,
735
                 viewing_azimuth_angle,
736
                 processing_quality_flags,
737
                 geolocation_flags,
738
                 surface_classification):
739

    
740
        if self.dimension is not None:
741
            dimension = self.dimension
742
        else:
743
            dimension = -1
744

    
745
        vcd = np.zeros(primary.shape, dtype=np.float32)
746
        dryair = 28.9644
747
        Ra = 6.0221367e23
748
        g = 9.80665
749

    
750
        primary[primary.mask] = np.nan
751
        p = lambda l : (self.tm5_constant_ai[l] + secondary*self.tm5_constant_bi[l])/100.0
752
        for level in range(primary.shape[dimension]):
753
            dp = p(level) - p(level+1)
754
            vcd[:, :, level] = (10.0*Ra*primary[:, :, level]*dp)/(g*dryair)
755

    
756
        primary = np.sum(vcd, axis=dimension)
757
        primary = np.ma.asarray(primary)
758
        primary.mask = np.logical_not(np.isfinite(primary))
759
        return primary*1.6605390404271641565e-20
760

    
761
## Translate a sun normalized radiance into a reflectance
762
#
763
#  Or do the reverse. The conversion (either way) uses solar zenith angles up to 89 degrees.
764
#  Larger angles are clipped.
765
#
766
class SunNormalizedRadiance(AbstractTransformer):
767
    ## The constructor
768
    #
769
    #  @param[in] inverse Do the inverse conversion (from reflectance to sun
770
    #         normalized radiance) when true.
771
    def __init__(self, **kwargs):
772
        super().__init__()
773
        ## Flag to indicate the direction of the translation.
774
        self.inverse = 'inverse' in kwargs and kwargs['inverse']
775

    
776
    ## The callable interface for the SunNormalizedRadiance translator.
777
    #
778
    # The arguments are as defined in the superclass.
779
    #
780
    def __call__(self,
781
                 primary,
782
                 secondary,
783
                 latitude,
784
                 longitude,
785
                 solar_zenith_angle,
786
                 viewing_zenith_angle,
787
                 solar_azimuth_angle,
788
                 viewing_azimuth_angle,
789
                 processing_quality_flags,
790
                 geolocation_flags,
791
                 surface_classification):
792

    
793
        limit = 88.5
794

    
795
        if self.inverse:
796
            try:
797
                mu0 = np.cos(np.where(solar_zenith_angle < limit, solar_zenith_angle, limit)*np.pi/180.0)
798
            except TypeError:
799
                mu0 = np.pi
800

    
801
            return primary*mu0/np.pi
802
        else:
803
            try:
804
                mu0 = np.cos(np.where(solar_zenith_angle < limit, solar_zenith_angle, limit)*np.pi/180.0)
805
            except TypeError:
806
                mu0 = np.pi
807

    
808
            return np.pi*primary/mu0
809

    
810

    
811
## Snow/Ice filter
812
class SnowIceFilter(AbstractTransformer):
813
    ## The constructor
814
    #
815
    #  @param filter What to remove ("snow", "no_snow", default "snow")
816
    #  @param threshold Threshold for the ice fraction before removal (0.2)
817
    def __init__(self, **kwargs):
818
        super().__init__()
819

    
820
        if 'filter' in kwargs:
821
            self.filter = kwargs['filter']
822
        else:
823
            self.filter = 'snow'
824

    
825
        if self.filter not in ('snow', 'no_snow', 'none'):
826
            self.logger.warning("Filter must be one of ('snow', 'no_snow', 'none'), resetting to 'snow'")
827
            self.filter = 'snow'
828

    
829
        if 'threshold' in kwargs:
830
            self.threshold = int(100*kwargs['threshold'])
831
        else:
832
            self.threshold = 20
833

    
834
        self.OMI_shift = 'OMI_shift' in kwargs and kwargs['OMI_shift']
835

    
836
    ## The callable
837
    #
838
    # The arguments are as defined in the superclass.
839
    def __call__(self,
840
                 primary,
841
                 secondary,
842
                 latitude,
843
                 longitude,
844
                 solar_zenith_angle,
845
                 viewing_zenith_angle,
846
                 solar_azimuth_angle,
847
                 viewing_azimuth_angle,
848
                 processing_quality_flags,
849
                 geolocation_flags,
850
                 surface_classification):
851

    
852
        if self.filter == 'none':
853
            return primary
854

    
855
        try:
856
            mask = primary.mask
857
            dummy = mask.shape[1]
858
        except (AttributeError, IndexError):
859
            primary = np.ma.asarray(primary)
860
            mask = np.zeros(primary.shape, dtype=np.bool)
861

    
862
        total = np.sum(mask)
863

    
864
        if self.OMI_shift:
865
            primary = (primary >> 8) & 0x7F
866
        else:
867
            primary = primary & 0x7F
868

    
869
        if self.filter == 'snow':
870
            mask = np.logical_or(mask, np.logical_and(primary > self.threshold, primary <= 103))
871
            new_total = np.sum(mask)
872
            primary.mask = mask
873
            self.logger.debug("{0}: Removing {1} pixels for snow/ice".format(self.__class__.__name__,
874
                              new_total-total))
875
        else:
876
            mask = np.logical_or(mask, primary < self.threshold)
877
            new_total = np.sum(mask)
878
            primary.mask = mask
879
            self.logger.debug("{0}: Removing {1} pixels for no snow/ice".format(self.__class__.__name__, new_total-total))
880

    
881

    
882
        return primary
883

    
884
## Resize the primary array to a new width of the swath.
885
#
886
class Resize(AbstractTransformer):
887
    ## The constructor
888
    #
889
    #  @param shape The new width (defaults to None, which will reuse the
890
    #               width of the latitude variable).
891
    #  @param skip  The stride in which to sample (default: 1, sample all).
892
    #  @param dimension  The dimension to reduce (default: 1, the `ground_pixel` dimension).
893
    #                    Can be up to 2.
894
    def __init__(self, **kwargs):
895
        super().__init__()
896

    
897
        if 'shape' in kwargs:
898
            ## The new number of ground pixels (rows).
899
            self.shape = kwargs['shape']
900
        else:
901
            self.shape = None
902

    
903
        if 'skip' in kwargs:
904
            ## The stride to use when filtering
905
            self.skip = kwargs['skip']
906
        else:
907
            self.skip = 1
908

    
909
        if 'dimension' in kwargs:
910
            ## The dimension to reduce.
911
            self.dimension = kwargs['dimension']
912
        else:
913
            self.dimension = 1
914

    
915
    ## The callable
916
    #
917
    # The arguments are as defined in the superclass.
918
    def __call__(self,
919
                 primary,
920
                 secondary,
921
                 latitude,
922
                 longitude,
923
                 solar_zenith_angle,
924
                 viewing_zenith_angle,
925
                 solar_azimuth_angle,
926
                 viewing_azimuth_angle,
927
                 processing_quality_flags,
928
                 geolocation_flags,
929
                 surface_classification):
930

    
931
        if self.shape is not None:
932
            shape = self.shape
933
        else:
934
            try:
935
                shape = latitude.shape
936
            except AttributeError:
937
                shape = 60
938

    
939
        if self.skip is not None:
940
            skip = self.skip
941
        else:
942
            skip = 1
943

    
944
        if self.dimension is not None:
945
            dimension = self.dimension
946
        else:
947
            dimension = 1
948

    
949
        if dimension == 0:
950
            if len(primary.shape) == 1:
951
                return np.repeat(primary.reshape((1, primary.shape[0])), shape, axis=0)
952
            elif primary.shape[0] == 1:
953
                return np.repeat(primary, shape, axis=0)
954
            else:
955
                return primary[0:shape:skip, ...]
956
        elif dimension == 1:
957
            if len(primary.shape) == 1:
958
                return np.repeat(primary.reshape((primary.shape[0], 1)), shape, axis=1)
959
            elif len(primary.shape) == 2:
960
                if primary.shape[1] == 1:
961
                    return np.repeat(primary, shape, axis=1)
962
                else:
963
                    return primary[:, 0:shape:skip]
964
            else:
965
                if primary.shape[1] == 1:
966
                    return np.repeat(primary, shape, axis=1)
967
                else:
968
                    return primary[:, 0:shape:skip, ...]
969
        elif dimension == 2:
970
            if len(primary.shape) < 2:
971
                raise CAMAException("{0}: Source data has not enough dimensions".format(self.__class__.__name__))
972
            elif len(primary.shape) == 2:
973
                return np.repeat(primary.reshape((primary.shape[0], primary.shape[1], 1)), shape, axis=2)
974
            elif len(primary.shape) == 3:
975
                if primary.shape[2] == 1:
976
                    return np.repeat(primary, shape, axis=2)
977
                else:
978
                    return primary[:, :, 0:shape:skip]
979
            else:
980
                if primary.shape[2] == 1:
981
                    return np.repeat(primary, shape, axis=2)
982
                else:
983
                    return primary[:, :, 0:shape:skip, ...]
984

    
985

    
986
## DLR pixel mapper
987
#
988
#  The cloud product uses both band 3, 4 and band 6.
989
#  This filter maps from band 3 to band 6 or vice versa, using a direct pixel
990
#  mapping scheme. A full colocation filter for PyCAMA is still on the to-do list.
991
#
992
class CloudProductPixelMapper(AbstractTransformer):
993
    ## The constructor
994
    #
995
    #  @param[in] target Target band for the mapping. Defaults to 6.
996
    def __init__(self, **kwargs):
997
        super().__init__()
998
        ## Flag to indicate the direction of the translation.
999
        if 'target' in kwargs and kwargs['target'] and kwargs['target'] in (3, 4, 5, 6):
1000
            self.target = 6 if kwargs['target'] in (5, 6) else 3
1001
        else:
1002
            self.target = 6
1003

    
1004
    ## The callable interface for the SunNormalizedRadiance translator.
1005
    #
1006
    # The arguments are as defined in the superclass.
1007
    #
1008
    def __call__(self,
1009
                 primary,
1010
                 secondary,
1011
                 latitude,
1012
                 longitude,
1013
                 solar_zenith_angle,
1014
                 viewing_zenith_angle,
1015
                 solar_azimuth_angle,
1016
                 viewing_azimuth_angle,
1017
                 processing_quality_flags,
1018
                 geolocation_flags,
1019
                 surface_classification):
1020

    
1021
        if self.target == 6:
1022
            # cut of slice
1023
            if primary is not None and primary.shape[0] == 450:
1024
                primary = primary[2:, ...]
1025
            if secondary is not None and secondary.shape[0] == 450:
1026
                secondary = secondary[2:, ...]
1027
            if latitude is not None and latitude.shape[0] == 450:
1028
                latitude = latitude[2:, ...]
1029
            if longitude is not None and longitude.shape[0] == 450:
1030
                longitude = longitude[2:, ...]
1031
            if solar_zenith_angle is not None and solar_zenith_angle.shape[0] == 450:
1032
                solar_zenith_angle = solar_zenith_angle[2:, ...]
1033
            if viewing_zenith_angle is not None and viewing_zenith_angle.shape[0] == 450:
1034
                viewing_zenith_angle = viewing_zenith_angle[2:, ...]
1035
            if solar_azimuth_angle is not None and solar_azimuth_angle.shape[0] == 450:
1036
                solar_azimuth_angle = solar_azimuth_angle[2:, ...]
1037
            if viewing_azimuth_angle is not None and viewing_azimuth_angle.shape[0] == 450:
1038
                viewing_azimuth_angle = viewing_azimuth_angle[2:, ...]
1039
            if processing_quality_flags is not None and processing_quality_flags.shape[0] == 450:
1040
                processing_quality_flags = processing_quality_flags[2:, ...]
1041
            if geolocation_flags is not None and geolocation_flags.shape[0] == 450:
1042
                geolocation_flags = geolocation_flags[2:, ...]
1043
            if surface_classification is not None and surface_classification.shape[0] == 450:
1044
                surface_classification = surface_classification[2:, ...]
1045
        elif self.target == 3:
1046
            # add slice
1047
            b = np.zeros(tuple([n if (i != 0) else 2 for i, n in enumerate(primary.shape)]), dtype=np.float64) + np.nan
1048
            if primary is not None and primary.shape[0] == 448:
1049
                primary = np.concatenate((b, primary))
1050
            if secondary is not None and secondary.shape[0] == 448:
1051
                secondary = np.concatenate((b, secondary))
1052
            if latitude is not None and latitude.shape[0] == 448:
1053
                latitude = np.concatenate((b, latitude))
1054
            if longitude is not None and longitude.shape[0] == 448:
1055
                longitude = np.concatenate((b, longitude))
1056
            if solar_zenith_angle is not None and solar_zenith_angle.shape[0] == 448:
1057
                solar_zenith_angle = np.concatenate((b, solar_zenith_angle))
1058
            if viewing_zenith_angle is not None and viewing_zenith_angle.shape[0] == 448:
1059
                viewing_zenith_angle = np.concatenate((b, viewing_zenith_angle))
1060
            if solar_azimuth_angle is not None and solar_azimuth_angle.shape[0] == 448:
1061
                solar_azimuth_angle = np.concatenate((b, solar_azimuth_angle))
1062
            if viewing_azimuth_angle is not None and viewing_azimuth_angle.shape[0] == 448:
1063
                viewing_azimuth_angle = np.concatenate((b, viewing_azimuth_angle))
1064
            if processing_quality_flags is not None and processing_quality_flags.shape[0] == 448:
1065
                processing_quality_flags = np.concatenate((b, processing_quality_flags))
1066
            if geolocation_flags is not None and geolocation_flags.shape[0] == 448:
1067
                geolocation_flags = np.concatenate((b, geolocation_flags))
1068
            if surface_classification is not None and surface_classification.shape[0] == 448:
1069
                surface_classification = np.concatenate((b, surface_classification))
1070

    
1071
        return primary