Project

General

Profile

Revision 815:c0e5b3712092

IDc0e5b3712092
Parent 814:8290334af1af
Child 816:d52976b6ccfc

Added by Maarten Sneep about 2 years ago

New TM5 profile extractor

View differences:

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

  
3
# Extract NO2 profile from CTMANA or CTMFCT file. 
4
# 
5
# There are two modes for this tool: 
6
#
7
# 1. extract for a given time and geolocation (single point)
8
# 2. extract all profiles for a level 2 granule and add it to the L2 data file.
9
#    This requires that you have write access to the L2 data file. 
10
#
11
# One of the species can be selected: NO2, SO2 or HCHO.
12
# 
13

  
14
import os
15
import sys
16
import argparse
17
import datetime
18
from requests.structures import CaseInsensitiveDict
19

  
20
import inspect
21

  
22
import numpy as np
23
import netCDF4
24
from scipy import constants
25

  
26
named_geolocation_points = CaseInsensitiveDict({
27
    'KNMI':         {'lon': 5.177999,   'lat': 52.101544},
28
    'Cabauw':       {'lon': 4.926318,   'lat': 51.970268},
29
    'BIRA':         {'lon': 50.796781,  'lat': 4.356972},
30
    'IUP':          {'lon': 8.84964,    'lat': 53.103718},
31
    'DLR':          {'lon': 11.276202,  'lat': 48.084865},
32
    'ESTEC':        {'lon': 4.42015,    'lat': 52.217917},
33
    'ESRIN':        {'lon': 12.671675,  'lat': 41.827141},
34
    'RAL':          {'lon': -1.310875,  'lat': 51.576541},
35
    'FMI':          {'lon': 24.961061,  'lat': 60.203779},
36
    'MPI':          {'lon': 8.228973,   'lat': 49.991104},
37
    'Sodankyla':    {'lon': 26.588973,  'lat': 67.415999},
38
    'GSFC':         {'lon': -76.839854, 'lat': 38.992783},
39
    'Paramaribo':   {'lon': -55.204717, 'lat': 5.449259},
40
    'Jungfraujoch': {'lon': 7.985249,   'lat': 46.547463},
41
    'Xianghe':      {'lon': 116.960,    'lat': 39.750}})
42

  
43
class CTM:
44
    def __init__(self, fname:str, species:str='no2'):
45
        self.fname = fname
46
        self.slice_index = None
47
        self.weights = None
48
        self.spatial_index_lon = None
49
        self.spatial_index_lat = None
50
        self.species = species
51
        self.final_output_shape = None
52
        self.ref = None
53
    
54
    def __enter__(self): 
55
        self.ref = netCDF4.Dataset(self.fname, 'r')
56
        return self.ref
57
  
58
    def __exit__(self, exception_type, exception_value, traceback): 
59
        if self.ref is not None:
60
            self.ref.close()
61
            self.ref = None
62
    
63
    def close(self):
64
        self.__exit__(None, None, None)
65
    
66
    def find_slice(self, t: datetime.datetime) -> int:
67
        """Find the correct time index"""
68
        time_ref = self.ref["time"]
69
        times = time_ref[:]
70
        time_units = time_ref.units
71
        
72
        try:
73
            f_times = netCDF4.num2date(times, time_units, 
74
                        only_use_cftime_datetimes=False, 
75
                        only_use_python_datetimes=True)
76
        except TypeError:
77
            f_times = netCDF4.num2date(times, time_units)
78
        
79
        if t < f_times[0] or t > f_times[-1]:
80
            raise RuntimeError("Requested time not covered in CTM file.")
81
        else:
82
            self.slice_index = np.argmin([np.fabs((t - tt).total_seconds()) for tt in f_times])
83
        return self.slice_index
84
        
85
    def find_indices_weights(self, latitude, longitude):
86
        """Find indices and weights for the spatial interpolation
87

  
88
Result applies to flattened array
89
"""
90
        if isinstance(latitude, list):
91
            latitude = np.asarray(latitude, dtype=np.float64)
92
            longitude = np.asarray(longitude, dtype=np.float64)
93
        elif isinstance(latitude, (int, float)):
94
            latitude = np.asarray([latitude], dtype=np.float64)
95
            longitude = np.asarray([longitude], dtype=np.float64)
96
        
97
        if np.any(latitude > 90.0) or np.any(latitude < -90.0):
98
            raise RuntimeError("Latitude out of bounds")
99
        if np.prod(latitude.shape) != np.prod(longitude.shape):
100
            raise RuntimeError("Size of latitude and longitude must match")
101
        
102
        self.final_output_shape = list(self.ref['hyam'].shape)
103
        self.final_output_shape.extend(list(latitude.shape))
104
        
105
        latitude = latitude.flatten()
106
        longitude = longitude.flatten()
107
        
108
        # size of TM5 array
109
        jm, im = self.ref[self.species].shape[-2:]
110
        
111
        longitude = np.where(longitude < -180.0, longitude + 360.0, longitude)
112
        longitude = np.where(longitude >= 180.0, longitude - 360.0, longitude)
113
        dlon = 360.0/im
114
        firstcellmid = 0.5*dlon-180.0
115
        rx = (longitude-firstcellmid)/dlon
116
        ix = np.asarray(rx, dtype=np.int32)
117
        gdx = rx - ix
118
        ix = np.where(ix == im, 0, ix)
119
        
120
        dlat = 180.0/jm
121
        firstcellmid = 0.5*dlat-90.0
122
        ry = (latitude-firstcellmid)/dlat
123
        iy = np.asarray(ry, dtype=np.int32)
124
        gdy = ry - iy
125
        
126
        # the other corner
127
        ix2 = ix + 1
128
        ix2 = np.where(ix2 == im, 0, ix2)
129
        iy2 = iy + 1
130
        if np.any(iy2 >= jm):
131
            idy = iy2 >= jm
132
            iy[idy] = jm-2
133
            iy2[idy] = jm-1
134
            gdy[idy] = 1.0
135
        if np.any(iy < 0):
136
            idy = iy < 0
137
            iy[idy] = 0
138
            iy2[idy] = 1
139
            gdy[idy] = 1.0
140
        
141
        output_shape = [4]
142
        output_shape.extend(list(latitude.shape))
143
        output_shape = tuple(output_shape)
144
        ix4a = np.concatenate((ix, ix2, ix, ix2)).reshape(output_shape)
145
        iy4a = np.concatenate((iy, iy, iy2, iy2)).reshape(output_shape)
146
        w4a = np.concatenate(((1-gdx)*(1-gdy), gdx*(1-gdy), (1-gdx)*gdy, gdx*gdy)).reshape(output_shape)
147
        
148
        self.weights = w4a
149
        self.spatial_index_lon = ix4a
150
        self.spatial_index_lat = iy4a
151

  
152
        return ix4a, iy4a, w4a
153
    
154
    def get_pressures(self, surface_pressure):
155
        hyai = self.ref['hyai'][:]
156
        hybi = self.ref['hybi'][:]
157
        return hyai + hyai*surface_pressure
158
        
159
    def get_pressures_interfaces(self):
160
        """Get pressure at the interfaces of the model"""
161
        field = self.ref[self.species]
162
        pressure_shape = list(field.shape[1:])
163
        pressure_shape[0] += 1
164
        p_interfaces = np.zeros(pressure_shape, dtype=np.float32)
165
        ps = self.ref['ps'][self.slice_index, ...]
166
        hyai = self.ref['hyai'][:]
167
        hybi = self.ref['hybi'][:]
168
        
169
        for i, (a, b) in enumerate(zip(hyai, hybi)):
170
            p_interfaces[i, ...] = a + b*ps
171
            
172
        return p_interfaces
173

  
174
    def get_partial_column_field(self, p_interfaces=None):
175
        """Translate input field in volume mixing ratio to mol/m² per layer"""
176
        xmair = 28.94 # mass of air (gram/mol)
177
        # conversion factor from mixing ratio, pressure drop (in Pa) to columns (in 10^15 molec cm^-2)
178
        pdiff2moleccm2 = 1.0e-15 * 0.1 * constants.Avogadro/(constants.g * xmair)
179
        # conversion factor from mixing ratio, pressure drop (in Pa) to columns in mol/m2.
180
        pdiffmolm2 =  1000 /(constants.g * xmair)
181
        
182
        field = self.ref[self.species]
183
        ps = self.ref['ps']
184
        hyai = self.ref['hyai'][:]
185
        hybi = self.ref['hybi'][:]
186
        
187
        # partial columns
188
        partial_column_field = np.zeros(field.shape[1:], dtype=np.float32)
189
        if p_interfaces is None:
190
            p_interfaces = self.get_pressures_interfaces()
191
        
192
        for level in range(self.ref[self.species].shape[1]):
193
            dp = p_interfaces[level, ...] - p_interfaces[level+1, ...]
194
            partial_column_field[level, ...] = pdiffmolm2 * dp * field[self.slice_index, level, ...]
195
        
196
        return partial_column_field
197
    
198
    def get_tropopause_pressure(self, surface_pressure_external=None):
199
        t_idx = self.slice_index
200
        hyai = self.ref['hyai'][:]
201
        hybi = self.ref['hybi'][:]
202
        surface_pressure = self.ref["ps"][t_idx, ...]
203
        tropopause_layer_index = self.ref["tropopause_layer_index"][t_idx, ...]
204
        p_interfaces = self.get_pressures_interfaces()
205
        dest_shape = (self.weights.shape[1], )
206
        p_tropo = np.zeros(dest_shape, dtype=np.float32)
207
        
208
        for pix_idx in range(self.weights.shape[1]):
209
            if surface_pressure_external is not None:
210
                ps = surface_pressure_external.flat[pix_idx]
211
            else:
212
                ps = 0.0
213
                for i in range(4):
214
                    ps += (self.weights[i, pix_idx] * 
215
                            surface_pressure[self.spatial_index_lat[i, pix_idx], 
216
                                             self.spatial_index_lon[i, pix_idx]])
217
                    
218
            for i in range(4):
219
                idx = tropopause_layer_index[self.spatial_index_lat[i, pix_idx], 
220
                                             self.spatial_index_lon[i, pix_idx]] + 1 # top of layer
221
                
222
                p_tropo[pix_idx] += self.weights[i, pix_idx] * (hyai[idx] + ps * hybi[idx])
223
        
224
        p_tropo = p_tropo.reshape(self.final_output_shape[1:])
225
        return p_tropo
226
        
227
    def get_temperature(self, surface_pressure_external=None):
228
        t_idx = self.slice_index
229
        p_interfaces = self.get_pressures_interfaces()
230
        field = self.ref["t"][t_idx, ...]
231
        surface_pressure = self.ref["ps"][t_idx, ...]
232
        hyam = self.ref['hyam'][:]
233
        hybm = self.ref['hybm'][:]
234
        
235
        dest_shape = (field.shape[0], self.weights.shape[1])
236
        result_field = np.zeros(dest_shape, dtype=np.float32)
237
        result_pressure = np.zeros(dest_shape, dtype=np.float32)
238
        
239
        for pix_idx in range(self.weights.shape[1]):
240
            if surface_pressure_external is not None:
241
                ps = surface_pressure_external.flat[pix_idx]
242
            else:
243
                ps = 0.0
244
            for i in range(4):
245
                if surface_pressure_external is None:
246
                    ps += (self.weights[i, pix_idx] * 
247
                            surface_pressure[self.spatial_index_lat[i, pix_idx], 
248
                                             self.spatial_index_lon[i, pix_idx]])
249
                result_field[:, pix_idx] += (self.weights[i, pix_idx] * 
250
                        field[:, self.spatial_index_lat[i, pix_idx], 
251
                                         self.spatial_index_lon[i, pix_idx]])
252
            
253
            result_pressure[:, pix_idx] = hyam + ps * hybm
254
        
255
        result_field = result_field.reshape(self.final_output_shape)
256
        result_pressure = result_pressure.reshape(self.final_output_shape)
257
        
258
        return result_field, result_pressure
259
        
260
    def get_profiles(self, vmr=True, surface_pressure_external=None):
261
        t_idx = self.slice_index
262
        
263
        p_interfaces = self.get_pressures_interfaces()
264
        if not vmr:
265
            field = self.get_partial_column_field(p_interfaces)
266
        else:
267
            field = self.ref[self.species][t_idx, ...]
268
        surface_pressure = self.ref["ps"][t_idx, ...]
269
        hyai = self.ref['hyai'][:]
270
        hybi = self.ref['hybi'][:]
271
        
272
        dest_shape = (field.shape[0], self.weights.shape[1])
273
        result_field = np.zeros(dest_shape, dtype=np.float32)
274
        pressure_shape = (field.shape[0]+1, self.weights.shape[1])
275
        result_pressure = np.zeros(pressure_shape, dtype=np.float32)
276
                
277
        for pix_idx in range(self.weights.shape[1]):
278
            if surface_pressure_external is not None:
279
                ps = surface_pressure_external.flat[pix_idx]
280
            else:
281
                ps = 0.0
282
            for i in range(4):
283
                if surface_pressure_external is None:
284
                    ps += (self.weights[i, pix_idx] * 
285
                            surface_pressure[self.spatial_index_lat[i, pix_idx], 
286
                                             self.spatial_index_lon[i, pix_idx]])
287
                result_field[:, pix_idx] += (self.weights[i, pix_idx] * 
288
                        field[:, self.spatial_index_lat[i, pix_idx], 
289
                                         self.spatial_index_lon[i, pix_idx]])
290
            
291
            result_pressure[:, pix_idx] = hyai + ps * hybi
292
        
293
        result_field = result_field.reshape(self.final_output_shape)
294
        pressure_shape = list(self.final_output_shape)
295
        pressure_shape[0] += 1
296
        result_pressure = result_pressure.reshape(tuple(pressure_shape))
297
        
298
        return result_field, result_pressure
299
        
300
        
301
        
302
def get_function_name():
303
   return inspect.stack()[1][3]
304

  
305

  
306
def extract_single(species:str='no2', ctm:str=None, vmr:bool=False, 
307
        time:datetime.datetime=None, geo_point:tuple=None, 
308
        pressure:float=None, verbose:int=0, **kwargs):
309
    ctmobj = CTM(ctm, species)
310
    with ctmobj:
311
        try:
312
            ctmobj.find_slice(time)
313
        except Exception as err:
314
            print(err, file=sys.stderr)
315
            sys.exit(1)
316
        ctmobj.find_indices_weights(geo_point[1], geo_point[0])
317
        if pressure is None:
318
            pres_array = None
319
        else:
320
            pres_array = np.asarray([100*pressure]) # hPa -> Pa.
321
        result_field, result_pressure  = ctmobj.get_profiles(vmr, surface_pressure_external=pres_array)
322
        result_temperature, result_pressure_mid = ctmobj.get_temperature(pres_array)
323
        
324
        tropopause_pressure = ctmobj.get_tropopause_pressure(pres_array)
325
    
326
    if vmr:
327
        print(f"pressure_bottom\tpressure_top\tpressure_mid\ttemperature\t{species}_vmr\n# hPa\t# hPa \t# hPa\t# K\t# 1")
328
    else:
329
        print(f"pressure_bottom\tpressure_top\tpressure_mid\ttemperature\t{species}_partial_columns\n# hPa\t# hPa \t# hPa\t# K\t# mol/m2 in layer")
330
    print(f"# time: {time:%Y-%m-%dT%H:%M:%SZ}, latitude: {geo_point[1]}, longitude: {geo_point[0]}")
331
    print(f"# source: {os.path.basename(ctm)}")
332
    print(f"# tropopause pressure: {tropopause_pressure[0]/100} hPa")
333
    
334
    for pres_bot, pres_top, pres_mid, t, val in zip(result_pressure[0:-1], result_pressure[1:], result_pressure_mid, result_temperature, result_field):
335
        print(f"{pres_bot[0]/100}\t{pres_top[0]/100}\t{pres_mid[0]/100}\t{t[0]}\t{val[0]}")
336
    
337
    
338
def extract_orbit(species:str='no2', ctm:str=None, vmr:bool=False,
339
        l2_file:str=None, use_TM5_pressure:bool=False, verbose:int=0, 
340
        add_temperature_profile:bool=False, **kwargs):
341
    
342
    if not os.access(l2_file, os.O_RDWR):
343
        print(f"""Error when opening '{os.path.basename(l2_file)}' in append mode. 
344
We need write permission to add the model profiles to the level 2 file.""", file=sys.stderr)
345
        return
346
    
347
    ctmobj = CTM(ctm, species)
348
    ref = netCDF4.Dataset(l2_file, 'a')
349
    with ref, ctmobj:
350
        tstart = datetime.datetime.strptime(ref.time_coverage_start[0:19], '%Y-%m-%dT%H:%M:%S')
351
        tend = datetime.datetime.strptime(ref.time_coverage_end[0:19], '%Y-%m-%dT%H:%M:%S')
352
        time = tstart + (tend-tstart)/2
353
        
354
        if verbose > 0:
355
            print("Reading L2 data", file=sys.stderr)
356
            
357
        latitudes = ref['/PRODUCT/latitude'][:]
358
        longitudes = ref['/PRODUCT/longitude'][:]
359
        if use_TM5_pressure:
360
            pressures = None
361
        else:
362
            pressures = ref['/PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_pressure'][:]
363
            
364
        try:
365
            ctmobj.find_slice(time)
366
        except Exception as err:
367
            print(err, file=sys.stderr)
368
            sys.exit(1)
369
        
370
        if verbose > 0:
371
            print("Calculating weights", file=sys.stderr)
372
        ctmobj.find_indices_weights(latitudes, longitudes)
373
        
374
        if verbose > 0:
375
            print("Retrieving profiles", file=sys.stderr)
376
        result_field, result_pressure = ctmobj.get_profiles(vmr,
377
            surface_pressure_external=pressures)
378
        
379
        if verbose > 1:
380
            print(f"""NO2 dimensions: {result_field.shape}
381
pressure bounds dimensions: {result_pressure.shape}""")
382
            print("Writing output", file=sys.stderr)
383

  
384
        name_map = {'no2': 'nitrogen_dioxide', 'so2': 'sulfur_dioxide', 'hcho': 'formaldehyde'}
385
        grp = ref['/PRODUCT/SUPPORT_DATA/DETAILED_RESULTS']
386
        vname = f"{species}_vmr" if vmr else f"{species}_partial_columns"
387
        try:
388
            var = grp.createVariable(vname, np.float32,
389
                                    dimensions=('time', 'scanline', 'ground_pixel', 'layer'),
390
                                    zlib=True, fletcher32=True)
391
        except RuntimeError:
392
            var = grp[vname]
393
            
394
        var[:] = result_field.transpose((1,2,3,0))
395
        var.units = "1" if vmr else "mol/m2"
396
        var.long_name = (f"TM5-MP a-priori profile of mixing ratios of {species.upper()} co-located with the observations" 
397
                         if vmr else 
398
                         f"TM5-MP a-priori profile of model grid-cell partial columns of {species.upper()} co-located with the observations")
399
        var.standard_name = (f"mole_fraction_of_{name_map[species]}_in_air" 
400
                             if vmr else 
401
                             f"mole_content_of_{name_map[species]}_in_atmosphere_layer")
402
        var.coordinates = "/PRODUCT/longitude /PRODUCT/latitude"
403
        if use_TM5_pressure:
404
            var.ancillary_variables = "/PRODUCT/tm5_constant_a /PRODUCT/tm5_constant_b /PRODUCT/tm5_tropopause_layer_index /PRODUCT/SUPPORT_DATA/INPUT_DATA/TM5_surface_pressure"
405
        else:
406
            var.ancillary_variables = "/PRODUCT/tm5_constant_a /PRODUCT/tm5_constant_b /PRODUCT/tm5_tropopause_layer_index /PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_pressure"
407
        
408
        if use_TM5_pressure:
409
            try:
410
                var = grp.createVariable("TM5_surface_pressure", np.float32,
411
                                    dimensions=('time', 'scanline', 'ground_pixel'),
412
                                    zlib=True, fletcher32=True)
413
            except RuntimeError:
414
                var = grp["TM5_surface_pressure"]
415
            var[:] = result_pressure[0, ...]
416
            var.units = "Pa"
417
            var.long_name = f"Surface pressure used to calculate pressure grid for {vname} profiles"
418
            var.standard_name = 'surface_air_pressure'
419
            var.coordinates = "/PRODUCT/longitude /PRODUCT/latitude"
420
        
421
        if add_temperature_profile:
422
            if verbose > 0:
423
                print("Extracting temperatures", file=sys.stderr)
424
            result_temperature, result_pressure_mid = ctmobj.get_temperature()
425
            if verbose > 1:
426
                print(f"""temperature dimensions: {result_temperature.shape}
427
mid pressure dimensions: {result_pressure_mid.shape}""", file=sys.stderr)
428
            if verbose > 0:
429
                print("Writing temperatures", file=sys.stderr)
430
            vname = "temperature"
431
            try:
432
                var = grp.createVariable(vname, np.float32,
433
                                    dimensions=('time', 'scanline', 'ground_pixel', 'layer'),
434
                                    zlib=True, fletcher32=True)
435
            except RuntimeError:
436
                var = grp[vname]
437
            var[:] = result_temperature.transpose((1,2,3,0))
438
            var.units = "K"
439
            var.long_name = "TM5-MP temperature profile"
440
            var.standard_name = "air_temperature"
441
            var.coordinates = "/PRODUCT/longitude /PRODUCT/latitude"
442
            if use_TM5_pressure:
443
                var.ancillary_variables = "/PRODUCT/tm5_constant_a /PRODUCT/tm5_constant_b /PRODUCT/tm5_tropopause_layer_index /PRODUCT/SUPPORT_DATA/INPUT_DATA/TM5_surface_pressure"
444
            else:
445
                var.ancillary_variables = "/PRODUCT/tm5_constant_a /PRODUCT/tm5_constant_b /PRODUCT/tm5_tropopause_layer_index /PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_pressure"
446
    
447
def argparse_datetime(s: str):
448
    try:
449
        value = datetime.datetime.strptime(s, "%Y-%m-%dT%H:%M:%S")
450
    except ValueError as err:
451
        msg = "{0} does not match format '%Y-%m-%dT%H:%M:%S'".format(s)
452
        raise argparse.ArgumentTypeError(msg)
453
    return value
454

  
455
def geospec_verify(s: str):
456
    global named_geolocation_points
457
    if s not in named_geolocation_points:
458
        names = "', '".join(sorted(list(named_geolocation_points.keys())))
459
        msg = "'{0}' not found. \nValid point names are '{1}'.".format(s, names)
460
        raise argparse.ArgumentTypeError(msg)
461
    return s
462

  
463
def geo_point_name(name: str):
464
    global named_geolocation_points
465
    # https://mynasadata.larc.nasa.gov/latitudelongitude-finder/
466
    if name in named_geolocation_points:
467
        return (named_geolocation_points[name]['lon'], named_geolocation_points[name]['lat'])
468
    else:
469
        names = "', '".join(sorted(list(named_geolocation_points.keys())))
470
        raise RuntimeError("Warning, '{0}' not found. \nValid point names are '{1}'.".format(name, names))
471
        
472

  
473
def main():
474
    global named_geolocation_points
475
    geolocation_point_names = "', '".join(sorted(list(named_geolocation_points.keys())))
476
    
477
    parser = argparse.ArgumentParser(description='Get profile(s) for NO2')
478
    parser.add_argument('-v', '--verbose', action='count', default=0, help="Be chatty")
479
    parser.add_argument('-s', '--species', metavar='MOLECULE', 
480
                        choices=('no2', 'so2', 'ch2o', 'hcho', 'NO2', 'SO2', 'CH2O', 'HCHO', 'NO₂', 'SO₂', 'CH₂O'),
481
                        default='no2', 
482
                        help='Trace gas to extract, default NO2, valid arguments: so2, no2, hcho, ch2o')
483
    parser.add_argument('-c', '--ctm', required=True, help='The CTM file with the profiles.')
484
    parser.add_argument('-V', '--vmr', action='store_true',
485
                        help="Store profile as vmr rather than partial columns")
486
    
487
    subparsers = parser.add_subparsers(help='Processing mode', required=True)
488
    
489
    parser_s = subparsers.add_parser('single', help='Extract the profile for a single time and location')
490
    parser_s.set_defaults(func=extract_single)
491
    parser_s.add_argument('-t', '--time', metavar='YYYY-MM-DDTHH:MM:SS', type=argparse_datetime,
492
                          required=True,
493
                          help='Date/time to extract for')
494
    group = parser_s.add_mutually_exclusive_group(required=True)
495
    group.add_argument('-g', '--geo-point', metavar="L", type=float, dest="geo_point", 
496
                          nargs=2, 
497
                          help="Specify the location to extract a profile for as a LAT LON pair")
498
    group.add_argument('-G', '--geo-name', metavar="NAME", type=geospec_verify, dest="geo_point_name", 
499
                          help="Specify the location to extract a profile for by name. Recognized names are: '{0}'.".format(geolocation_point_names))
500
    parser_s.add_argument('-p', '--pressure', type=float, default=None, 
501
                          help="Surface pressure for profile. Use TM5 surface pressure if not given.")
502
    
503
    parser_m = subparsers.add_parser('level2', help='Extract profiles for a full level 2 file. Profiles will be added to the file, read/write access is required.')    
504
    parser_m.set_defaults(func=extract_orbit)
505
    parser_m.add_argument('-i', '--input', metavar='FILE', dest='l2_file', required=True, help='L2 NO2 file to add data to')
506
    parser_m.add_argument('--use-TM5-pressure', action='store_true',
507
                        help="Use the surface pressure from TM5 rather than the L2 surface pressure")
508
    parser_m.add_argument('--add-temperature-profile', action='store_true',
509
                        help="Also add the temperature profile from TM5 to the output.")
510

  
511
    args = vars(parser.parse_args())
512
    
513
    processor = args["func"]
514
    del args["func"]
515
    
516
    if "geo_point_name" in args and args["geo_point_name"]:
517
        args["geo_point"] = geo_point_name(args["geo_point_name"])
518
        del args["geo_point_name"]
519
    elif "geo_point" in args and args["geo_point"]:
520
        args["geo_point"] = args["geo_point"][1], args["geo_point"][0]
521
    
522
    if args["species"] in ("NO2", "no2", "NO₂"):
523
        args["species"] = "no2"
524
    elif args["species"] in ("SO2", "so2", "SO₂"):
525
        args["species"] = "so2"
526
    else:
527
        args["species"] = "ch2o"
528
        
529
    processor(**args)
530
    
531
if __name__ == "__main__":
532
    main()
src/TM5_profile_extractor.py
1
#!/usr/bin/env python3
2

  
3
import os
4
import sys
5
import datetime
6
import glob
7

  
8
import netCDF4
9
import numpy as np
10

  
11

  
12
class NO2Profile:
13
    """Extract NO2 (or SO2, CH2O) profiles from TM5 file.
14

  
15
matching the time & location of a L2 observation.
16

  
17
After initialization, you need to select the correct input profile using the
18
select_input_profiles() method.
19

  
20
The actual extraction can then be performed using the extract_partial_columns() method.
21

  
22
The resulting profiles can then be found in profiles instance variable,
23
and integrated columns in the columns instance variable. If the surface pressures
24
from TM5 are used, the interpolated surface pressure can be found in the pressures
25
instance variable. To reconstruct the pressure levels (either mid-point of interfaces)
26
you can use p = a + ps * b, with a and b from the hyai/hybi instance variables for
27
the interfaces, or from hyam/hybm instance variables for the mid points.
28

  
29
The last dimension contains the profiles, all dimensions before that follow the
30
input dimensions, the columns and surface pressure follow the shape of the input
31
parameters.
32

  
33
The units and (tentavive) standard_name values can be found in the columns_units,
34
columns_standard_name, profiles_units, profiles_standard_name, pressure_units and
35
pressure_standard_name properties.
36

  
37
The other methods are for internal use only.
38

  
39
The command-line tool included at the end demonstrates useage of this class, and
40
doubles as a tool to add profiles to an existing NO2 level 2 file.
41
The tool has a --help argument.
42

  
43
Note that using the surface pressure from the level 2 file will add considerable
44
time to the profile calculation.
45
"""
46
    def __init__(self, aux_files):
47
        """Initializer.
48

  
49
The single argument is the list of candidate auxiliary files from which the
50
profiles should be extracted
51
"""
52
        self.aux_files = aux_files
53
        self.active_file = None
54
        self.active_slice = None
55
        self.offset = datetime.datetime.strptime("1950-01-01T00:00:00", "%Y-%m-%dT%H:%M:%S")
56
        self.species = 'no2'
57

  
58
    def select_input_profiles(self, time=None, level2File=None):
59
        """Find auxiliary input file.
60

  
61
Either use the time range in a L2 file, or use the time (ISO string
62
"YYYY-MM-DDTHH:MM:SS" or python datetime object).
63

  
64
Sets self.active_file to the file in self.aux_files that matches,
65
and sets self.active_slice to the actual desired index.
66

  
67
The correct time slice for NO2 is the slice closest to the mid-time
68
of the time-range covered by the full orbit. This time is extracted
69
here, but only in the case of an offline product. It can be
70
overridden by supplying a time.
71
        """
72

  
73
        if level2File is not None:
74
            files = []
75
            ref = netCDF4.Dataset(level2File, 'r')
76
            with ref:
77
                start, stop = [datetime.datetime.strptime(tt[:19], "%Y-%m-%dT%H:%M:%S") for tt in (ref.time_coverage_start, ref.time_coverage_end)]
78
            time = start + (stop-start)/2
79
        elif time is not None:
80
            try:
81
                t = datetime.datetime.strptime(time[:19], "%Y-%m-%dT%H:%M:%S")
82
            except TypeError:
83
                t = time
84
            except ValueError:
85
                raise RuntimeError("Input time does not match template")
86

  
87
            for f in self.aux_files:
88
                ref = netCDF4.Dataset(f, 'r')
89
                with ref:
90
                    start, stop = [datetime.datetime.strptime(tt, "%Y%m%dT%H%M%S") for tt in (ref.validity_start, ref.validity_stop)]
91
                    if start < t < stop:
92
                        self.active_file = f
93
                        break
94
            if self.active_file is None:
95
                raise RuntimeError("No suitable file was found for time {0:%Y-%m-%dT%H:%M:%S}".format(t))
96
            self.select_time(t)
97
            return self.active_file
98
        else:
99
            raise RuntimeError("Must provide time or level 2 file (not both)")
100

  
101
    def select_time(self, theTime):
102
        """Determine the time slice to use based on the supplied time."""
103
        ref = netCDF4.Dataset(self.active_file, 'r')
104
        with ref:
105
            self.time_axis = [self.offset + datetime.timedelta(days=tt) for tt in ref.variables['time'][:]]
106
            dt = np.asarray([np.fabs((t-theTime).total_seconds()) for t in self.time_axis], dtype=np.float64)
107
            self.active_slice = np.argmin(dt)
108
            self.time = self.time_axis[self.active_slice]
109

  
110
    def extract_vmr(self, latitude, longitude, trace_gas="no2"):
111
        """Interpolate the profiles for the requested latitude/longitude values as volume mixing ratios.
112

  
113
The latitude and longitude can be numpy arrays, even with higher dimensions,
114
such as the latitude/longitude parameters as found in a L2 file.
115

  
116
The trace gas parameter selects the species to extract from the TM5 file,
117
can be one of no2, so2 or ch2o (hcho). Default is no2.
118

  
119
This method checks that the aux file and time-slice have been properly selected.
120

  
121
After calling this method the profiles can be read from the self.profiles variable.
122
The last dimension contains the profiles, all dimensions before that follow the
123
input dimensions.
124

  
125
In this method the total colum is not calculated.
126
"""
127
        if self.active_file is None or self.active_slice is None:
128
            raise RuntimeError("No active file or active slice available.")
129

  
130
        self.species = trace_gas if trace_gas != "hcho" else "ch2o"
131
        self.prof_type = 'volume_mixing_ratio'
132

  
133
        latitude = np.asarray(latitude)
134
        longitude = np.asarray(longitude)
135

  
136
        try:
137
            npix = list(latitude.shape)
138
        except:
139
            npix = [1]
140

  
141
        ref = netCDF4.Dataset(self.active_file, 'r')
142
        with ref:
143
            self.latitude = ref.variables['lat'][:]
144
            self.longitude = ref.variables['lon'][:]
145
            self.hyai = ref.variables['hyai'][:]
146
            self.hybi = ref.variables['hybi'][:]
147
            self.hyam = ref.variables['hyam'][:]
148
            self.hybm = ref.variables['hybm'][:]
149

  
150
            ps = ref.variables['ps'][self.active_slice, :, :]
151
            vmr = ref.variables[self.species][self.active_slice, :, :, :]
152
            vmr = np.transpose(vmr, (1,2,0))
153

  
154
        nlev = self.hyai.shape[0]-1
155
        profile_shape = npix[:]
156
        profile_shape.append(nlev)
157
        self.profiles = np.ma.zeros(tuple(profile_shape), dtype=np.float32)
158
        self.pressure = np.zeros(tuple(npix), dtype=np.float32)
159

  
160
        ix4a, iy4a, w4a = self.getCellIndices(latitude, longitude)
161

  
162

  
163
        for jj in range(np.prod(npix)): # loop over coordinates
164
            for ix, iy, w in zip(ix4a, iy4a, w4a): # loop over 4 surrounding pixels
165
                self.pressure.flat[jj] += ps[iy.flat[jj], ix.flat[jj]] * w.flat[jj]
166
                self.profiles.flat[nlev*jj:nlev*(jj+1)] += w.flat[jj]*vmr[iy.flat[jj], ix.flat[jj], :]
167

  
168
        self.columns = None
169
        self.pressure = None
170

  
171
    def convert_vmr_to_partial_columns(self, surface_pressure=None):
172
        """Convert the colocated vmr profiles to partial columns"""
173
        grav = 9.80665
174
        Navog = 6.02214076e23
175
        xmair = 28.94
176
        pdiff2molm2 = 1.66054e-21 * Navog / ( grav * xmair)
177

  
178
        npix = list(self.profiles.shape[0:-1])
179
        nlev = self.hyai.shape[0]-1
180
        profile_shape = npix[:]
181
        profile_shape.append(nlev)
182
        pp_shape = npix[:]
183
        pp_shape.append(nlev+1)
184

  
185
        pp = np.zeros(tuple(pp_shape), dtype=np.float32)
186
        no2pcf = np.ma.zeros(tuple(profile_shape), dtype=np.float32)
187
        no2pcf.mask = np.zeros(tuple(profile_shape), dtype=np.bool)
188

  
189
        if surface_pressure is None:
190
            surface_pressure = self.pressure
191
            mask_data = False
192
        else:
193
            mask_data = True
194
        for i, (a, b) in enumerate(zip(self.hyai, self.hybi)):
195
            pp[..., i] = surface_pressure * b + a
196
        for i in range(self.hybi.shape[0]-1):
197
            no2pcf[..., i] = pdiff2molm2 * (pp[..., i] - pp[..., i+1]) * self.profiles[..., i]
198
            if mask_data:
199
                no2pcf.mask[..., i] = surface_pressure.mask
200

  
201
        self.profiles = no2pcf
202
        self.columns = self.profiles.sum(axis=-1)
203
        self.prof_type = "partial_columns"
204

  
205
    def extract_partial_columns(self, latitude, longitude, surface_pressure=None, trace_gas="no2"):
206
        """Interpolate the profiles for the requested latitude/longitude values as partial columns.
207

  
208
The latitude and longitude can be numpy arrays, even with higher dimensions,
209
such as the latitude/longitude parameters as found in a L2 file. Actual
210
surface pressure can be provided, otherwise the TM5 surface pressure will be
211
interpolated. When provided the surface pressure is assumed to be in Pa.
212

  
213
The trace gas parameter selects the species to extract from the TM5 file,
214
can be one of no2, so2 or ch2o. Default is no2.
215

  
216
This method checks that the aux file and time-slice have been properly selected.
217

  
218
After calling this method the profiles can be read from the self.profiles variable.
219
The last dimension contains the profiles, all dimensions before that follow the
220
input dimensions.
221

  
222
Similarly the columns can be read from self.columns, and the TM5 surface pressures
223
can be read from self.pressure (only if no input pressures were provided,
224
otherwise it will be a copy of the input pressures).
225
"""
226
        if self.active_file is None or self.active_slice is None:
227
            raise RuntimeError("No active file or active slice available.")
228

  
229
        grav = 9.80665
230
        Navog = 6.02214076e23
231
        xmair = 28.94
232
        pdiff2molm2 = 1.66054e-21 * Navog / ( grav * xmair)
233

  
234
        self.species = trace_gas if trace_gas != "hcho" else "ch2o"
235
        self.prof_type = 'partial_columns'
236

  
237
        latitude = np.asarray(latitude)
238
        longitude = np.asarray(longitude)
239

  
240
        try:
241
            npix = list(latitude.shape)
242
        except:
243
            npix = [1]
244

  
245
        ref = netCDF4.Dataset(self.active_file, 'r')
246
        with ref:
247
            self.latitude = ref.variables['lat'][:]
248
            self.longitude = ref.variables['lon'][:]
249
            self.hyai = ref.variables['hyai'][:]
250
            self.hybi = ref.variables['hybi'][:]
251
            self.hyam = ref.variables['hyam'][:]
252
            self.hybm = ref.variables['hybm'][:]
253

  
254
            ps = ref.variables['ps'][self.active_slice, :, :]
255
            vmr = ref.variables[self.species][self.active_slice, :, :, :]
256

  
257
        nlev = self.hyai.shape[0]-1
258
        profile_shape = npix[:]
259
        profile_shape.append(nlev)
260
        pp_shape = npix[:]
261
        pp_shape.append(nlev+1)
262
        self.profiles = np.ma.zeros(tuple(profile_shape), dtype=np.float32)
263

  
264
        if surface_pressure is None:
265
            self.pressure = np.zeros(tuple(npix), dtype=np.float32)
266
            calc_surface_pressure = True
267
        else:
268
            self.pressure = surface_pressure
269
            calc_surface_pressure = False
270
            self.profiles.mask = np.zeros(tuple(profile_shape), dtype=np.bool)
271

  
272
        ix4a, iy4a, w4a = self.getCellIndices(latitude, longitude)
273

  
274

  
275
        if calc_surface_pressure:
276
            pp = np.zeros(tuple(pp_shape), dtype=np.float32)
277
            no2pcf = np.zeros(tuple(profile_shape), dtype=np.float32)
278

  
279
            for i, (a,b) in enumerate(zip(self.hyai, self.hybi)):
280
                pp[:, :, i] = ps * b + a
281
            for i in range(self.hybi.shape[0]-1):
282
                no2pcf[:, :, i] = pdiff2molm2 * (pp[:, :,i] - pp[:, :, i+1]) * vmr[i, :, :]
283

  
284
            pp = None
285
        else:
286
            # contiguous access, saves ~5 minutes on 22 total when converting a whole orbit
287
            vmr = np.transpose(vmr, (1,2,0))
288

  
289
        for jj in range(np.prod(npix)): # loop over coordinates
290
            try:
291
                # skip fill values in the surface pressure.
292
                if (not calc_surface_pressure) and surface_pressure.mask.flat[jj]:
293
                    self.profiles.mask.flat[nlev*jj:nlev*(jj+1)] = True
294
                    continue
295
            except:
296
                pass
297

  
298
            if not calc_surface_pressure:
299
                pp = surface_pressure.flat[jj] * self.hybi + self.hyai
300

  
301
            for ix, iy, w in zip(ix4a, iy4a, w4a): # loop over 4 surrounding pixels
302
                if calc_surface_pressure:
303
                    self.pressure.flat[jj] += ps[iy.flat[jj], ix.flat[jj]] * w.flat[jj]
304
                    self.profiles.flat[nlev*jj:nlev*(jj+1)] += w.flat[jj]*no2pcf[iy.flat[jj], ix.flat[jj], :]
305
                else:
306
                    no2pcf = pdiff2molm2 * (pp[0:-1] - pp[1:]) * vmr[iy.flat[jj], ix.flat[jj], :]
307
                    self.profiles.flat[nlev*jj:nlev*(jj+1)] += w.flat[jj]*no2pcf
308

  
309
        self.columns = self.profiles.sum(axis=-1)
310

  
311
    @property
312
    def columns_units(self):
313
        """Units of the columns"""
314
        if self.prof_type == "partial_columns":
315
            return "mol/m2"
316
        else:
317
            return None
318

  
319
    @property
320
    def profiles_units(self):
321
        """Units of the profiles"""
322
        if self.prof_type == "partial_columns":
323
            return "mol/m2"
324
        else:
325
            return "1"
326

  
327
    @property
328
    def columns_standard_name(self):
329
        """Standard name of the columns. Note that some are 'extrapolated' following existing names"""
330
        if self.prof_type == "partial_columns":
331
            if self.species == 'no2':
332
                return "atmosphere_mole_content_of_nitrogen_dioxide"
333
            elif self.species == 'so2':
334
                return "atmosphere_mole_content_of_sulfur_dioxide"
335
            else:
336
                return "atmosphere_mole_content_of_formaldehyde"
337
        else:
338
            return None
339

  
340
    @property
341
    def profiles_standard_name(self):
342
        """Standard name of the profiles. Note that some are 'extrapolated' following existing names"""
343
        if self.prof_type == "partial_columns":
344
            if self.species == 'no2':
345
                return "mole_content_of_nitrogen_dioxide_in_atmosphere_layer"
346
            elif self.species == 'so2':
347
                return "mole_content_of_sulfur_dioxide_in_atmosphere_layer"
348
            else:
349
                return "mole_content_of_formaldehyde_in_atmosphere_layer"
350
        else:
351
            if self.species == 'no2':
352
                return "mole_fraction_of_nitrogen_dioxide_in_air"
353
            elif self.species == 'so2':
354
                return "mole_fraction_of_sulfur_dioxide_in_air"
355
            else:
356
                return "mole_fraction_of_formaldehyde_in_air"
357

  
358
    @property
359
    def profiles_long_name(self):
360
        """Long name of the profiles. """
361
        if self.prof_type == "partial_columns":
362
            if self.species == 'no2':
363
                return "TM5-MP a-priori profile of model grid-cell partial columns of NO2 co-located with the observations"
364
            elif self.species == 'so2':
365
                return "TM5-MP a-priori profile of model grid-cell partial columns of SO2 co-located with the observations"
366
            else:
367
                return "TM5-MP a-priori profile of model grid-cell partial columns of HCHO co-located with the observations"
368
        else:
369
            if self.species == 'no2':
370
                return "TM5-MP a-priori profile of mixing ratios of NO2 co-located with the observations"
371
            elif self.species == 'so2':
372
                return "TM5-MP a-priori profile of mixing ratios of SO2 co-located with the observations"
373
            else:
374
                return "TM5-MP a-priori profile of mixing ratios of HCHO co-located with the observations"
375

  
376

  
377
    @property
378
    def pressure_units(self):
379
        """Units of the surface pressure"""
380
        return "Pa"
381

  
382
    @property
383
    def pressure_standard_name(self):
384
        """Standard name of the surface pressure"""
385
        return "surface_air_pressure"
386

  
387
    def getCellIndices(self, latitude, longitude):
388
        """Get the indices and weights
389

  
390
Returns tuple of ix4a, iy4a, w4a: index in TM5 longitude dimension,
391
index in TM5 latitude dimension and fractional weights.
392
Each element is a list of 4 values, each value is an array with the
393
same dimensions as the input.
394
"""
395
        xlon2 = longitude[:]
396
        xlon2 = np.where(xlon2 < -180.0, xlon2 + 360.0, xlon2)
397
        xlon2 = np.where(xlon2 >= 180.0, xlon2 - 360.0, xlon2)
398

  
399
        dlon = self.longitude[1] - self.longitude[0]
400
        rx = (xlon2-self.longitude[0])/dlon
401
        ix = np.asarray(np.floor(rx), dtype=np.int32)
402
        gdx = rx - ix
403
        ix = np.where(ix == self.longitude.shape[0], 0, ix)
404
        ix = np.where(ix == 0, self.longitude.shape[0] - 1, ix)
405

  
406
        dlat = self.latitude[1] - self.latitude[0]
407
        ry = (latitude-self.latitude[0])/dlat
408
        iy = np.asarray(np.floor(ry), dtype=np.int32)
409
        gdy = ry - iy
410

  
411
        ix2 = np.where(ix + 1 == self.longitude.shape[0], 0, ix + 1)
412

  
413
        if np.any(iy == self.latitude.shape[0] - 1):
414
            iy2 = np.where(iy + 1 >= self.latitude.shape[0], self.latitude.shape[0]-1, iy + 1)
415
            iy = np.where(iy + 1 >= self.latitude.shape[0], self.latitude.shape[0]-2, iy)
416
            gdy = np.where(iy + 1 == self.latitude.shape[0], 1.0, gdy)
417
        elif np.any(iy < 0):
418
            iy2 = np.where(iy <= 0, 1, iy+1)
419
            gdy = np.where(iy <= 0, 1.0, gdy)
420
            iy = np.where(iy <= 0, 0, iy)
421
        else:
422
            iy2 = iy + 1
423

  
424
        ix4a = [ix, ix2, ix, ix2]
425
        iy4a = [iy, iy, iy2, iy2]
426
        w4a = [(1.0-gdx)*(1.0-gdy), gdx*(1.0-gdy), (1.0-gdx)*gdy, gdx*gdy]
427

  
428
        return ix4a, iy4a, w4a
429

  
430

  
431
if __name__ == "__main__":
432
    import argparse
433

  
434
    def argparse_datetime(s):
435
        try:
436
            value = datetime.datetime.strptime(s, "%Y-%m-%dT%H:%M:%S")
437
        except ValueError as err:
438
            msg = "{0} does not match format '%Y-%m-%dT%H:%M:%S'".format(s)
439
            raise argparse.ArgumentTypeError(msg)
440
        return value
441

  
442
    parser = argparse.ArgumentParser(description='Get profile(s) for NO2')
443
    parser.add_argument('-v', '--verbose', action='count', default=0, help="Be chatty")
444
    parser.add_argument('-i', '--input', metavar='FILE', help='L2 NO2 file for time information')
445
    parser.add_argument('-s', '--species', metavar='MOLECULE', choices=('no2', 'so2', 'ch2o', 'hcho'),
446
                        default='no2', help='Trace gas to extract, default NO2, valid arguments: so2, no2, hcho, ch2o')
447
    parser.add_argument('-t', '--time', metavar='YYYY-MM-DDTHH:MM:SS', type=argparse_datetime,
448
                        help='Date/time to extract for')
449
    parser.add_argument('-a', '--aux', nargs="+", metavar="AUXFILE",
450
                        help="Aux files to search for the correct profiles")
451
    parser.add_argument('--latitude', nargs="+", metavar="LAT", type=float,
452
                        help="Latitude(s) to extract (will be meshed with the longitude argument)")
453
    parser.add_argument('--longitude', nargs="+", metavar="LON", type=float,
454
                        help="Longitude(s) to extract (will be meshed with the latitude argument)")
455
    parser.add_argument('--scanline', nargs="+", metavar="INDEX", type=int,
456
                        help="Scanline(s) to extract")
457
    parser.add_argument('--ground-pixel', nargs="+", metavar="INDEX", type=int,
458
                        help="ground pixel(s) to extract")
459
    parser.add_argument('-A', '--add', action='store_true',
460
                        help="Add profiles to L2 file (requires read/write access to the level 2 file)")
461
    parser.add_argument('-p', '--use-TM5-pressure', action='store_true',
462
                        help="Use the surface pressure from TM5 rather than the L2 surface pressure")
463
    parser.add_argument('-V', '--vmr', action='store_true',
464
                        help="Store profile as vmr rather than partial columns")
465
    parser.add_argument('-q', '--alternative', action="store_false",
466
                        help="Use aternative (direct but slow) interpolation scheme for partial columns.")
467

  
468
    args = parser.parse_args()
469

  
470
    aux_list = args.aux
471
    L2File = args.input
472
    time = args.time
473

  
474
    nox = NO2Profile(aux_list)
475

  
476
    try:
477
        profileFile = nox.select_input_profiles(time=args.time, level2File=args.input)
478
    except RuntimeError:
479
        print("Require one of time or input arguments.")
480
        parser.print_usage()
481
        sys.exit(1)
482

  
483
    if args.verbose > 0:
484
        print("Found profile '{0}', using slice {1}, time {2:%Y-%m-%dT%H:%M:%S}".format(os.path.basename(profileFile), nox.active_slice, nox.time))
485

  
486
    if L2File is not None and args.add and args.alternative:
487
        ref = netCDF4.Dataset(L2File, 'r')
488
        with ref:
489
            lat = ref.groups['PRODUCT'].variables['latitude'][0, :, :]
490
            lon = ref.groups['PRODUCT'].variables['longitude'][0, :, :]
491
            if not args.use_TM5_pressure:
492
                try:
493
                    ps = ref.groups['PRODUCT'].groups['SUPPORT_DATA'].groups['INPUT_DATA'].variables['surface_pressure'][0, :, :]
494
                except:
495
                    ps = None
496
        if args.verbose > 0:
497
            start = datetime.datetime.now()
498
            print("{0:%H:%M:%S.%f}: Starting to extract profiles".format(start))
499

  
500
        nox.extract_vmr(lat, lon, trace_gas=args.species)
501

  
502
        if args.verbose > 0:
503
            end = datetime.datetime.now()
504
            print("{0:%H:%M:%S.%f}: Finished extracting vmr profiles, required {1} seconds".format(end, (end-start).total_seconds()))
505
            start = datetime.datetime.now()
506

  
507
        if not args.vmr:
508
            nox.convert_vmr_to_partial_columns(ps)
509

  
510
            if args.verbose > 0:
511
                end = datetime.datetime.now()
512
                print("{0:%H:%M:%S.%f}: Finished converting vmr profiles, required {1} seconds".format(end, (end-start).total_seconds()))
513
                start = datetime.datetime.now()
514

  
515
        if args.species == 'no2':
516
            output_varname = 'nitrogendioxide_profile_apriori'
517
        elif args.species == 'so2':
518
            output_varname = 'sulfurdioxide_profile_apriori'
519
        else:
520
            output_varname = 'formaldehyde_profile_apriori'
521

  
522
        try:
523
            ref = netCDF4.Dataset(L2File, 'a')
524
        except OSError:
525
            print("Could not open level 2 file in read/write mode to add the profiles")
526
            sys.exit(1)
527

  
528
        with ref:
529
            grp = ref.groups['PRODUCT'].groups['SUPPORT_DATA'].groups['INPUT_DATA']
530
            try:
531
                var = grp.createVariable(output_varname, np.float32,
532
                                dimensions=('time', 'scanline', 'ground_pixel', 'layer'),
533
                                zlib=True, fletcher32=True)
534
            except:
535
                var = grp.variables[output_varname]
536
            var[0, :, :, :] = nox.profiles
537
            var.units = nox.profiles_units
538
            var.standard_name = nox.profiles_standard_name
539
            var.long_name = nox.profiles_long_name
540
            var.coordinates = "/PRODUCT/longitude /PRODUCT/latitude"
541

  
542
        if args.verbose > 0:
543
            end = datetime.datetime.now()
544
            print("{0:%H:%M:%S.%f}: Finished writing profiles, required {1} seconds".format(end, (end-start).total_seconds()))
545

  
546
    elif L2File is not None and args.add:
547
        ref = netCDF4.Dataset(L2File, 'r')
548
        if args.vmr:
549
            print("volume mixing ratio can only be extracted using the 'alternative' method. ")
550
            sys.exit(1)
551

  
552
        with ref:
553
            lat = ref.groups['PRODUCT'].variables['latitude'][0, :, :]
554
            lon = ref.groups['PRODUCT'].variables['longitude'][0, :, :]
555
            if not args.use_TM5_pressure:
556
                try:
557
                    ps = ref.groups['PRODUCT'].groups['SUPPORT_DATA'].groups['INPUT_DATA'].variables['surface_pressure'][0, :, :]
558
                except KeyError:
559
                    ps = None
560
            else:
561
                ps = None
562

  
563
        if args.verbose > 0:
564
            start = datetime.datetime.now()
565
            print("{0:%H:%M:%S.%f}: Starting to extract profiles".format(start))
566

  
567
        nox.extract_partial_columns(lat, lon, ps, trace_gas=args.species)
568

  
569
        if args.verbose > 0:
570
            end = datetime.datetime.now()
571
            print("{0:%H:%M:%S.%f}: Finished extracting profiles, required {1} seconds".format(end, (end-start).total_seconds()))
572
            start = datetime.datetime.now()
573

  
574
        if args.species == 'no2':
575
            output_varname = 'nitrogendioxide_profile'
576
        elif args.species == 'so2':
577
            output_varname = 'sulfurdioxide_profile'
578
        else:
579
            output_varname = 'formaldehyde_profile'
580

  
581
        try:
582
            ref = netCDF4.Dataset(L2File, 'a')
583
        except OSError:
584
            print("Could not open level 2 file in read/write mode to add the profiles")
585
            sys.exit(1)
586

  
587
        with ref:
588
            grp = ref.groups['PRODUCT'].groups['SUPPORT_DATA'].groups['INPUT_DATA']
589
            try:
590
                var = grp.createVariable(output_varname, np.float32,
591
                                dimensions=('time', 'scanline', 'ground_pixel', 'layer'),
592
                                zlib=True, fletcher32=True)
593
            except:
594
                var = grp.variables[output_varname]
595
            var[0, :, :, :] = nox.profiles
596
            var.units = nox.profiles_units
597
            var.standard_name = nox.profiles_standard_name
598
            var.long_name = nox.profiles_long_name
599
            var.coordinates = "/PRODUCT/longitude /PRODUCT/latitude"
600

  
601
        if args.verbose > 0:
602
            end = datetime.datetime.now()
603
            print("{0:%H:%M:%S.%f}: Finished writing profiles, required {1} seconds".format(end, (end-start).total_seconds()))
604
    elif L2File is not None and args.scanline is not None and args.ground_pixel is not None:
605
        ref = netCDF4.Dataset(L2File, 'r')
606
        with ref:
607
            lat = ref.groups['PRODUCT'].variables['latitude'][0, args.scanline, args.ground_pixel]
608
            lon = ref.groups['PRODUCT'].variables['longitude'][0, args.scanline, args.ground_pixel]
609
            l2col = ref.groups['PRODUCT'].groups['SUPPORT_DATA'].groups['DETAILED_RESULTS'].variables['nitrogendioxide_total_column'][0, args.scanline, args.ground_pixel]
610
            psurf = ref.groups['PRODUCT'].groups['SUPPORT_DATA'].groups['INPUT_DATA'].variables['surface_pressure'][0, args.scanline, args.ground_pixel]
611

  
612
        if args.verbose > 0:
613
            print("{0:%H:%M:%S.%f}: Starting to extract profiles".format(datetime.datetime.utcnow()))
614

  
615
        if args.vmr:
616
            nox.extract_vmr(lat, lon, trace_gas=args.species)
617
        else:
618
            nox.extract_partial_columns(lat, lon, psurf, trace_gas=args.species)
619

  
620
        if args.verbose > 0:
621
            print("{0:%H:%M:%S.%f}: Finished extracting profiles".format(datetime.datetime.utcnow()))
622

  
623
        if not args.vmr:
624
            cols = nox.columns
625
        profiles = nox.profiles
626

  
627
        print("Standard name: {0}, units: {1}".format(nox.profiles_standard_name, nox.profiles_units))
628
        for ii in range(l2col.shape[0]):
629
            for jj in range(l2col.shape[1]):
630
                print("Pixel [{0}, {1}] @ ({2:.5f}, {3:.5f})".format(args.scanline[ii], args.ground_pixel[jj], lon[ii,jj], lat[ii,jj]))
631
                if not args.vmr:
632
                    print("    L2 total column: {0:.5e} [{2}], TM5 total column: {1:.5e} [{2}]".format(l2col[ii,jj], cols[ii,jj], nox.columns_units))
633
                print("    L2 surface pressure: {0:.1f} [{1}]".format(psurf[ii,jj], nox.pressure_units))
634
                print("    Profile: [{0}]".format(", ".join(["{0:.3e}".format(v) for v in profiles[ii,jj,:]])))
635
                ps = nox.pressure[ii,jj]
636
                pp = ps * nox.hybm + nox.hyam # mid-level, interfaces are available as well.
637
                print("    pressures: [{0}]".format(", ".join(["{0:.3e}".format(v) for v in pp])))
638
    elif args.latitude is not None and args.longitude is not None:
639
        lat, lon = np.meshgrid(args.latitude, args.longitude)
640
        if args.verbose > 0:
641
            print("{0:%H:%M:%S.%f}: Starting to extract profiles".format(datetime.datetime.utcnow()))
642

  
643
        if args.vmr:
644
            nox.extract_vmr(lat, lon, trace_gas=args.species)
645
        else:
646
            nox.extract_partial_columns(lat, lon, trace_gas=args.species)
647

  
648
        if args.verbose > 0:
649
            print("{0:%H:%M:%S.%f}: Finished extracting profiles".format(datetime.datetime.utcnow()))
650

  
651
        if not args.vmr:
652
            cols = nox.columns
653
        profiles = nox.profiles
654

  
655
        print("Standard name: {0}, units: {1}".format(nox.profiles_standard_name, nox.profiles_units))
656
        for ii, lat in enumerate(args.latitude):
657
            lon = args.longitude[ii]
658
            print("Pixel [{0}] @ ({1:.5f}, {2:.5f})".format(ii, lon, lat))
659
            if not args.vmr:
660
                print("    TM5 total column: {0:.5e} [{1}]".format(cols[ii], nox.columns_units))
661
            print("    TM5 surface pressure: {0:.1f} [{1}]".format(nox.pressure[ii], nox.pressure_units))
662
            print("    Profile: [{0}]".format(", ".join(["{0:.3e}".format(v) for v in profiles[ii,:]])))
663
            ps = nox.pressure[ii]
664
            pp = ps * nox.hybm + nox.hyam # mid-level, interfaces are available as well.
665
            print("    pressures: [{0}]".format(", ".join(["{0:.3e}".format(v) for v in pp])))
666
    else:
667
        print("The combination of parameters is illogical of inconsistent")
668
        parser.print_usage()
669
        sys.exit(0)

Also available in: Unified diff