source: palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_camx.py @ 4866

Last change on this file since 4866 was 4843, checked in by raasch, 4 years ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

File size: 7.4 KB
Line 
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3
4#------------------------------------------------------------------------------#
5#
6# Scripts for processing of WRF and CAMx files to PALM dynamic driver
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the GNU General Public License
19# along with this program.  If not, see <https://www.gnu.org/licenses/>.
20#
21# Copyright 2018-2021 Institute of Computer Science
22#                     of the Czech Academy of Sciences, Prague
23# Authors: Krystof Eben, Jaroslav Resler, Pavel Krc
24#
25#------------------------------------------------------------------------------#
26'''
27Load and process initial and boundary conditions from CAMx NetCDF files
28'''
29
30import sys
31import datetime
32import re
33import numpy as np
34import netCDF4
35from metpy.interpolate import interpolate_1d
36
37import palm_wrf_utils
38
39_na = np.newaxis
40re_num = re.compile(r'[0-9\.]+')
41
42class HelperNotFound(Exception):
43    def __init__(self, name):
44        self.name = name
45
46class Helpers(object):
47    '''Pre-loaded helper variables from CAMx'''
48    def __getattr__(self, name):
49        try:
50            return self.__dict__[name]
51        except KeyError:
52            raise HelperNotFound(name)
53
54def tflag(data, req_dts):
55    assert len(data.shape) == 3 and data.shape[2] == 2
56    xdate = data[:,0,0]
57    xtime = data[:,0,1]
58
59    # Verify that dates are equal for each variable
60    assert (data[:,:,0] == xdate[:,_na]).all()
61    assert (data[:,:,1] == xtime[:,_na]).all()
62
63    dts = []
64    for i in range(len(xdate)):
65        dt = datetime.datetime.strptime(
66            '{0:07d} {1:06d}'.format(xdate[i], xtime[i]), '%Y%j %H%M%S')
67        try:
68            ireq = req_dts[dt]
69        except KeyError:
70            continue
71        dts.append((ireq, i))
72    return dts
73
74def load_conversion(spc, f, sl, hlp, camx_vars, camx_units, formula, **kwargs):
75    loaded_vars = []
76    for varname, unit in zip(camx_vars, camx_units):
77        try:
78            v = f.variables[varname]
79        except KeyError:
80            print('Skipping {0}, because input variable {1}({2}) is missing.'.format(
81                spc, varname, unit))
82            return None
83        if getattr(v, 'units', None) != unit:
84            print('Skipping {0}, because input variable {1} has wrong unit ({2} <> {3}).'.format(
85                spc, varname, getattr(v, 'units', None), unit))
86            return None
87
88        print('Loading variable {0}.'.format(varname))
89        loaded_vars.append(v[sl])
90    try:
91        val = formula(*(loaded_vars + [hlp]))
92    except HelperNotFound as e:
93        print('Skipping {0} - missing helper variable {1}.'.format(spc, e.name))
94        return None
95
96    return val
97
98def process_tstep(f, itf, regridder, lay_height, fout, itout, z_levels,
99        vars_remaining, filled, conversions, helpers):
100
101    # Load helper vars for this timestep
102    hlp = Helpers()
103    for helper_name, helper in helpers:
104        data = load_conversion(helper_name, f,
105                (itf,slice(None),regridder.ys,regridder.xs), hlp, **helper)
106        if data is not None:
107            setattr(hlp, helper_name, data)
108
109    # Load all usable vars for this timestep, regrid horizontally
110    varmeta = []
111    vardata = []
112    for spc in list(vars_remaining):
113        conv = conversions[spc]
114        data = load_conversion(spc, f, (itf,slice(None),regridder.ys,regridder.xs),
115                hlp, **conv)
116        if data is None:
117            continue
118
119        data = regridder.regrid(data)
120        vardata.append(np.r_[data[0:1], data]) #add peg below
121        varmeta.append((spc, conv['output_unit']))
122
123    # Perform vertical interpolation on all currently loaded vars at once
124    print('Interpolating vertically...')
125    vinterp = interpolate_1d(z_levels, lay_height, *vardata)
126    if len(vardata) == 1:
127        # return_list_always=True argument is only in later versions of MetPy
128        vinterp = [vinterp]
129    del vardata
130    for (vn, vu), vd in zip(varmeta, vinterp):
131        v = fout.variables[vn]
132        v[itout] = vd
133        v.units = vu
134        filled[vn][itout] = True
135
136def process_files(camx_file_list, camx_interp_fname, palm_grid_lat,
137        palm_grid_lon, terrain_rel, z_levels, times, species_names,
138        conversions, helpers):
139
140    terrain_shift = terrain_rel[_na,:,:]
141    lowest_layer = np.zeros(((1,) + palm_grid_lat.shape), dtype='f4')
142    lowest_layer[:] = -999.
143    tindex = dict((dt, i) for i, dt in enumerate(times))
144    filled = {}
145
146    with netCDF4.Dataset(camx_interp_fname, 'w', format='NETCDF4') as fout:
147        fout.createDimension('time', len(times))
148        fout.createDimension('z', len(z_levels))
149        fout.createDimension('y', palm_grid_lat.shape[0])
150        fout.createDimension('x', palm_grid_lat.shape[1])
151        for vn in species_names:
152            fout.createVariable(vn, 'f4', ('time', 'z', 'y', 'x'))
153            filled[vn] = [False] * len(times)
154
155        for fname in sorted(camx_file_list):
156            with netCDF4.Dataset(fname) as f:
157                dts = tflag(f.variables['TFLAG'][:], tindex)
158                if dts:
159                    print('Processing CAMx file {0}.'.format(fname))
160
161                    # preprare projection
162                    trans = palm_wrf_utils.CAMxCoordTransform(f)
163                    palm_in_camx_y, palm_in_camx_x = trans.latlon_to_ji(
164                                                    palm_grid_lat, palm_grid_lon)
165                    regridder = palm_wrf_utils.BilinearRegridder(
166                            palm_in_camx_x, palm_in_camx_y, preloaded=True)
167
168                    # locate layer heights
169                    try:
170                        vz = f.variables['z']
171                    except KeyError:
172                        print('Loading heights from separate file')
173                        with open(fname+'.heights') as fh:
174                            fix_hgt = np.array(list(map(float, re_num.findall(fh.read())))) * 1000. #orig in km
175                            fix_hgt = fix_hgt[:,_na,_na]
176                    else:
177                        print('Loading heights from variable z')
178                        fix_hgt = None
179
180                    for itout, itf in dts:
181                        print('Timestep {0}'.format(itout))
182                        vars_remaining = [vn for vn, vf in filled.items() if not vf[itout]]
183
184                        lay_height = fix_hgt if fix_hgt is not None else (
185                                regridder.regrid(vz[itf,:,regridder.ys,regridder.xs]))
186                        lay_height = np.r_[lowest_layer, lay_height + terrain_shift]
187                                                #add 1 pegging layer always below
188                        process_tstep(f, itf, regridder, lay_height, fout, itout, z_levels,
189                                vars_remaining, filled, conversions, helpers)
190                else:
191                    print('Skipping CAMx file {0} - no required times.'.format(fname))
192
193    if not all(all(vf) for vf in filled.values()):
194        sys.exit('CAMx data not complete - missing some variables/timesteps: {0}'
195                .format(filled))
Note: See TracBrowser for help on using the repository browser.