Project

General

Profile

Revision 817:b1891bf5727b

IDb1891bf5727b
Parent 816:d52976b6ccfc
Child 818:8316d57368b6

Added by Maarten Sneep about 2 years ago

New script to check irradiance files for corruption

View differences:

src/check_irradiance.py
1
#!/usr/bin/env python3
2

  
3
import os
4
import sys
5
import warnings
6
import datetime, time
7
import warnings
8
import logging
9
import logging.config
10
import argparse
11
import socket
12

  
13
import netCDF4
14
import numpy as np
15

  
16

  
17
__version__ = "1.0.0"
18

  
19
## Extra logging level for progress messages.
20
#
21
#  Injected at the same level as expected by the thin layer interface.
22
PROGRESS_LEVEL_NUM = 35
23

  
24
## Create a new logging level with the name "PROGRESS"
25
#
26
# This code is executed at load time.
27
# \code
28
# logging.addLevelName(PROGRESS_LEVEL_NUM, "PROGRESS")
29
# \endcode
30
#
31
# The progress_method() is the actual logger, similar to the \c warning and \c debug methods of the logger.
32
# The function is injected as the method \c progress into the <tt>logging.Logger</tt> class.
33
#
34
# \code
35
# logging.Logger.progress = progress_method
36
# \endcode
37
#
38
# we're not calling the function here, just inserting the method into the class object.
39
# Gotta love python for its monkey-patching.
40
logging.addLevelName(PROGRESS_LEVEL_NUM, "PROGRESS")
41

  
42
def progress_method(self, message, *args, **kws):
43
    # Yes, logger takes its '*args' as 'args'.
44
    if self.isEnabledFor(PROGRESS_LEVEL_NUM):
45
        self._log(PROGRESS_LEVEL_NUM, message, args, **kws)
46

  
47
logging.Logger.progress = progress_method
48

  
49
## A formatter class specifically for the production logger.
50
#
51
#  The main aim for this subclass of `logging.Formatter` is to format the timestamp
52
#  accoding to the thin layer interface, i.e. including microseconds. This is not
53
#  supported out of the box by the `logging` framework of Python 3.
54
#
55
#  The base class is documented
56
#  [on the Python 3 documentation site](https://docs.python.org/3/library/logging.html#formatter-objects).
57
#
58
class MyFormatter(logging.Formatter):
59
    ## Convert the creation date for easier formatting
60
    #
61
    #  @param[in] r Creation timestamp in unix epoch-seconds.
62
    #  @return  `datetime` object in the `utc` timezone.
63
    def converter(self, r):
64
        return datetime.datetime.fromtimestamp(r, tz=datetime.timezone.utc)
65

  
66
    ## Format the time into an ISO 8601 string
67
    #
68
    # Override the [Python standard library `logging.Formatter.formatTime` method to include microseconds
69
    # @param[in] record a `logging.LogRecord` object. Only the `record.created` attribute is used.
70
    # @param[in] datefmt A specific format to use for the date formatting
71
    # @return a formatted string.
72
    def formatTime(self, record, datefmt=None):
73
        ct = self.converter(record.created)
74
        if datefmt:
75
            s = ct.strftime(datefmt)
76
        else:
77
            s = ct.strftime("%Y-%m-%dT%H:%M:%S.%f")
78
        return s
79

  
80
## A filter object for the logging facility
81
#
82
#  [See the python documentation for a full specification](https://docs.python.org/3/library/logging.html#filter-objects).
83
class ErrorLogFilter(object):
84
    ## Initialization
85
    #
86
    #  set up the filter
87
    #  @param[in] outlevel Level for standard output logging
88
    #  @param[in] errlevel Level for standard error logging
89
    def __init__(self, outlevel, errlevel):
90
        self.outlevel = outlevel if outlevel >= logging.DEBUG else logging.DEBUG
91
        self.errlevel = errlevel if errlevel <= logging.ERROR else logging.ERROR
92
    ## Is the specified record to be logged?
93
    #
94
    #  @param[in] record The record to operate on
95
    #  @returns 0 for no, non-zero for yes.
96
    #
97
    #  As all logging records pass through this method, we also count the number of errors and warnings here.
98
    def filter(self, record):
99
        return self.outlevel <= record.levelno < self.errlevel
100

  
101
def setup_logging(production_logger=True, loglevel="warning", errlevel="error"):
102
    global __version__
103
    warnings.simplefilter("ignore")
104

  
105
    try:
106
        num_levels = [int(getattr(logging, loglevel.upper(), None)), int(getattr(logging, errlevel.upper(), None))]
107
    except TypeError:
108
        raise ValueError("Invalid log level(s) in configuration: {0}, {1}.\nValid levels are 'DEBUG', 'INFO', 'WARNING', 'PROGRESS' or 'ERROR'".format(loglevel.upper(), errlevel.upper()))
109

  
110
    for i in range(2):
111
        if num_levels[i] > logging.ERROR:
112
            num_levels[i] = logging.ERROR
113

  
114
    if num_levels[0] == num_levels[1]:
115
        if num_levels[0] == logging.ERROR:
116
            num_levels = (logging.WARNING, logging.ERROR)
117
        else:
118
            num_levels = (num_levels[0], logging.ERROR)
119

  
120
    # create logger
121
    logger = logging.getLogger(os.path.basename(sys.argv[0]))
122

  
123
    if len(logger.handlers) == 0:
124
        # create formatter
125
        # 2016-06-14T09:44:15.482164 bhltrdl2 TROPNLL2DP 0.10.1 [0000006984]: [I] Starting TROPNLL2DP version 0.10.1
126
        # 2016-06-14T12:13:52.396889 bhw465.knmi.nl PyCAMA 0.1.0 [0000020929]: [I] Opening file 'G2A_TEST_L2__AER_AI_20100319T112929_20100319T131047_17713_02_001000_20160502T150326.nc'
127
        if production_logger:
128
            fmt='%(asctime)s {machine} %(name)s {version} [%(process)010d]: [%(levelname)-.1s] %(message)s'
129
            formatter = MyFormatter(fmt=fmt.format(version=__version__, machine=socket.gethostname()),
130
                                    datefmt='%Y-%m-%dT%H:%M:%S.%f')
131
        else:
132
            formatter = logging.Formatter(fmt='[%(levelname)8s] %(filename)s [%(lineno)4d]: %(message)s',
133
                                          datefmt='%Y-%m-%dT%H:%M:%S')
134

  
135
        logger.setLevel(min(num_levels))
136
        for output_channel in (sys.stdout, sys.stderr):
137
            # create console handler and set level to
138
            ch = logging.StreamHandler(stream=output_channel)
139

  
140
            if output_channel is sys.stdout:
141
                ch.setLevel(logging.DEBUG)
142
                # add filter to ignore error messages
143
                ch.addFilter(ErrorLogFilter(num_levels[0], num_levels[1]))
144
            else:
145
                ch.setLevel(num_levels[1])
146

  
147
            # add formatter to ch
148
            ch.setFormatter(formatter)
149

  
150
            # add ch to logger
151
            logger.addHandler(ch)
152

  
153
    return logger
154

  
155
## simplify a list of indices for reporting purposes.
156
#
157
# Reduce a list of indices to a sequence of continuous ranges, making long ranges more human readable.
158
# @param r The list to simplify.
159
# @returns A string for printing.
160
#
161
# Example:
162
# \code
163
# r = [0, 1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 14, 16]
164
# print(row_list_simplification(r))
165
# \endcode
166
# prints `0-5, 7-9, 11, 13-14, 16`
167
def row_list_simplification(r):
168
    rows = sorted(r)
169
    groups = []
170
    i = 1
171
    while i < len(rows):
172
        group = [rows[i-1]]
173
        while i < len(rows) and rows[i] <= rows[i-1]+1:
174
            i+=1
175
        group.append(rows[i-1])
176
        groups.append(group)
177
        i+=1
178
    if groups[-1][-1] != rows[-1]:
179
        groups.append([rows[-1], rows[-1]])
180

  
181
    s = []
182
    for group in groups:
183
        if group[0] != group[1]:
184
            s.append("{0[0]:d}-{0[1]:d}".format(group))
185
        else:
186
            s.append("{0[0]:d}".format(group))
187
    return ", ".join(s)
188

  
189
def cmdline():
190
    global __version__
191
    descr = """Check the quality of a set of irradiance files. 
192
An exit status of 0 indicates a good irradiance set, 127 indicates a suspect dataset, and 255 indicates a bad dataset. 
193
In case of a suspect or bad dataset, a pdf file with plots will be produced in the indicated directory, with a file name based on the input data. 
194
This file should be attached to a new Redmine issue for manual inspection by MPC (Maarten Sneep, Mark ter Linden, Antje Ludewig, Pascal Hedelt, Fabian Rohman). Please set the priority of the issue to urgent, assign to Maarten, and add the others as watchers.
195
For suspect of bad data-sets please prevent usage in any production stream (NRTI, OFFL). """
196
    parser = argparse.ArgumentParser(description=descr,
197
                                     prog=os.path.basename(sys.argv[0]))
198
    parser.add_argument('--version', action='version',
199
                        version='%(prog)s {0}'.format(__version__))
200
    parser.add_argument('-l', '--log', metavar='LEVEL', 
201
                        choices=('debug', 'info', 'progress', 'warning', 'error'),
202
                        help="Set the log level ('debug', 'info', 'progress', 'warning', 'error')", default="info")
203
    # input files
204
    parser.add_argument('-u', '--uvn', default=None, required=True,
205
                        help='UVN irradiance file to check (L1B_IR_UVN)')
206
    parser.add_argument('-s', '--swir', default=None, required=True,
207
                        help='SWIR irradiance file to check (L1B_IR_SIR)')
208
    # output options
209
    parser.add_argument('-o', "--out", metavar="DIR", help="Output directory", default=os.getcwd())
210
    parser.add_argument('--keep', action='store_true', help='Keep pdf file, regardless of test result')
211
    parser.add_argument('--no-fig', action='store_true', help="No not make pdf output")
212
    args = parser.parse_args()
213

  
214
    return args, parser
215

  
216

  
217
class dummy_with(object):
218
    def __init__(self):
219
        pass
220
    def __enter__(self):
221
        pass
222
    def __exit__(self, t, value, tb):
223
        pass
224
    def savefig(self):
225
        pass
226
    def infodict(self):
227
        return {}
228

  
229
class Weight: 
230
    def __init__(self, startwvl, endwvl, rampup): 
231
        self.startwvl = startwvl 
232
        self.endwvl = endwvl 
233
        self.delta = rampup 
234
    def __call__(self, wvl): 
235
        result = np.zeros(wvl.shape, dtype=np.float64) 
236
        for i, w in enumerate(wvl): 
237
            if w < self.startwvl or w > self.endwvl: 
238
                result[i] = 0.0 
239
            elif self.startwvl < w < self.startwvl + self.delta: 
240
                result[i] = (w-self.startwvl)/self.delta 
241
            elif self.endwvl-self.delta < w < self.endwvl: 
242
                result[i] = (self.endwvl-w)/self.delta 
243
            else: 
244
                result[i] = 1.0 
245
        return result 
246

  
247
def integrate_irrad(ref, band):
248
    irrad = ref[f'BAND{band:1d}_IRRADIANCE/STANDARD_MODE/OBSERVATIONS/irradiance'][0, 0, :, :]
249
    wvl = ref[f'BAND{band:1d}_IRRADIANCE/STANDARD_MODE/INSTRUMENT/calibrated_wavelength'][0, :, :]
250
    weight = Weight(max(wvl.min(axis=1))+1, min(wvl.max(axis=1))-1, 4.0)
251
    
252
    result = np.zeros((wvl.shape[0],), dtype=np.float64)
253
    for i in range(wvl.shape[0]):
254
        result[i] = np.sum(weight(wvl[i, 0:-1])*irrad[i, 0:-1]/(wvl[i, 1:] - wvl[i, 0:-1]))
255
    return result
256

  
257

  
258
def main():
259
    global __version__
260
    
261
    args, parser = cmdline()
262
    logger = setup_logging(production_logger=True, loglevel=args.log, errlevel="error")
263

  
264
    refuvn = netCDF4.Dataset(args.uvn)
265
    refsir = netCDF4.Dataset(args.swir)
266
    
267
    orbit = refuvn.orbit
268
    cov_start = datetime.datetime.strptime(refuvn.time_coverage_start, "%Y-%m-%dT%H:%M:%SZ")
269
    collection = os.path.basename(args.uvn)[58:60]
270
    
271
    vstring = "".join(["{0:02d}".format(int(s)) for s in __version__.split('.')])
272
    outfname = os.path.basename(args.uvn)[:52].replace('L1B_IR_UVN', 'MPC_ICHECK') + f"{collection}_{orbit:05d}_{vstring}_{datetime.datetime.utcnow():%Y%m%dT%H%M%S}.pdf"
273
        
274
    results = {}
275
    thresholds = [[0.01, 0.03], # 1
276
                  [0.06, 0.08], # 2
277
                  [0.02, 0.05], # 3
278
                  [0.01, 0.03], # 4
279
                  [0.01, 0.03], # 5
280
                  [0.01, 0.03], # 6
281
                  [0.025, 0.04], # 7
282
                  [0.01, 0.03]] # 8
283
    
284
    logger.info(f"Checking irradiance for orbit {orbit}")
285
    
286
    with refuvn, refsir:
287
        for band in range(1, 9):
288
            logger.progress(f"Processing band {band}")
289
            if 1 <= band <= 6:
290
                ref = refuvn
291
            else:
292
                ref = refsir
293
            grp = ref[f'BAND{band:1d}_IRRADIANCE/STANDARD_MODE/OBSERVATIONS']
294
            integrated = integrate_irrad(ref, band)
295
            d_integrated = 2*(integrated[1:] - integrated[0:-1])/(integrated[1:] + integrated[0:-1])
296
            max_deviation = np.max(np.fabs(d_integrated[1:-1]))
297
            irrad = grp['irradiance'][0, 0, :, :]
298
            results[f'band_{band:1d}'] = {'integrated': integrated, 
299
                                          'd_integrated': d_integrated,
300
                                          'max_deviation': max_deviation,
301
                                          'irrad': irrad}
302
    have_warning = False
303
    have_error = False
304
    bands_with_issues = []
305
    
306
    for band in range(1, 9):
307
        max_dev = results[f'band_{band:1d}']['max_deviation']
308
        
309
        if max_dev < thresholds[band-1][0]:
310
            logger.info(f"Max deviation for band {band} is {max_dev:.4f}")
311
        elif thresholds[band-1][0] <= max_dev < thresholds[band-1][1]:
312
            logger.warning(f"Max deviation for band {band} is {max_dev:.4f}")
313
            have_warning = True
314
            bands_with_issues.append(band)
315
        else:
316
            logger.error(f"Max deviation for band {band} is {max_dev:.4f}")
317
            have_error = True
318
            bands_with_issues.append(band)
319
        
320

  
321
    if (not args.no_fig) and (have_warning or have_error or args.keep):
322
        logger.info("Creating graphical output")
323
        from matplotlib import pyplot as plt
324
        from matplotlib.backends.backend_pdf import PdfPages
325
        import matplotlib as mpl
326

  
327
        mpl.use('pdf')
328
        pdf = PdfPages(os.path.join(args.out, outfname), keep_empty=False)
329
        with pdf:
330
            for band in range(1,9):
331
                # add figures
332
                logger.debug(f"Plotting integrated irradiance for band {band}")
333
                fig = plt.figure()
334
                rows = np.arange(results[f'band_{band:1d}']['integrated'].shape[0])
335
                plt.plot(rows, results[f'band_{band:1d}']['integrated'], 'k-', label=f'Orbit {orbit}, band {band}')
336
                plt.xlabel('Row')
337
                plt.ylabel('Integrated irradiance')
338
                plt.legend()
339
                plt.tight_layout()
340
                pdf.savefig()
341
                plt.close()
342
            
343
                logger.debug(f"Plotting row difference of integrated irradiance for band {band}")
344
                fig = plt.figure()
345
                ax = plt.gca()
346
                plt.plot(rows[:-1]+0.5, results[f'band_{band:1d}']['d_integrated'], 'k-', label=f'Orbit {orbit}, band {band}')
347
                plt.xlabel('Row')
348
                plt.ylabel('Δ Integrated irradiance')
349
            
350
                plt.text(0.05, 0.15, f"Maximum relative row-to-row\nvariation {results[f'band_{band:1d}']['max_deviation']:.4f}", 
351
                    transform=ax.transAxes, horizontalalignment='left', verticalalignment='center', 
352
                    fontsize=8, color=('r' if results[f'band_{band:1d}']['max_deviation'] > thresholds[band-1][0] else 'k'))
353
                plt.legend()
354
                plt.tight_layout()
355
                pdf.savefig()
356
                plt.close()
357
            
358
                logger.debug(f"Plotting irradiance for band {band}")
359
                fig = plt.figure()
360
                ax = plt.gca()
361
                vmin, vmax = np.percentile(results[f'band_{band:1d}']['irrad'], [2, 98])
362
                norm = mpl.colors.Normalize(vmin, vmax, clip=True)
363
                cmap = mpl.cm.viridis
364
                im = ax.imshow(results[f'band_{band:1d}']['irrad'], cmap=cmap, norm=norm, 
365
                    aspect='auto', interpolation='antialiased', origin='lower')
366
                plt.xlabel('Column')
367
                plt.ylabel('Row')
368
                cbar = fig.colorbar(im, ax=ax)
369
                cbar.ax.set_ylabel("Irradiance [mol·nm⁻¹·m⁻²·s⁻¹]", rotation=-90, va="bottom")
370
                fig.tight_layout()
371
                pdf.savefig()
372
                plt.close()
373
        
374
            metadata = pdf.infodict()
375
            metadata['Title'] = f"Irradiance check report for orbit {orbit:05d} on {cov_start:%Y-%m-%d}"
376
            metadata['Author'] = "Maarten Sneep"
377
            metadata['Subject'] = "Level 2 anomaly prevention"
378
            metadata['Keywords'] = ", ".join(['Sentinel 5 precursor', "Anomaly", "Error prevention",  "Quality control", f"{orbit:05d}"])
379
            metadata['Creator'] = "PDGS"
380
            metadata['Producer'] =  os.path.basename(sys.argv[0])
381
            metadata['CreationDate'] = f"{datetime.datetime.utcnow():%Y-%m-%dT%H:%M:%S}"
382
            metadata['ModDate'] = None
383
    
384
    if len(bands_with_issues) > 0:
385
        logger.info("Tests completed, found issues in orbit {0} for band(s) {1},".format(orbit, ", ".join([str(s) for s in bands_with_issues])))
386
    else:
387
        logger.info(f"Tests completed, no issues found in {orbit}.")
388
    
389
    if have_error:
390
        exitval = 255
391
    elif have_warning:
392
        exitval = 127
393
    else:
394
        exitval = 0
395
        if not args.keep:
396
            try:
397
                os.unlink(os.path.join(args.out, outfname))
398
            except:
399
                pass
400
    
401
    if exitval == 255:
402
        logger.error("Terminating with exit value 255")
403
    elif exitval == 127:
404
        logger.warning("Terminating with exit value 127")
405
    else:
406
        logger.info("Terminating with exit value 0")
407
    sys.exit(exitval)
408
    
409
if __name__ == "__main__":
410
    main()

Also available in: Unified diff