"""Module for reading raw cloud radar data."""
import datetime
import logging
import os
import re
from collections import OrderedDict
from collections.abc import Sequence
from os import PathLike
from tempfile import NamedTemporaryFile, TemporaryDirectory
from uuid import UUID
import numpy as np
from numpy import ma
from cloudnetpy import concat_lib, output, utils
from cloudnetpy.cloudnetarray import CloudnetArray
from cloudnetpy.constants import HZ_TO_GHZ
from cloudnetpy.instruments.instruments import MIRA10, MIRA35, MIRA35C, MIRA35S
from cloudnetpy.instruments.nc_radar import NcRadar
from cloudnetpy.metadata import MetaData
[docs]
def mira2nc(
raw_mira: str | PathLike | Sequence[str | PathLike],
output_file: str | PathLike,
site_meta: dict,
uuid: str | UUID | None = None,
date: str | datetime.date | None = None,
) -> UUID:
"""Converts METEK MIRA cloud radar data into Cloudnet Level 1b netCDF file.
This function converts raw MIRA file(s) into a much smaller file that
contains only the relevant data and can be used in further processing
steps.
Args:
raw_mira: Filename of a daily MIRA ``.mmclx`` or ``.znc`` file. Can also
be a folder containing several non-concatenated ``.mmclx`` or
``.znc`` files from one day, or a list of files. ``.znc`` files take
precedence because they are the newer filetype.
output_file: Output filename.
site_meta: Dictionary containing information about the site. Required
key-value pair is `name`. Optional `model` takes `mira-35` (default),
`mira-35s`, `mira-35c`, or `mira-10`.
uuid: Set specific UUID for the file.
date: Expected date as YYYY-MM-DD of all profiles in the file.
Returns:
UUID of the generated file.
Raises:
ValidTimeStampError: No valid timestamps found.
FileNotFoundError: No suitable input files found.
ValueError: Wrong suffix in input file(s).
TypeError: Mixed mmclx and znc files.
Examples:
>>> from cloudnetpy.instruments import mira2nc
>>> site_meta = {'name': 'Vehmasmaki'}
>>> mira2nc('raw_radar.mmclx', 'radar.nc', site_meta)
>>> mira2nc('raw_radar.znc', 'radar.nc', site_meta)
>>> mira2nc('/one/day/of/mira/mmclx/files/', 'radar.nc', site_meta)
>>> mira2nc('/one/day/of/mira/znc/files/', 'radar.nc', site_meta)
"""
if isinstance(date, str):
date = datetime.date.fromisoformat(date)
uuid = utils.get_uuid(uuid)
with TemporaryDirectory() as temp_dir:
input_filename, keymap = _parse_input_files(raw_mira, temp_dir)
with Mira(input_filename, site_meta) as mira:
mira.init_data(keymap)
if date is not None:
mira.screen_by_date(date)
mira.date = date
mira.sort_timestamps()
mira.remove_duplicate_timestamps()
mira.linear_to_db(("Zh", "ldr", "SNR"))
mira.screen_low_power()
if "snr_limit" in site_meta and site_meta["snr_limit"] is not None:
snr_limit = site_meta["snr_limit"]
else:
# Empirical values, should be checked
snr_limit = -30 if mira.instrument == MIRA10 else -17
# Old MIRA files don't have angle variables.
if "elevation" not in mira.data:
mira.append_data(ma.masked_all_like(mira.time.data), "elevation")
if "azimuth_angle" not in mira.data:
mira.append_data(ma.masked_all_like(mira.time.data), "azimuth_angle")
mira.screen_by_snr(snr_limit)
mira.screen_invalid_ldr()
mira.mask_invalid_data()
mira.mask_bad_angles()
mira.add_time_and_range()
mira.add_site_geolocation()
mira.add_radar_frequency()
mira.broadcast_nyquist_velocity()
valid_indices = mira.add_zenith_and_azimuth_angles(
elevation_threshold=1.1,
elevation_diff_threshold=1e-6,
azimuth_diff_threshold=1e-3,
zenith_offset=site_meta.get("zenith_offset"),
azimuth_offset=site_meta.get("azimuth_offset"),
)
mira.screen_time_indices(valid_indices)
mira.add_height()
mira.test_if_all_masked()
attributes = output.add_time_attribute(ATTRIBUTES, mira.date)
output.update_attributes(mira.data, attributes)
output.save_level1b(mira, output_file, uuid)
return uuid
class Mira(NcRadar):
"""Class for MIRA raw radar data. Child of NcRadar().
Args:
full_path: Filename of a daily MIRA .mmclx NetCDF file.
site_meta: Site properties in a dictionary. Required keys are: `name`.
"""
epoch = datetime.datetime(1970, 1, 1, tzinfo=datetime.timezone.utc)
def __init__(self, full_path: str | PathLike, site_meta: dict) -> None:
super().__init__(full_path, site_meta)
self.date = self._init_mira_date()
self.hrd = self._read_hrd()
if "model" not in site_meta or site_meta["model"] == "mira-35":
self.instrument = MIRA35
elif site_meta["model"] == "mira-35s":
self.instrument = MIRA35S
elif site_meta["model"] == "mira-35c":
self.instrument = MIRA35C
elif site_meta["model"] == "mira-10":
self.instrument = MIRA10
else:
msg = f"Invalid model: {site_meta['model']}"
raise ValueError(msg)
def _read_hrd(self) -> dict[str, str]:
lines = self.dataset.hrd.split("\n")
output = {}
for line in lines:
if ":" not in line:
continue
key, value = line.split(":", maxsplit=1)
output[key] = value.strip()
return output
def add_radar_frequency(self) -> None:
# XMT should always be present. In newer files, there's also
# LO_frequency variable which doesn't match XMT exactly and sometimes
# contains only zeros.
key = "radar_frequency"
frequency = float(self.hrd["XMT"]) * HZ_TO_GHZ
self.data[key] = CloudnetArray(frequency, key)
def broadcast_nyquist_velocity(self) -> None:
# Scalar in original file, 1d array in concatenated file.
key = "nyquist_velocity"
nv = self.data[key].data
if len(nv.shape) == 1:
self.data[key].data = ma.getdata(nv)[:, np.newaxis]
def screen_by_date(self, expected_date: datetime.date) -> None:
"""Screens incorrect time stamps."""
time_stamps = self.getvar("time")
valid_indices = []
for ind, timestamp in enumerate(time_stamps):
if not timestamp:
continue
date = utils.seconds2date(timestamp, self.epoch).date()
if date == expected_date:
valid_indices.append(ind)
self.screen_time_indices(valid_indices)
def _init_mira_date(self) -> datetime.date:
time_stamps = self.getvar("time")
return utils.seconds2date(float(time_stamps[0]), self.epoch).date()
def screen_low_power(self) -> None:
"""Screen times with low average transmit power (tpow).
The expected average power depends on the specific model:
- MIRA-35 / 35S: average power should be between 30 and 60 W according
to the data sheet. Based on a random sample, tpow values range from 15
to 25 W.
- MIRA-35C: average power should be around 3 W according to the data
sheet. Based on a random sample, tpow values range from 1.5 to 2 W.
- MIRA-10: average power is up to 40 W (website) or 50 W (data sheet).
Please note that the raw tpow values cannot always be trusted:
- Files with SN:fzk and without FZK100 field in 'hrd' global attribute
should be multiplied by 100.
- It looks like tpow could contain the peak power instead of the average
power. As of now, this has only be seen in Munich MIRA-10, but because
this is the only MIRA-10 in Cloudnet, let's leave it as it is.
Example case:
- Limassol 2024-10-20: https://hdl.handle.net/21.12132/1.159fe518fe5b403d
"""
if "tpow" not in self.data:
logging.warning("Variable tpow is missing")
return
self._correct_fzk_tpow()
tpow = self.data["tpow"][:]
low_threshold = 0.5 if self.instrument == MIRA35C else 5
is_low = tpow < low_threshold
n_removed = np.count_nonzero(is_low)
if n_removed > 0:
logging.warning(
"Filtering %s profiles due to low average transmit power", n_removed
)
self.screen_time_indices(~is_low)
def _correct_fzk_tpow(self) -> None:
"""Corrects tpow for old instruments missing the FZK100 fix.
In old MIRA-35 files with 'SN:fzk', tpow was reported 100x too small
(e.g. Schneefernerhaus 2012-01-01). The tpow was corrected in a firmware
update as indicated by 'FZK100: 1.00000' in the 'hrd' global
attribute (e.g. Schneefernerhaus 2018-01-01).
Please note that this doesn't affect all old MIRA-35 instruments (e.g.
'SN:dwd', Lindenberg 2007-01-01).
"""
if self.hrd["SN"] != "fzk" or "FZK100" in self.hrd:
return
logging.info("Correcting tpow by factor 100 for old FZK instrument firmware")
self.data["tpow"].data[:] *= 100
self.data["tpow"].correction_factor = 100
def screen_invalid_ldr(self) -> None:
"""Masks LDR in MIRA STSR mode data.
Is there a better way to identify this mode?
"""
if "ldr" not in self.data:
return
# Delete LDR if polarization is off. The LDR variable exists but
# contains no data. At least Munich MIRA-10 doesn't have polarization.
if self.hrd["POL"] == "0":
del self.data["ldr"]
return
ldr = self.data["ldr"][:]
if ma.mean(ldr) > 0:
logging.warning(
"LDR values suspiciously high. Mira in STSR mode? "
"Screening all LDR for now.",
)
self.data["ldr"].data[:] = ma.masked
def mask_bad_angles(self) -> None:
"""Masks clearly bad elevation and azimuth angles."""
limits = {
"elevation": (0, 180),
"azimuth_angle": (-360, 360),
}
for key, (lower, upper) in limits.items():
if (array := self.data[key].data) is not None:
margin = (upper - lower) * 0.05
array[array < (lower - margin)] = ma.masked
array[array > (upper + margin)] = ma.masked
def _parse_input_files(
input_files: str | PathLike | Sequence[str | PathLike], temp_dir: str
) -> tuple[str | PathLike, dict[str, str]]:
input_filename: str | PathLike
if (
not isinstance(input_files, str) and isinstance(input_files, Sequence)
) or os.path.isdir(input_files):
with NamedTemporaryFile(
dir=temp_dir,
suffix=".nc",
delete=False,
) as temp_file:
input_filename = temp_file.name
if not isinstance(input_files, str) and isinstance(input_files, Sequence):
valid_files = sorted(map(str, input_files))
else:
valid_files = utils.get_sorted_filenames(input_files, ".znc")
if not valid_files:
valid_files = utils.get_sorted_filenames(input_files, ".mmclx")
if not valid_files:
msg = (
(
f"Neither znc nor mmclx files found {input_files}. "
f"Please check your input."
),
)
raise FileNotFoundError(msg)
filetypes = list({_get_suffix(f) for f in valid_files})
if len(filetypes) > 1:
err_msg = "Mixed mmclx and znc files. Please use only one filetype."
raise TypeError(err_msg)
keymap = _get_keymap(filetypes[0])
variables = list(keymap.keys())
concat_lib.concatenate_files(
valid_files,
input_filename,
variables=variables,
ignore=_get_ignored_variables(filetypes[0]),
)
else:
input_filename = input_files
keymap = _get_keymap(_get_suffix(input_filename))
return input_filename, keymap
def _get_ignored_variables(filetype: str) -> list | None:
"""Returns variables to ignore for METEK MIRA-35 cloud radar concat."""
_check_file_type(filetype)
# Ignore spectral variables for now
keymaps = {
"znc": ["DropSize", "SPCco", "SPCcx", "SPCcocxRe", "SPCcocxIm", "doppler"],
"mmclx": None,
}
return keymaps.get(filetype.lower(), keymaps.get("mmclx"))
def _get_suffix(filename: str | PathLike) -> str:
m = re.search(r"\.(\w+)(\.\d+)?$", str(filename))
if m is None:
return ""
return m[1].lower()
def _get_keymap(filetype: str) -> dict[str, str]:
"""Returns a dictionary mapping the variables in the raw data to the processed
Cloudnet file.
"""
_check_file_type(filetype)
# Order is relevant with the new znc files from STSR radar
keymaps = {
"znc": OrderedDict(
[
("Zg", "Zh"), # fallback
("Zh2l", "Zh"),
("VELg", "v"), # fallback
("VELh2l", "v"),
("RMSg", "width"), # fallback
("RMSh2l", "width"),
("LDRg", "ldr"), # fallback
("LDRh2l", "ldr"),
("SNRg", "SNR"), # fallback
("SNRh2l", "SNR"),
("elv", "elevation"),
("azi", "azimuth_angle"),
("nfft", "nfft"),
("nave", "nave"),
("prf", "prf"),
("rg0", "rg0"),
("NyquistVelocity", "nyquist_velocity"),
("tpow", "tpow"),
],
),
"mmclx": OrderedDict(
[
("Ze", "Zh"), # fallback for old mmclx files
("Zg", "Zh"),
("VELg", "v"),
("RMSg", "width"),
("LDRg", "ldr"),
("SNRg", "SNR"),
("elv", "elevation"),
("azi", "azimuth_angle"),
("nfft", "nfft"),
("nave", "nave"),
("prf", "prf"),
("rg0", "rg0"),
("NyquistVelocity", "nyquist_velocity"),
("tpow", "tpow"),
]
),
}
return keymaps.get(filetype.lower(), keymaps["mmclx"])
def _check_file_type(filetype: str) -> None:
known_filetypes = ["znc", "mmclx"]
if filetype.lower() not in known_filetypes:
msg = f"Filetype must be one of {known_filetypes}"
raise ValueError(msg)
ATTRIBUTES = {
"nfft": MetaData(long_name="Number of FFT points", units="1", dimensions=("time",)),
"nave": MetaData(
long_name="Number of spectral averages (not accounting for overlapping FFTs)",
units="1",
dimensions=("time",),
),
"rg0": MetaData(
long_name="Number of lowest range gates", units="1", dimensions=("time",)
),
"prf": MetaData(
long_name="Pulse Repetition Frequency", units="Hz", dimensions=("time",)
),
"tpow": MetaData(
long_name="Average Transmit Power", units="W", dimensions=("time",)
),
"zenith_offset": MetaData(
long_name="Zenith offset of the instrument",
units="degrees",
comment="Zenith offset applied.",
dimensions=None,
),
"azimuth_offset": MetaData(
long_name="Azimuth offset of the instrument (positive clockwise from north)",
units="degrees",
comment="Azimuth offset applied.",
dimensions=None,
),
}