import numpy as np
from logging import info
[docs]
def vcoord(obsvect, **kwargs):
"""Assign a model vertical level to each observation.
Dispatches to one of two strategies based on whether a pre-defined
station level file is available:
* :func:`vcoordfromfile` — when ``obsvect.file_statlev`` is set, reads
a text table mapping station names to fixed level indices.
* :func:`vcoordfrommeteo` — otherwise, derives the level by matching
each observation's altitude against a hard-coded mid-layer altitude
array (LMDZ-specific, provisional implementation).
The ``level`` column of ``obsvect.datastore["data"]`` is filled in-place.
Args:
obsvect (Plugin): obsvect plugin instance (carries ``datastore``,
optionally ``file_statlev``, and ``workdir``).
**kwargs: forwarded to the chosen sub-function.
Returns:
Plugin: the updated *obsvect* (also modified in-place).
"""
info("Finding model levels corresponding to observations")
# Don't do anything if the datastore is empty
if len(obsvect.datastore["data"]) == 0:
return obsvect
# If a file with fixed vertical coordinates is specified, use it
if hasattr(obsvect, "file_statlev"):
file_statlev = obsvect.file_statlev
info(
f"Using pre-defined vertical coordinates for stations: {file_statlev}"
)
obsvect.datastore["data"] = vcoordfromfile(
obsvect.datastore["data"], file_statlev, **kwargs
)
# Else compute vertical coordinates from meteo files
# To be coded from pyCIF-CHIMERE
# To be generalized by using meteo plugin
else:
obsvect.datastore["data"] = vcoordfrommeteo(
obsvect.workdir, obsvect.datastore["data"], **kwargs
)
return obsvect
[docs]
def vcoordfromfile(datastore, file_lev, **kwargs):
"""Assign model levels to observations using a pre-defined station table.
Reads a whitespace-delimited file (header line skipped) whose columns are:
.. code-block:: text
STAT LAT LON LEV alt(m)
and sets ``datastore["level"]`` to the model level index (column 4,
0-based Python indexing) for each matching station name. Observations
from stations not listed in the file retain ``NaN``.
Args:
datastore (pd.DataFrame): observation dataframe with a ``"station"``
column (station names, lower-case).
file_lev (str): path to the station-level file.
**kwargs: unused; accepted for interface consistency.
Returns:
pd.DataFrame: *datastore* with the ``"level"`` column updated
in-place.
"""
lev_stats = np.genfromtxt(file_lev, skip_header=1, usecols=0, dtype=str)
lev_infos = np.genfromtxt(file_lev, skip_header=1)[:, 1:]
datastore.loc[:, "level"] = np.nan
for s, linfo in zip(lev_stats, lev_infos):
mask = datastore["station"] == s.lower()
datastore.loc[mask, "level"] = linfo[2]
return datastore
[docs]
def vcoordfrommeteo(workdir, datastore, **kwargs):
"""Assign model levels to observations by matching altitudes to LMDZ mid-layers.
Uses a hard-coded array of LMDZ mid-layer altitudes (39 levels, metres
above the surface) and assigns each observation to the nearest level by
minimising :math:`|\\text{alt}_{obs} - \\text{alt}_{level}|`.
.. warning::
This function is provisional and LMDZ-specific. The altitude array
is hard-coded and the meteo directory (``$workdir/meteo/``) is not
actually used in the current implementation. A commented-out block
shows the intended xarray-based generalisation.
Args:
workdir (str): pyCIF working directory (currently unused; reserved
for the future generalised implementation).
datastore (pd.DataFrame): observation dataframe with an ``"alt"``
column (altitude in metres above the surface).
**kwargs: unused.
Returns:
pd.DataFrame: *datastore* with the ``"level"`` column filled with
1-based LMDZ level indices (closest level wins).
"""
meteodir = workdir + "/meteo/"
datastore["level"] = np.nan
# ds = xr.open_dataset(meteodir + 'fluxstoke.an2012.m01.nc')
# alt = ds['phi'] / 9.81
# alt = alt.mean(dim=['time_counter', 'lat', 'lon'])
alt = np.array(
[
419.947858,
491.729716,
590.03352,
731.746909,
933.003129,
1209.481853,
1577.072424,
2052.229597,
2652.034343,
3392.864731,
4287.454615,
5338.23227,
6527.101582,
7802.693814,
9078.312502,
10258.393071,
11298.825452,
12251.020931,
13215.014871,
14244.018004,
15347.801437,
16530.826455,
17803.910146,
19179.506405,
20667.282816,
22274.191707,
24012.661236,
25898.57798,
27949.867586,
30191.221244,
32660.342496,
35402.092974,
38490.531506,
41997.093232,
46044.013939,
50511.472241,
55800.118528,
61192.883408,
72735.110678,
]
)
for index, row in datastore.iterrows():
idx = (np.abs(alt - row["alt"])).argmin()
datastore.at[index, "level"] = idx + 1
return datastore
