Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / plot_s5p_spectrum.py @ 808:9dda39dd51cb

History | View | Annotate | Download (16 KB)

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

    
5
# PyCAMA Copyright (c) 2016, KNMI (Maarten Sneep).
6
#
7
# Permission is hereby granted, free of charge, to any person obtaining a
8
# copy of this software and associated documentation files (the "Software"),
9
# to deal in the Software without restriction, including without limitation
10
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
11
# and/or sell copies of the Software, and to permit persons to whom the
12
# Software is furnished to do so, subject to the following conditions:
13
#
14
# - The above copyright notice and this permission notice shall be included
15
#   in all copies or substantial portions of the Software.
16
#
17
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
18
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
20
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
21
# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
22
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
23
# OTHER DEALINGS IN THE SOFTWARE.
24

    
25
import os, sys, argparse
26
import matplotlib.pyplot as plt
27
import datetime
28
import warnings
29

    
30
import netCDF4
31
import numpy as np
32
from scipy.interpolate import InterpolatedUnivariateSpline as Spline
33

    
34
warnings.filterwarnings("ignore", category=UserWarning)
35

    
36
def nc4_path(ref, p):
37
    grp = ref
38
    for gname in p.split('/')[1:]:
39
        grp = grp.groups[gname]
40
    return grp
41

    
42
def read_var_filter(var, pixel, scanline=None):
43
    var.set_auto_maskandscale(False)
44
    try:
45
        fv = var._FilleValue
46
    except:
47
        fv = float.fromhex('1.ep122')
48
    if scanline is None:
49
        if len(var.shape) == 3:
50
            data = var[0, pixel, :]
51
            data = np.where(data == fv, np.nan, data)
52
        elif len(var.shape) == 2:
53
            data = var[0, pixel]
54
    else:
55
        if len(var.shape) == 4:
56
            data = var[0, scanline, pixel, :]
57
            data = np.where(data == fv, np.nan, data)
58
        elif len(var.shape) == 3:
59
            data = var[0, scanline, pixel]
60
    return data
61

    
62
def get_spectrum(fname, band, pixel, scanline=0, radiance=True):
63
    with netCDF4.Dataset(fname, 'r') as ref:
64
        try:
65
            p = '/BAND{0}_{1}/STANDARD_MODE/INSTRUMENT'.format(band, 'RADIANCE' if radiance else 'IRRADIANCE')
66
            grp = nc4_path(ref, p)
67
        except:
68
            raise
69

    
70
        try:
71
            var = grp.variables['calibrated_wavelength']
72
            wvl = read_var_filter(var, pixel)
73
        except:
74
            try:
75
                var = grp.variables['nominal_wavelength']
76
                wvl = read_var_filter(var, pixel)
77
            except:
78
                print("Could not extract wavelengths", file=sys.stderr)
79
                raise
80

    
81
        try:
82
            p = '/BAND{0}_{1}/STANDARD_MODE/OBSERVATIONS'.format(band, 'RADIANCE' if radiance else 'IRRADIANCE')
83
            grp = nc4_path(ref, p)
84
        except:
85
            raise
86

    
87
        try:
88
            var = grp.variables['radiance' if radiance else 'irradiance']
89
            r = read_var_filter(var, pixel, scanline)
90
            var = grp.variables['radiance_error' if radiance else 'irradiance_error']
91
            r_err = read_var_filter(var, pixel, scanline)
92
            var = grp.variables['radiance_noise' if radiance else 'irradiance_noise']
93
            r_noise = read_var_filter(var, pixel, scanline)
94
        except:
95
            print("Could not extract {0}".format('radiances' if radiance else 'irradiances'), file=sys.stderr)
96
            raise
97

    
98
        if radiance:
99
            try:
100
                p = '/BAND{0}_{1}/STANDARD_MODE/GEODATA'.format(band, 'RADIANCE' if radiance else 'IRRADIANCE')
101
                grp = nc4_path(ref, p)
102
            except:
103
                raise
104

    
105
            try:
106
                sza = grp.variables['solar_zenith_angle'][0, scanline, pixel]
107
            except:
108
                print("Could not extract solar zenith angle", file=sys.stderr)
109
                raise
110

    
111
            try:
112
                lat = grp.variables['latitude'][0, scanline, pixel]
113
            except:
114
                print("Could not extract latitude", file=sys.stderr)
115
                raise
116

    
117
            try:
118
                lon = grp.variables['longitude'][0, scanline, pixel]
119
            except:
120
                print("Could not extract longitude", file=sys.stderr)
121
                raise
122

    
123
            try:
124
                p = '/BAND{0}_{1}/STANDARD_MODE/OBSERVATIONS'.format(band, 'RADIANCE' if radiance else 'IRRADIANCE')
125
                grp = nc4_path(ref, p)
126
            except:
127
                raise
128

    
129
            try:
130
                dtime = grp.variables['delta_time']
131
            except:
132
                print("Could not extract delta_time", file=sys.stderr)
133
                raise
134

    
135
            reftime = datetime.datetime.strptime(dtime.units, 'milliseconds since %Y-%m-%d %H:%M:%S')
136
            time = reftime + datetime.timedelta(microseconds=int(1000*dtime[0, scanline]))
137
        else:
138
            sza = None
139
            lat = None
140
            lon = None
141
            time = None
142

    
143
    noise = r/10.0**(r_noise/10.0)
144
    error = r/10.0**(r_err/10.0)
145

    
146
    return wvl, r, noise, error, sza, lat, lon, time
147

    
148
def get_fresco(fresco_file, pixel, scanline):
149
    data = {}
150
    ref = netCDF4.Dataset(fresco_file, 'r')
151
    with ref:
152
        grp = nc4_path(ref, '/PRODUCT')
153
        data['fresco_fraction'] = grp.variables['cloud_fraction_crb'][0, scanline, pixel]
154
        data['fresco_pressure'] = grp.variables['cloud_pressure_crb'][0, scanline, pixel]
155
        data['fresco_height'] = grp.variables['cloud_height_crb'][0, scanline, pixel]
156
        data['fresco_fraction_precision'] = grp.variables['cloud_fraction_crb_precision'][0, scanline, pixel]
157
        data['fresco_pressure_precision'] = grp.variables['cloud_pressure_crb_precision'][0, scanline, pixel]
158
        data['fresco_height_precision'] = grp.variables['cloud_height_crb_precision'][0, scanline, pixel]
159
        data['fresco_scene_albedo'] = grp.variables['scene_albedo'][0, scanline, pixel]
160
        data['fresco_scene_albedo_precision'] = grp.variables['scene_albedo_precision'][0, scanline, pixel]
161
        data['fresco_scene_pressure'] = grp.variables['apparent_scene_pressure'][0, scanline, pixel]
162
        data['fresco_scene_pressure_precision'] = grp.variables['apparent_scene_pressure_precision'][0, scanline, pixel]
163

    
164
        grp = nc4_path(ref, '/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS')
165
        data['chi_square'] = grp.variables['chi_square'][0, scanline, pixel]
166
        data['error_covariance_matrix_element'] = grp.variables['error_covariance_matrix_element'][0, scanline, pixel]
167
        data['number_of_iterations'] = grp.variables['number_of_iterations'][0, scanline, pixel]
168
        data['processing_quality_flags'] = grp.variables['processing_quality_flags'][0, scanline, pixel]
169
        data['number_of_spectral_points_in_retrieval'] = grp.variables['number_of_spectral_points_in_retrieval'][0, scanline, pixel]
170
        data['wavelength_calibration_offset'] = grp.variables['wavelength_calibration_offset'][0, scanline, pixel]
171
        data['wavelength_calibration_irradiance_offset'] = grp.variables['wavelength_calibration_irradiance_offset'][0, pixel]
172

    
173
        grp = nc4_path(ref, '/PRODUCT/SUPPORT_DATA/INPUT_DATA')
174
        data['surface_albedo_assumed'] = grp.variables['surface_albedo_assumed'][0, scanline, pixel]
175
        data['snow_ice_flag'] = grp.variables['snow_ice_flag'][0, scanline, pixel]
176
        data['surface_pressure'] = grp.variables['surface_pressure'][0, scanline, pixel]
177
        data['surface_altitude'] = grp.variables['surface_altitude'][0, scanline, pixel]
178
        data['surface_altitude_precision'] = grp.variables['surface_altitude_precision'][0, scanline, pixel]
179

    
180
    return data
181

    
182
def write_spectrum(fname, wvl_refl=None, wvl_irr=None, wvl_rad=None, irr=None, rad=None,
183
                   refl=None, irr_noise=None, rad_noise=None, refl_noise=None,
184
                   sza=None, lat=None, lon=None, time=None, fresco=None):
185
    data = {}
186
    length = -1
187
    for s, n in zip((wvl_refl, wvl_irr, wvl_rad, irr, rad, refl, irr_noise, rad_noise, refl_noise),
188
                    ("wvl_refl", "wvl_irr", "wvl_rad", "irr", "rad", "refl", "irr_noise", "rad_noise", "refl_noise")):
189
        if s is not None:
190
            data[n] = s
191
            if len(s) > length:
192
                length = len(s)
193

    
194
    ref = open(fname, 'w')
195
    with ref:
196
        keys = list(data.keys())
197
        print(",".join(keys), file=ref)
198
        if sza is not None:
199
            print("# sza: {0}".format(sza), file=ref)
200
        if lat is not None:
201
            print("# latitude: {0}".format(lat), file=ref)
202
        if lon is not None:
203
            print("# longitude: {0}".format(lon), file=ref)
204
        if time is not None:
205
            print("# time: {0:%Y-%m-%dT%H:%M:%S.%f}".format(time), file=ref)
206
        if fresco is not None:
207
            for k, v in fresco.items():
208
                print("# {0}: {1}".format(k, v), file=ref)
209
        for idx in range(length):
210
            print(",".join([str(data[k][idx]) for k in keys]), file=ref)
211

    
212
def main(argv):
213
    parser = argparse.ArgumentParser(description="Plot spectrum")
214
    parser.add_argument('-V', '--version', action='version', version='%(prog)s 1.0')
215
    parser.add_argument('-p', '--pixel', dest='pixel', type=int,
216
                        help="pixel to show")
217
    parser.add_argument('-s', '--scanline', dest='scanline', type=int, default=0,
218
                        help='Scanline to show (for radiances)')
219
    parser.add_argument('-b', '--band', dest='band', type=int, default=1, choices=(1,2,3,4,5,6,7,8),
220
                        help='Band ID (1-8)')
221
    parser.add_argument('-r', '--reflectance', dest='refl', action='store_true',
222
                        help='Show reflectance (if both radiance & irradiance are given, otherwise show both radiance & irradiance separately)')
223
    parser.add_argument('-R', '--range', dest='wvlrange', type=float, nargs=2,
224
                        help='Wavelength selection to show')
225
    parser.add_argument('-n', '--noise', dest='noise', action='store_true',
226
                        help='Show noise spectrum')
227
    parser.add_argument('-e', '--error', dest='error', action='store_true',
228
                        help='Show error spectrum')
229
    parser.add_argument('-f', '--figure', default=None, help='Plot to file')
230
    parser.add_argument('-o', '--out', dest='outfile', default=None, help="Write spectrum to ascii file.")
231

    
232
    parser.add_argument('files', nargs="+", metavar="FILE",
233
                        help="The files to process")
234
    args = parser.parse_args(argv)
235

    
236
    radiances = [f for f in args.files if "L1B_RA_BD{0:1d}".format(args.band) in f]
237
    irradiances = [f for f in args.files if "L1B_IR" in f]
238
    fresco_files = [f for f in args.files if "L2__FRESCO" in f]
239

    
240
    irr_wvl, irr, irr_noise, irr_error, irr_sza = None, None, None, None, None
241
    if len(irradiances) > 0:
242
        for irr_file in irradiances:
243
            try:
244
                irr_wvl, irr, irr_noise, irr_error, irr_sza, lat, lon, time = get_spectrum(irradiances[0], args.band, args.pixel, scanline=0, radiance=False)
245
                print('Read irradiance')
246
                break
247
            except:
248
                pass
249

    
250
    rad_wvl, rad, rad_noise, rad_error, rad_sza = None, None, None, None, None
251
    if len(radiances) > 0:
252
        for rad_file in radiances:
253
            try:
254
                rad_wvl, rad, rad_noise, rad_error, rad_sza, lat, lon, time = get_spectrum(rad_file, args.band, args.pixel, scanline=args.scanline, radiance=True)
255
                print('Read data for band {0}'.format(args.band))
256
                break
257
            except:
258
                pass
259

    
260
    fresco_data = None
261
    if len(fresco_files) > 0:
262
        for fresco_file in fresco_files:
263
            try:
264
                fresco_data = get_fresco(fresco_file, args.pixel, args.scanline)
265
                print('read FRESCO results')
266
                break
267
            except:
268
                break
269

    
270
    if rad_wvl is not None and irr_wvl is not None:
271
        idx = np.logical_and(np.logical_and(np.logical_and(np.isfinite(rad_wvl), np.isfinite(rad)),
272
                                            np.logical_and(np.isfinite(rad_noise), np.isfinite(rad_error))),
273
                             np.logical_and(np.logical_and(np.isfinite(irr_wvl), np.isfinite(irr)),
274
                                            np.logical_and(np.isfinite(irr_noise), np.isfinite(irr_error))))
275
    elif rad_wvl is not None:
276
        idx = np.logical_and(np.logical_and(np.isfinite(rad_wvl), np.isfinite(rad)),
277
                             np.logical_and(np.isfinite(rad_noise), np.isfinite(rad_error)))
278
    elif irr_wvl is not None:
279
        idx = np.logical_and(np.logical_and(np.isfinite(irr_wvl), np.isfinite(irr)),
280
                             np.logical_and(np.isfinite(irr_noise), np.isfinite(irr_error)))
281
    else: # rad is None and irr is None:
282
        print("No radiance or irradiance was found", file=sys.stderr)
283
        sys.exit(1)
284

    
285
    if args.refl and irr_wvl is None:
286
        print("No irradiance was found to calculate reflectance", file=sys.stderr)
287
        sys.exit(1)
288

    
289
    if args.noise and not args.error:
290
        irr_dat = irr_noise
291
        rad_dat = rad_noise
292
        title_suffix = " Noise"
293
    elif args.error and not args.noise:
294
        irr_dat = irr_error
295
        rad_dat = rad_error
296
        title_suffix = " Error"
297
    elif args.error and args.noise:
298
        try:
299
            irr_dat = np.sqrt(irr_error**2 + irr_noise**2)
300
        except:
301
            irr_dat = None
302
        try:
303
            rad_dat = np.sqrt(rad_error**2 + rad_noise**2)
304
        except:
305
            rad_dat = None
306
        title_suffix = " Noise + Error"
307
    else:
308
        if args.refl:
309
            rad_spline = Spline(rad_wvl[idx], rad[idx])
310
            rad_at_irr = rad_spline(irr_wvl[idx])
311
            rad_dat = (rad_at_irr/irr[idx])*(np.pi/np.cos(rad_sza*np.pi/180.0))
312
            irr_dat = None
313
            title_suffix = ""
314
        else:
315
            irr_dat = irr[idx]
316
            rad_dat = rad[idx]
317
            title_suffix = ""
318

    
319

    
320
    if irr_dat is None:
321
        plt.plot(rad_wvl[idx], rad_dat, 'k-')
322
        if args.refl:
323
            plt.title("Reflectance ({0}, {1})".format(args.pixel, args.scanline))
324
        else:
325
            plt.title("Radiance" + title_suffix + " ({0}, {1})".format(args.pixel, args.scanline))
326
        plt.xlabel("Wavelength [nm]")
327
        plt.ylabel("Reflectance [1]" if args.refl else "Signal [mol/nm/m^2/s]")
328
        if args.wvlrange is not None:
329
            plt.xlim(*args.wvlrange)
330
    else:
331
        if rad_dat is None:
332
            plt.plot(irr_wvl[idx], irr_dat, 'k-')
333
            plt.title("Irradiance" + title_suffix + " ({0}, {1})".format(args.pixel, args.scanline))
334
            plt.xlabel("Wavelength [nm]")
335
            plt.ylabel("Reflectance [1]" if args.refl else "Signal [mol/nm/m^2/s]")
336
            if args.wvlrange is not None:
337
                plt.xlim(*args.wvlrange)
338
        else:
339
            f, axarr = plt.subplots(2, sharex=True)
340
            axarr[0].plot(rad_wvl[idx], rad_dat)
341
            axarr[0].set_ylabel('Radiance')
342
            axarr[0].set_title(title_suffix[1:] + " ({0}, {1})".format(args.pixel, args.scanline))
343
            axarr[1].plot(irr_wvl[idx], irr_dat)
344
            axarr[1].set_xlabel("Wavelength [nm]")
345
            axarr[1].set_ylabel('Irradiance')
346
            if args.wvlrange is not None:
347
                axarr[0].xlim(*args.wvlrange)
348
                axarr[1].xlim(*args.wvlrange)
349

    
350
    if args.figure is None:
351
        plt.show()
352
    else:
353
        plt.savefig(args.figure)
354

    
355
    if args.outfile is not None:
356
        rad_spline = Spline(rad_wvl[idx], rad[idx])
357
        rad_at_irr = rad_spline(irr_wvl[idx])
358
        refl = (rad_at_irr/irr[idx])*(np.pi/np.cos(rad_sza*np.pi/180.0))
359
        refl_noise = refl * np.sqrt((irr_noise[idx]/irr[idx])**2 + (rad_noise[idx]/rad[idx])**2)
360
        write_spectrum(args.outfile, wvl_refl=irr_wvl[idx], wvl_irr=irr_wvl[idx], wvl_rad=rad_wvl[idx],
361
                       irr=irr[idx], rad=rad[idx], refl=refl,
362
                       irr_noise=irr_noise[idx], rad_noise=rad_noise[idx], refl_noise=refl_noise,
363
                       sza=rad_sza, lat=lat, lon=lon, time=time, fresco=fresco_data)
364
if __name__ == "__main__":
365
    main(sys.argv)