Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / TM5_profile.py @ 929:0ad059426223

History | View | Annotate | Download (24.7 KB)

1
#!/usr/bin/env python3
2

    
3
# Extract trace gas profile(s) from a 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
__version__ = "1.0.1"
14
__date__ = "February 2022" 
15
DownloadURL = "http://www.tropomi.eu/data-products/nitrogen-dioxide/"
16
 
17
# =====================================        
18

    
19
import os
20
import sys
21
import argparse
22
import datetime
23
from requests.structures import CaseInsensitiveDict
24

    
25
import warnings
26
warnings.filterwarnings("ignore", category=DeprecationWarning)
27

    
28
import inspect
29

    
30
import numpy as np
31
import netCDF4
32
from scipy import constants
33

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

    
51
# =====================================        
52

    
53
class CTM:
54
    def __init__(self, fname:str, species:str='no2'):
55
        self.fname = fname
56
        self.slice_index = None
57
        self.weights = None
58
        self.spatial_index_lon = None
59
        self.spatial_index_lat = None
60
        self.species = species
61
        self.final_output_shape = None
62
        self.ref = None
63
    
64
    def __enter__(self): 
65
        self.ref = netCDF4.Dataset(self.fname, 'r')
66
        return self.ref
67
  
68
    def __exit__(self, exception_type, exception_value, traceback): 
69
        if self.ref is not None:
70
            self.ref.close()
71
            self.ref = None
72
    
73
    def close(self):
74
        self.__exit__(None, None, None)
75
    
76
    def find_slice(self, t: datetime.datetime) -> int:
77
        """Find the correct time index"""
78
        time_ref = self.ref["time"]
79
        times = time_ref[:]
80
        time_units = time_ref.units
81
        
82
        try:
83
            f_times = netCDF4.num2date(times, time_units, 
84
                        only_use_cftime_datetimes=False, 
85
                        only_use_python_datetimes=True)
86
        except TypeError:
87
            f_times = netCDF4.num2date(times, time_units)
88
        
89
        quarter = datetime.timedelta(minutes=15)
90
        second  = datetime.timedelta(seconds=1)
91
        extra = quarter - second
92

    
93
        if t < f_times[0] - extra or t > f_times[-1] + extra:
94
            raise RuntimeError("Requested time not covered in CTM file.")
95
        else:
96
            self.slice_index = np.argmin([np.fabs((t - tt).total_seconds()) for tt in f_times])
97
        self.slice_time = f_times[self.slice_index]
98
        return self.slice_index,self.slice_time
99
        
100
    def find_indices_weights(self, latitude, longitude):
101
        """Find indices and weights for the spatial interpolation
102

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

    
167
        return ix4a, iy4a, w4a
168
    
169
    def get_pressures(self, surface_pressure):
170
        hyai = self.ref['hyai'][:]
171
        hybi = self.ref['hybi'][:]
172
        return hyai + hyai*surface_pressure
173
        
174
    def get_pressures_interfaces(self):
175
        """Get pressure at the interfaces of the model"""
176
        field = self.ref[self.species]
177
        pressure_shape = list(field.shape[1:])
178
        pressure_shape[0] += 1
179
        p_interfaces = np.zeros(pressure_shape, dtype=np.float32)
180
        ps = self.ref['ps'][self.slice_index, ...]
181
        hyai = self.ref['hyai'][:]
182
        hybi = self.ref['hybi'][:]
183
        
184
        for i, (a, b) in enumerate(zip(hyai, hybi)):
185
            p_interfaces[i, ...] = a + b*ps
186
            
187
        return p_interfaces
188

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

    
320

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

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

    
471
def geospec_verify(s: str):
472
    global named_geolocation_points
473
    if s not in named_geolocation_points:
474
        names = "', '".join(sorted(list(named_geolocation_points.keys())))
475
        msg = "'{0}' not found. \nValid point names are '{1}'.".format(s, names)
476
        raise argparse.ArgumentTypeError(msg)
477
    return s
478

    
479
def geo_point_name(name: str):
480
    global named_geolocation_points
481
    # https://mynasadata.larc.nasa.gov/latitudelongitude-finder/
482
    if name in named_geolocation_points:
483
        return (named_geolocation_points[name]['lon'], named_geolocation_points[name]['lat'])
484
    else:
485
        names = "', '".join(sorted(list(named_geolocation_points.keys())))
486
        raise RuntimeError("Warning, '{0}' not found. \nValid point names are '{1}'.".format(name, names))
487

    
488
# =====================================        
489

    
490
def main():
491
    global named_geolocation_points
492
    geolocation_point_name_list = list(named_geolocation_points.keys())
493
    geolocation_point_name_list.sort()
494
    geolocation_point_names = "', '".join(geolocation_point_name_list)
495
    
496
    parser = argparse.ArgumentParser(description='Get CTM profile(s) for NO2, SO2 or HCHO.',
497
                                     epilog="This program can be downloaded from {0}".format(DownloadURL))
498
    
499
    # Main parser arguments
500
    
501
    parser.add_argument('--version', action='version', version='%(prog)s {0} ({1})'.format(__version__, __date__), 
502
                        help='show version info and exit')
503
    
504
    parser.add_argument('-v', '--verbose', action='count', default=0, help="be chatty")
505
    
506
    parser.add_argument('-s', '--species', metavar='MOLECULE', 
507
                        choices=('no2', 'so2', 'ch2o', 'hcho', 'NO2', 'SO2', 'CH2O', 'HCHO', 'NO???', 'SO???', 'CH???O'),
508
                        default='no2', 
509
                        help='trace gas for which to extract the profil, default: NO2; valid arguments: NO2, SO2, HCHO, CH2O')
510
    
511
    parser.add_argument('-c', '--ctm', required=True, help='the CTM file with the profiles')
512
    
513
    parser.add_argument('-V', '--vmr', action='store_true',
514
                        help="store the profile as vmr rather than partial columns")
515
    
516
    subparsers = parser.add_subparsers(help='Processing mode (these have their own help)', required=True)
517
    
518
    # Subparser arguments for "single" mode
519
     
520
    parser_s = subparsers.add_parser('single', help='Extract the profile for a single time and location.')
521
    
522
    parser_s.set_defaults(func=extract_single)
523
    parser_s.add_argument('-t', '--time', metavar='YYYY-MM-DDTHH:MM:SS', type=argparse_datetime,
524
                          required=True,
525
                          help='date and time to extract the profile for')
526
    group = parser_s.add_mutually_exclusive_group(required=True)
527
    group.add_argument('-g', '--geo-point', metavar="L", type=float, dest="geo_point", 
528
                          nargs=2, 
529
                          help="location to extract the profile for as a LAT LON pair")
530
    group.add_argument('-G', '--geo-name', metavar="NAME", type=geospec_verify, dest="geo_point_name", 
531
                          help="location to extract a profile for by name; recognized names are: '{0}'.".format(geolocation_point_names))
532
    parser_s.add_argument('-p', '--pressure', type=float, default=None, 
533
                          help="surface pressure in hPa to use for the profile; if not given the TM5 surface pressure is used")
534
    
535
    # Subparser arguments for "level2" mode
536
    
537
    parser_m = subparsers.add_parser('level2', help='Extract profiles for a full level-2 file; the profiles will be added to the file (for this mode read/write access to the level-2 file is required).')
538
    
539
    parser_m.set_defaults(func=extract_orbit)
540
    parser_m.add_argument('-i', '--input', metavar='L2_FILE', dest='l2_file', required=True, help='level-2 file to add the profile data to')
541
    parser_m.add_argument('--use-TM5-pressure', action='store_true',
542
                          help="use the surface pressure from TM5 rather than the surface pressure given in the level-2 file")
543
    parser_m.add_argument('--add-temperature-profile', action='store_true',
544
                          help="also add the temperature profile from TM5 to the level-2 file")
545
    
546
    # Check the parser argements
547
    
548
    args = vars(parser.parse_args())
549
    
550
    processor = args["func"]
551
    del args["func"]
552
        
553
    if "geo_point_name" in args and args["geo_point_name"]:
554
        args["geo_point"] = geo_point_name(args["geo_point_name"])
555
        del args["geo_point_name"]
556
    elif "geo_point" in args and args["geo_point"]:
557
        args["geo_point"] = args["geo_point"][1], args["geo_point"][0]
558
    
559
    if args["species"] in ("NO2", "no2", "NO???"):
560
        args["species"] = "no2"
561
    elif args["species"] in ("SO2", "so2", "SO???"):
562
        args["species"] = "so2"
563
    else:
564
        args["species"] = "ch2o"
565
    
566
    # Call the processor
567
    
568
    processor(**args)
569
    
570
    
571
if __name__ == "__main__":
572
    main()
573