diff --git a/cloudnetpy/categorize/model.py b/cloudnetpy/categorize/model.py index 805421b8..0387f9f7 100644 --- a/cloudnetpy/categorize/model.py +++ b/cloudnetpy/categorize/model.py @@ -11,7 +11,6 @@ from cloudnetpy import utils from cloudnetpy.cloudnetarray import CloudnetArray -from cloudnetpy.constants import G from cloudnetpy.datasource import DataSource from cloudnetpy.exceptions import ModelDataError @@ -146,27 +145,9 @@ def _get_model_heights(self, alt_site: float) -> npt.NDArray: except KeyError as err: msg = "No 'height' variable in the model file." raise ModelDataError(msg) from err - surface_altitude = self._get_model_surface_altitude(alt_site) + surface_altitude = utils.get_model_surface_altitude(self.dataset, alt_site) return self.to_m(model_heights) + surface_altitude - def _get_model_surface_altitude(self, alt_site: float) -> float | npt.NDArray: - """Returns model surface altitude from geopotential if available. - - For sites in complex terrain (e.g. mountains), the model grid cell - surface height can differ significantly from the actual site altitude. - Using the model's own surface height ensures that thermodynamic fields - are placed at their physically correct absolute heights. - - Note: Model surface altitude might be higher than the site altitude. - """ - if "sfc_height" in self.dataset.variables: - sfc_height = self.dataset.variables["sfc_height"][:] - return sfc_height[:, np.newaxis] - if "sfc_geopotential" in self.dataset.variables: - geopotential = self.dataset.variables["sfc_geopotential"][:] - return geopotential[:, np.newaxis] / G - return alt_site - def calc_attenuations(self, frequency: float) -> None: temperature = self.getvar("temperature") pressure = self.getvar("pressure") diff --git a/cloudnetpy/cli.py b/cloudnetpy/cli.py index 8e485b9e..eecf282a 100644 --- a/cloudnetpy/cli.py +++ b/cloudnetpy/cli.py @@ -90,6 +90,33 @@ def run(args: argparse.Namespace, tmpdir: str, client: APIClient) -> None: l2_filename = _process_mwrpy_product(product, mwrpy_filepath, args) _plot(l2_filename, product, args) + # Radar + model based products (e.g. epsilon-radar) + if "epsilon-radar" in args.products: + epsilon_filepath = _process_epsilon_radar(cat_files, args, client) + if epsilon_filepath is not None: + _plot(epsilon_filepath, "epsilon-radar", args) + + +def _process_epsilon_radar( + cat_files: dict, args: argparse.Namespace, client: APIClient +) -> str | None: + radar_filepath = cat_files.get("radar") or _fetch_product(args, "radar", client) + if radar_filepath is None: + logging.info("No radar data available for epsilon-radar") + return None + model_filepath = _fetch_model(args, client) + if model_filepath is None: + logging.info("No model data available for epsilon-radar") + return None + filename = f"{args.date.replace('-', '')}_{args.site}_epsilon-radar.nc" + output_file = _create_output_folder("geophysical", args) / filename + products = importlib.import_module("cloudnetpy.products") + products.generate_epsilon_from_radar( + str(radar_filepath), model_filepath, str(output_file) + ) + logging.info("Processed epsilon-radar: %s", output_file) + return str(output_file) + def _process_categorize( input_files: dict, diff --git a/cloudnetpy/output.py b/cloudnetpy/output.py index e203d877..b032da38 100644 --- a/cloudnetpy/output.py +++ b/cloudnetpy/output.py @@ -74,6 +74,7 @@ def save_product_file( file_name: str | PathLike, uuid: UUID, copy_from_cat: tuple = (), + extra_sources: tuple[DataSource, ...] = (), ) -> None: """Saves a standard Cloudnet product file. @@ -83,13 +84,18 @@ def save_product_file( file_name: Name of the output file to be generated. uuid: Set specific UUID for the file. copy_from_cat: Variables to be copied from the categorize file. + extra_sources: Additional input DataSources whose ``file_uuid`` should + also be listed in ``source_file_uuids`` (e.g. a model file used + alongside the primary L1b input). """ human_readable_file_type = _get_identifier(short_id) - dimensions = { - "time": len(obj.time), - "height": len(obj.dataset.variables["height"]), - } + height_size = ( + len(obj.data["height"][:]) + if "height" in obj.data + else len(obj.dataset.variables["height"]) + ) + dimensions = {"time": len(obj.time), "height": height_size} with init_file(file_name, dimensions, obj.data, uuid) as nc: nc.cloudnet_file_type = short_id vars_from_source = ( @@ -105,7 +111,7 @@ def save_product_file( f"{human_readable_file_type.capitalize()} products from" f" {obj.dataset.location}" ) - nc.source_file_uuids = get_source_uuids([nc, obj]) + nc.source_file_uuids = get_source_uuids([nc, obj, *extra_sources]) copy_global( obj.dataset, nc, @@ -116,13 +122,28 @@ def save_product_file( "year", "source", "source_instrument_pids", + "instrument_pid", "voodoonet_version", ), ) - merge_history(nc, human_readable_file_type, obj) + _append_extra_sources(nc, extra_sources) + merge_history(nc, human_readable_file_type, obj, extra_sources=extra_sources) nc.references = get_references(short_id) +def _append_extra_sources( + nc: netCDF4.Dataset, extra_sources: tuple[DataSource, ...] +) -> None: + """Merges ``source`` strings from auxiliary input files into ``nc.source``.""" + if not extra_sources: + return + existing = nc.source.split("\n") if "source" in nc.ncattrs() else [] + extras = [src.dataset.source for src in extra_sources if src.dataset.source] + merged = list(dict.fromkeys([*existing, *extras])) + if merged: + nc.source = "\n".join(merged) + + def get_l1b_source(instrument: Instrument) -> str: """Returns level 1b file source.""" parts = [ @@ -173,6 +194,11 @@ def get_references(identifier: str | None = None, extra: list | None = None) -> references += ", https://doi.org/10.1175/JAM2340.1" case "drizzle": references += ", https://doi.org/10.1175/JAM-2181.1" + case "epsilon-radar": + references += ( + ", https://doi.org/10.5194/amt-13-5335-2020" + ", https://doi.org/10.1002/2015JD024543" + ) if extra is not None: for reference in extra: references += f", {reference}" @@ -203,14 +229,21 @@ def get_source_uuids(data: Observations | list[netCDF4.Dataset | DataSource]) -> def merge_history( - nc: netCDF4.Dataset, file_type: str, data: Observations | DataSource + nc: netCDF4.Dataset, + file_type: str, + data: Observations | DataSource, + extra_sources: tuple[DataSource, ...] = (), ) -> None: """Merges history fields from one or several files and creates a new record.""" def extract_history(obj: DataSource | Observations) -> list[str]: if hasattr(obj, "dataset") and hasattr(obj.dataset, "history"): history = obj.dataset.history - if isinstance(obj, Model): + is_model_file = ( + isinstance(obj, Model) + or getattr(obj.dataset, "cloudnet_file_type", "") == "model" + ) + if is_model_file: return [history.split("\n")[-1]] return history.split("\n") return [] @@ -221,6 +254,8 @@ def extract_history(obj: DataSource | Observations) -> list[str]: elif isinstance(data, Observations): for field in fields(data): histories.extend(extract_history(getattr(data, field.name))) + for src in extra_sources: + histories.extend(extract_history(src)) # Remove duplicates histories = list(dict.fromkeys(histories)) @@ -298,7 +333,7 @@ def copy_variables( """ for key in keys: - if key in source.variables: + if key in source.variables and key not in target.variables: fill_value = getattr(source.variables[key], "_FillValue", False) variable = source.variables[key] var_out = target.createVariable( @@ -436,6 +471,7 @@ def _get_identifier(short_id: str) -> str: "der", "ier", "classification-voodoo", + "epsilon-radar", ) if short_id not in valid_ids: msg = f"Invalid file identifier: {short_id}" @@ -448,6 +484,8 @@ def _get_identifier(short_id: str) -> str: return "ice effective radius" if short_id == "der": return "droplet effective radius" + if short_id == "epsilon-radar": + return "dissipation rate of turbulent kinetic energy" return short_id diff --git a/cloudnetpy/plotting/plot_meta.py b/cloudnetpy/plotting/plot_meta.py index 56234d35..375687a4 100644 --- a/cloudnetpy/plotting/plot_meta.py +++ b/cloudnetpy/plotting/plot_meta.py @@ -667,5 +667,10 @@ class PlotMeta(NamedTuple): plot_range=(1e-7, 1e-1), log_scale=True, ), + "epsilon_error": PlotMeta( + cmap="inferno", + plot_range=(1e-7, 1e-1), + log_scale=True, + ), }, } diff --git a/cloudnetpy/products/__init__.py b/cloudnetpy/products/__init__.py index 3d742d04..65cf087f 100644 --- a/cloudnetpy/products/__init__.py +++ b/cloudnetpy/products/__init__.py @@ -1,6 +1,7 @@ from .classification import generate_classification from .der import generate_der from .drizzle import generate_drizzle +from .epsilon_radar import generate_epsilon_from_radar from .ier import generate_ier from .iwc import generate_iwc from .lwc import generate_lwc diff --git a/cloudnetpy/products/epsilon_radar.py b/cloudnetpy/products/epsilon_radar.py new file mode 100644 index 00000000..04fd8be8 --- /dev/null +++ b/cloudnetpy/products/epsilon_radar.py @@ -0,0 +1,373 @@ +"""Module for creating Cloudnet eddy dissipation rate product, based on the +pipeline of Griesche et al. (2020) with the inertial-subrange slope-acceptance +criterion of Borque et al. (2016). +""" + +from collections.abc import Callable +from os import PathLike +from uuid import UUID + +import numpy as np +import numpy.typing as npt +from numpy import ma +from scipy.interpolate import CubicSpline, RectBivariateSpline + +from cloudnetpy import output, utils +from cloudnetpy.datasource import DataSource +from cloudnetpy.metadata import COMMON_ATTRIBUTES, MetaData +from cloudnetpy.utils import get_uuid + + +def generate_epsilon_from_radar( + radar_file: str | PathLike, + model_file: str | PathLike, + output_file: str | PathLike, + uuid: str | UUID | None = None, +) -> UUID: + """Generates Cloudnet radar-based dissipation rate of turbulent kinetic + energy product. + + Based on the pipeline of Griesche et al. (2020) with the inertial-subrange + slope-acceptance criterion of Borque et al. (2016). + + Args: + radar_file: Cloud radar L1b file name (provides Doppler velocity). + model_file: Cloudnet model file name (provides horizontal wind). + output_file: Output file name. + uuid: Set specific UUID for the file. + + Returns: + UUID of the generated file. + + Examples: + >>> from cloudnetpy.products import generate_epsilon_from_radar + >>> generate_epsilon_from_radar('radar.nc', 'ecmwf.nc', 'epsilon.nc') + + References: + Griesche, H. J., Seifert, P., Ansmann, A., Baars, H., Barrientos + Velasco, C., Bühl, J., Engelmann, R., Radenz, M., Zhenping, Y., and + Macke, A. (2020): Application of the shipborne remote sensing + supersite OCEANET for profiling of Arctic aerosols and clouds during + Polarstern cruise PS106, Atmos. Meas. Tech., 13, 5335-5358, + https://doi.org/10.5194/amt-13-5335-2020. + + Borque, P., Luke, E., and Kollias, P. (2016): On the unified + estimation of turbulence eddy dissipation rate using Doppler cloud + radars and lidars, J. Geophys. Res. Atmos., 120, 5972-5989, + https://doi.org/10.1002/2015JD024543. + """ + uuid = get_uuid(uuid) + with ( + EpsilonRadarSource(radar_file) as epsilon_source, + DataSource(model_file) as model_source, + ): + if epsilon_source.altitude is None: + msg = "Radar file is missing 'altitude' attribute." + raise ValueError(msg) + wind_interp = _get_wind_interpolator( + model_source, alt_site=float(epsilon_source.altitude) + ) + epsilon_source.append_epsilon(wind_interp) + epsilon_source.append_grid_variables() + date = epsilon_source.get_date() + attributes = output.add_time_attribute(EPSILON_RADAR_ATTRIBUTES, date) + output.update_attributes(epsilon_source.data, attributes) + output.save_product_file( + "epsilon-radar", + epsilon_source, + output_file, + uuid, + extra_sources=(model_source,), + ) + return uuid + + +def _get_wind_interpolator( + model: DataSource, alt_site: float +) -> Callable[[npt.NDArray, npt.NDArray], npt.NDArray]: + """Returns a bilinear interpolator for horizontal wind speed on + (time, height_amsl) where height is meters above MSL. + + Model heights are above the model's own surface; shift to absolute MSL so + they share the radar's height frame. In complex terrain the model grid + cell surface can differ from the site altitude, so prefer the model's own + surface field when available. + """ + uwind = model.getvar("uwind") + vwind = model.getvar("vwind") + surface_altitude = utils.get_model_surface_altitude(model.dataset, alt_site) + heights = model.to_m(model.dataset.variables["height"]) + surface_altitude + wind_speed = np.hypot(np.asarray(uwind), np.asarray(vwind)) + mean_height = np.asarray(ma.mean(heights, axis=0)) + common = np.empty((wind_speed.shape[0], mean_height.size)) + for i in range(wind_speed.shape[0]): + common[i] = np.interp(mean_height, np.asarray(heights[i]), wind_speed[i]) + return RectBivariateSpline(model.time, mean_height, common, kx=1, ky=1) + + +class EpsilonRadarSource(DataSource): + """Reads radar Doppler velocity and computes turbulent kinetic energy + dissipation rate (epsilon) on the product grid. + """ + + height: npt.NDArray + + def append_epsilon( + self, wind_interp: Callable[[npt.NDArray, npt.NDArray], npt.NDArray] + ) -> None: + """Estimate the dissipation rate (epsilon) on the 30 s product grid.""" + radar_time = np.asarray(self.getvar("time")) + self._radar_time = radar_time + height = self.height + radar_dt_sec = float(np.round(np.median(np.diff(radar_time)) * 3600.0)) + + v = self.getvar("v") + if isinstance(v, ma.MaskedArray): + v = v.filled(np.nan) + v = np.asarray(v, dtype=np.float64) + + product_time = utils.time_grid(time_step=PRODUCT_TIME_STEP_SEC) + self.time = product_time + wind_at_product = wind_interp(product_time, height) + + n_time = product_time.size + n_height = height.size + epsilon = np.full((n_time, n_height), EPSILON_INVALID, dtype=np.float64) + epsilon_error = np.full((n_time, n_height), EPSILON_INVALID, dtype=np.float64) + + starts = np.searchsorted(radar_time, product_time, side="right") + stops = np.searchsorted( + radar_time, product_time + AVERAGING_TIME_HR, side="right" + ) + + for time_idx in range(n_time): + after, stop = int(starts[time_idx]), int(stops[time_idx]) + if stop <= after: + continue + time_window = radar_time[after:stop] + min_valid = int(MIN_VALID_FRACTION * time_window.size) + + for height_idx in range(n_height): + vel = v[after:stop, height_idx] + nan_mask = np.isnan(vel) + if nan_mask[:2].all() or nan_mask[-2:].all(): + continue + if (vel.size - int(nan_mask.sum())) < min_valid: + continue + if nan_mask.any(): + valid = ~nan_mask + vel = CubicSpline(time_window[valid], vel[valid])(time_window) + + wind_speed = float(wind_at_product[time_idx, height_idx]) + if not np.isfinite(wind_speed) or wind_speed <= 0: + continue + + freq_sp, power_sp = _periodogram(vel, radar_dt_sec, wind_speed) + if freq_sp.size < 2: + continue + + result = _epsilon_from_spectrum(freq_sp, power_sp) + if result is not None: + eps, eps_err = result + epsilon[time_idx, height_idx] = eps + epsilon_error[time_idx, height_idx] = eps_err + + self.append_data( + ma.masked_where(epsilon == EPSILON_INVALID, epsilon), "epsilon" + ) + self.append_data( + ma.masked_where( + (epsilon_error == EPSILON_INVALID) | ~np.isfinite(epsilon_error), + epsilon_error, + ), + "epsilon_error", + ) + + def append_grid_variables(self) -> None: + """Adds time/height/altitude/latitude/longitude on the product grid. + + altitude/latitude/longitude are always written as 1-D arrays on the + product time grid: scalars from stationary radars are broadcast, + time-varying fields from moving platforms are rebinned. + """ + self.append_data(self.time, "time", dtype="f8") + self.append_data(np.asarray(self.height, dtype=np.float32), "height") + for key in ("altitude", "latitude", "longitude"): + if key not in self.dataset.variables: + continue + src = np.asarray(self.dataset.variables[key][:]) + if src.ndim == 0: + values = np.full(self.time.size, float(src), dtype=np.float32) + else: + values = utils.rebin_1d(self._radar_time, src, self.time) + self.append_data(values, key) + + +def _periodogram( + vel: npt.NDArray, delta_t_sec: float, adv_vel: float +) -> tuple[npt.NDArray, npt.NDArray]: + """Compute angular wavenumber spectrum from a Doppler velocity series.""" + n = vel.size + delta_x = delta_t_sec * adv_vel + fft_result = np.fft.rfft(vel)[1:] + power_sp = (delta_x / n) * np.abs(fft_result) ** 2 / np.pi + freq_sp = np.fft.rfftfreq(n, d=delta_x)[1:] * 2.0 * np.pi + keep = power_sp > 0 + return freq_sp[keep], power_sp[keep] + + +def _epsilon_from_spectrum( + freq_sp: npt.NDArray, + power_sp: npt.NDArray, +) -> tuple[float, float] | None: + """Vectorised slope fit on log-log spectrum across all frequency bands. + + Uses cumulative sums over a sorted spectrum so each band's least-squares + fit reduces to a constant-time index lookup. Returns ``(mean, std)`` of + epsilon across accepted bands -- the band-spread used by Griesche et al. + (2020) as the random retrieval uncertainty -- or ``None`` if no band + passes the slope filter. ``std`` is NaN when only a single band passes. + """ + log_f = np.log10(freq_sp) + log_p = np.log10(power_sp) + + csx = np.concatenate(([0.0], np.cumsum(log_f))) + csy = np.concatenate(([0.0], np.cumsum(log_p))) + csxx = np.concatenate(([0.0], np.cumsum(log_f * log_f))) + csxy = np.concatenate(([0.0], np.cumsum(log_f * log_p))) + + lo = np.searchsorted(freq_sp, _FMIN_ARR, side="right") + hi = np.searchsorted(freq_sp, _FMAX_ARR, side="left") + n = hi - lo + + sx = csx[hi] - csx[lo] + sy = csy[hi] - csy[lo] + sxx = csxx[hi] - csxx[lo] + sxy = csxy[hi] - csxy[lo] + + with np.errstate(invalid="ignore", divide="ignore"): + slope = (n * sxy - sx * sy) / (n * sxx - sx * sx) + intercept = (sy - slope * sx) / n + + valid = ( + (n >= 2) + & np.isfinite(slope) + & (slope > THRESHOLD_SLOPE_MIN) + & (slope < THRESHOLD_SLOPE_MAX) + ) + if not valid.any(): + return None + epsilon = (10.0 ** intercept[valid] / KOLMOGOROV_CONSTANT) ** 1.5 + std = float(np.std(epsilon, ddof=1)) if epsilon.size >= 2 else float("nan") + return float(np.mean(epsilon)), std + + +KOLMOGOROV_CONSTANT = 0.5 +THRESHOLD_SLOPE_MIN = -2.0 # -5/3 - 20% (Borque, 2016) +THRESHOLD_SLOPE_MAX = -1.33 # -5/3 + 20% (Borque, 2016) +AVERAGING_TIME_HR = 5.0 / 60.0 +PRODUCT_TIME_STEP_SEC = 30 +MIN_VALID_FRACTION = 0.9 +EPSILON_INVALID = -999.0 + + +COMMENTS = { + "general": ( + "This dataset contains the dissipation rate of turbulent kinetic\n" + "energy calculated using the pipeline of Griesche et al. (2020)\n" + "with the inertial-subrange slope-acceptance criterion of Borque\n" + "et al. (2016). The turbulent energy spectrum is derived as the\n" + "power spectrum of cloud radar Doppler velocity over 5-minute\n" + "windows and converted to a wavenumber spectrum using model\n" + "horizontal wind via Taylor's hypothesis. The dissipation rate is\n" + "recovered from the intercept of a log-log fit across multiple\n" + "frequency bands; bands whose slope falls outside -5/3 +/- 20% are\n" + "rejected and the accepted bands are averaged.\n" + ), + "epsilon": ( + "This variable was calculated for profiles where 5 minutes of\n" + "continuous cloud radar Doppler velocity was available. With less\n" + "than 10% missing data, the missing values were filled by cubic\n" + "spline interpolation.\n" + ), + "epsilon_error": ( + "Random uncertainty estimated as the standard deviation of epsilon\n" + "across the frequency bands accepted by the inertial-subrange\n" + "slope criterion (Griesche et al. 2020). It captures the scatter\n" + "introduced by the multi-band fit but does not include systematic\n" + "contributions from horizontal wind uncertainty or the Kolmogorov\n" + "constant. The variable is masked where only a single band passed\n" + "the slope filter and no spread can be estimated.\n" + ), +} + + +EPSILON_RADAR_ATTRIBUTES = { + "comment": COMMENTS["general"], + "epsilon": MetaData( + long_name="Dissipation rate of turbulent kinetic energy", + units="m2 s-3", + ancillary_variables="epsilon_error", + comment=COMMENTS["epsilon"], + dimensions=("time", "height"), + ), + "epsilon_error": MetaData( + long_name="Absolute error in dissipation rate of turbulent kinetic energy", + units="m2 s-3", + comment=COMMENTS["epsilon_error"], + dimensions=("time", "height"), + ), + "height": COMMON_ATTRIBUTES["height"]._replace(dimensions=("height",)), +} + + +# Frequency-range pairs (rad/m) used to fit the inertial subrange slope. +# Adapted from Griesche et al. 2020 with a multi-decade extension. +freq_array = np.array( + [ + [5e-3, 1.5e-2], + [5e-3, 3e-2], + [5e-3, 5e-2], + [8e-3, 2e-2], + [8e-3, 4e-2], + [8e-3, 6e-2], + [2e-2, 7e-2], + [2e-2, 1e-1], + [2e-2, 1.5e-1], + [3e-2, 1e-1], + [3e-2, 2e-1], + [3e-2, 3e-1], + [5e-2, 1.5e-1], + [5e-2, 3e-1], + [5e-2, 5e-1], + [8e-2, 2e-1], + [8e-2, 4e-1], + [8e-2, 6e-1], + [2e-1, 7e-1], + [2e-1, 1e0], + [2e-1, 1.5e0], + [3e-1, 1e0], + [3e-1, 2e0], + [3e-1, 3e0], + [5e-1, 1.5e0], + [5e-1, 3e0], + [5e-1, 5e0], + [8e-1, 2e0], + [8e-1, 4e0], + [8e-1, 6e0], + [2e0, 7e0], + [2e0, 1e1], + [2e0, 1.5e1], + [3e0, 1e1], + [3e0, 2e1], + [3e0, 3e1], + [5e0, 1.5e1], + [5e0, 3e1], + [5e0, 5e1], + [8e0, 2e1], + [8e0, 4e1], + [8e0, 6e1], + ] +) +_FMIN_ARR = freq_array[:, 0] +_FMAX_ARR = freq_array[:, 1] diff --git a/cloudnetpy/utils.py b/cloudnetpy/utils.py index 1bd9ecde..8bf54f18 100644 --- a/cloudnetpy/utils.py +++ b/cloudnetpy/utils.py @@ -25,7 +25,7 @@ ) from cloudnetpy.cloudnetarray import CloudnetArray -from cloudnetpy.constants import SEC_IN_DAY, SEC_IN_HOUR, SEC_IN_MINUTE +from cloudnetpy.constants import SEC_IN_DAY, SEC_IN_HOUR, SEC_IN_MINUTE, G from cloudnetpy.exceptions import ValidTimeStampError @@ -875,6 +875,27 @@ def is_empty_line(line: str) -> bool: return line in ("\n", "\r\n") +def get_model_surface_altitude( + dataset: netCDF4.Dataset, alt_site: float +) -> float | npt.NDArray: + """Returns model surface altitude (m), preferring the model's own field. + + For sites in complex terrain (e.g. mountains), the model grid cell + surface height can differ significantly from the actual site altitude. + Using the model's own surface height ensures that thermodynamic and + wind fields are placed at their physically correct absolute heights. + + Returns a (n_time, 1) array when the model provides a per-timestep + surface field, otherwise the scalar ``alt_site`` fallback. Both shapes + broadcast against (n_time, n_height) model arrays. + """ + if "sfc_height" in dataset.variables: + return dataset.variables["sfc_height"][:][:, np.newaxis] + if "sfc_geopotential" in dataset.variables: + return dataset.variables["sfc_geopotential"][:][:, np.newaxis] / G + return alt_site + + def is_timestamp(timestamp: str) -> bool: """Tests if the input string is formatted as -yyyy-mm-dd hh:mm:ss.""" reg_exp = re.compile(r"-\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}")