Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / TM5_profile.py @ 815:c0e5b3712092

History | View | Annotate | Download (23.3 KB)

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()