import numpy as np
import xarray as xr
import netCDF4 as nc
import datetime as dtm
from .....utils.check.errclass import CifError, CifNotImplementedError, CifValueError
[docs]
def read(
self,
name,
varnames,
dates,
files,
interpol_flx=False,
comp_type=None,
ddi=None,
**kwargs
):
"""Get initial and boundary conditions from pre-computed WRF-Chem
files and load them into a pyCIF variable.
Reading initial conditions (``comp_type="inicond"``) is not
implemented yet and always raises. Reading lateral conditions
(``comp_type="latcond"``) compares the requested timestamps to those
found in the ``wrfbdy`` file: either they match exactly, or exactly
one extra timestamp is requested beyond the file's last one, in
which case its value is extrapolated from the ``_BTxx`` trend
variables. The four boundary-side variables
(``_BXS``/``_BXE``/``_BYS``/``_BYE``) are then read and assembled
into a single perimeter dimension matching
``self.domain.zlon_side``.
Args:
self: the wrfchem_icbc plugin
name: the name of the component
varnames: unused, accepted for interface compatibility
dates: list of ``[start, end]`` date pairs to extract
files: list with the single ``wrfbdy``/``wrfinput`` file path
interpol_flx (bool): unused, accepted for interface compatibility
comp_type: ``'inicond'`` or ``'latcond'``
ddi: unused, accepted for interface compatibility
Returns:
xarray.DataArray with dims ``(time, lev, lat, lon)``, where
``lon`` spans the domain perimeter.
Raises:
CifError: if ``comp_type`` is ``None`` or not recognized.
CifNotImplementedError: if ``comp_type == "inicond"``.
CifValueError: if more than one file is given for ``latcond``,
or if the file's timestamps don't match the expected
pattern.
"""
# Check the type of limit condition to check
if comp_type is None:
raise CifError(
"Trying to read limit conditions for wrfchem, "
"but did not specify the type"
)
# Read initial conditions
if comp_type == "inicond":
raise CifNotImplementedError("To do: implement 'read' for inicond")
# ic_file = files[0]
# ds = xr.open_dataset(ic_file)
# data = ds[var2extract][:]
#
# # -- Add fake lat dimensions at the end (TODO: remove 1e9)
# # data = data.values[:, :, np.newaxis, :] / 1e9
# data = data.values[:, :, np.newaxis, :]
#
# xmod = xr.DataArray(
# data,
# coords={"time": [np.min(dates)]},
# dims=("time", "lev", "lat", "lon"),
# )
# Read Lateral boundary conditions
elif comp_type == "latcond":
if len(files)>1:
raise CifValueError("There should only be one latcond file")
ncf = nc.Dataset(files[0])
# "datef" requested? If not in file (which is how I set
# up the wrfchem model mapper), have to compute from bt.
# Unpack requested times
times_req = sorted(set([x for y in dates for x in y]))
# Get times in file
times_char = ncf["Times"][:].astype(str)
times_str = np.apply_along_axis(lambda x: "".join(x), 1, times_char)
times_file = [dtm.datetime.strptime(t, "%Y-%m-%d_%H:%M:%S") for t in times_str]
# Compare requested times to those in file
# Can do the following cases:
# Case 1: All requested times are in the file
if set(times_req)==set(times_file):
ind_copy = range(len(times_file))
tend = False
# Case 2: Requesting one more date than in file
# -> Compute from tendency variables
# (This should be the usual case)
elif (set(times_req[:-1])==set(times_file)) & \
((times_file[-1] + (times_file[-1]-times_file[-2]))==times_req[-1]):
ind_copy = range(len(times_file))
tend = True
dt = (times_file[-1]-times_file[-2]).total_seconds()
else:
raise CifValueError("Times available in latcond file " + \
"are not as expected")
varname = getattr(self.chemistry.acspecies, name).varname
# Read all values
b = dict(
XS = ncf[varname + "_BXS"][:],
XE = ncf[varname + "_BXE"][:],
YS = ncf[varname + "_BYS"][:],
YE = ncf[varname + "_BYE"][:])
if tend:
btf = dict(
XS = ncf[varname + "_BTXS"][-1, ...],
XE = ncf[varname + "_BTXE"][-1, ...],
YS = ncf[varname + "_BTYS"][-1, ...],
YE = ncf[varname + "_BTYE"][-1, ...])
bdy_width = len(ncf.dimensions["bdy_width"])
ncf.close()
# Arrange the data into a single variable as defined
# in domains/wrfchem/get_sides.py (this is much the
# same as in write.py).
icbc = np.zeros((len(dates), ) + \
(len(self.domain.sigma_a)-1, ) + \
(self.domain.zlon_side.shape))
icbc = icbc*np.nan
nlon2d = self.domain.zlon2d.shape[1]
nlat2d = self.domain.zlon2d.shape[0]
# Fill variable
indi = 0
for nw in range(bdy_width):
inde = indi + nlon2d
icbc[ind_copy,:,0,indi:inde] = b["YS"][:,nw,:,:]
if tend:
icbc[-1, :,0,indi:inde] = b["YS"][-1,nw,:,:] + \
dt*btf["YS"][nw,:,:]
indi = inde
inde = indi + nlon2d
icbc[ind_copy,:,0,indi:inde] = b["YE"][:,nw,:,:]
if tend:
icbc[-1,:,0,indi:inde] = b["YE"][-1,nw,:,:] + \
dt*btf["YE"][nw,:,:]
indi = inde
inde = indi + nlat2d
icbc[ind_copy,:,0,indi:inde] = b["XS"][:,nw,:,:]
if tend:
icbc[-1,:,0,indi:inde] = b["XS"][-1,nw,:,:] + \
dt*btf["XS"][nw,:,:]
indi = inde
inde = indi + nlat2d
icbc[ind_copy,:,0,indi:inde] = b["XE"][:,nw,:,:]
if tend:
icbc[-1,:,0,indi:inde] = b["XE"][-1,nw,:,:] + \
dt*btf["XE"][nw,:,:]
indi = inde
# Put the data into an xarray
xmod = xr.DataArray(
icbc, coords={"time": times_req}, dims=("time", "lev", "lat", "lon")
)
else:
raise CifError(
f"Could not recognize the type of file to read in wrfchem: {comp_type}"
)
return xmod
