"""
Tools
=====
.. autosummary::
:nosignatures:
:toctree: generated/
{}
"""
__all__ = [
"add_axis",
"concat_sweep_angle",
"create_filelist",
"create_lineplot",
"data_path",
"estimate_cmap",
"filter_data",
"filter_radar_bins",
"filter_stations",
"gauss_kernel",
"get_beamblockage",
"get_bins_from_range",
"get_discrete_cmap",
"get_discrete_norm",
"get_polar_ax",
"get_qbbf",
"get_range_step",
"get_temperature_profile",
"get_xpol_path",
"interpolate_dem",
"load_single_sweep",
"load_timeseries_sweep",
"load_volume",
"merge_radar_profile",
"noise_correction",
"noise_hist",
"noise_level",
"phase_offset",
"phase_zphi",
"plot_2by2",
"plot_beamprop",
"plot_matched_refl",
"plot_noise_hist",
"plot_unmatched_refl",
"plot_moment",
"plot_reference_colorbar",
"process_regression",
"ProgressBar",
"read_dwd_stations",
"rpol_data",
"save_dwd_station",
"set_sweep_mode_as_coords",
"smooth_data",
"timeseries_sweep_to_zarr",
"visdict14",
"volume_matching",
"xr_rolling",
"x_cumtrapz",
"zphi",
]
__doc__ = __doc__.format("\n ".join(__all__))
import datetime as dt
import fnmatch
import glob
import os
import pathlib
import re
import shutil
from contextlib import nullcontext
from functools import partial
from html.parser import HTMLParser
import certifi
import dask
import fsspec
import geopandas as gpd
import matplotlib as mpl
import matplotlib.pyplot as plt
import numcodecs
import numpy as np
import pandas as pd
import pyproj
import scipy
import shapely
import urllib3
import wradlib as wrl
import xarray as xr
import xradar as xd
import zarr
from dask.diagnostics import ProgressBar as _ProgressBar
from matplotlib.colors import LogNorm
from matplotlib.path import Path
from osgeo import osr
from scipy.spatial import KDTree
from tqdm.autonotebook import tqdm
from xhistogram.xarray import histogram
os.environ["WRADLIB_DATA"] = "/automount/ftp/radar/wradlib-data"
os.environ["CPL_DEBUG"] = "OFF"
ProgressBar = nullcontext if os.environ.get("SPHINX_BUILD") else _ProgressBar
###
## Data IO
def repo_path(start=None):
path = pathlib.Path(start or pathlib.Path.cwd()).resolve()
for parent in [path, *path.parents]:
if (parent / "pyproject.toml").exists() or (parent / ".git").exists():
return parent
# if we do not find anything, we just take the current path
return path
[docs]
def data_path(start=None):
env = os.environ.get("RP_DATA")
if env:
return pathlib.Path(env).resolve()
path = repo_path(start)
return path / "data"
os.environ["RP_DATA"] = str(data_path())
rp_mode = os.environ.get("RP_MODE", "ftp")
if rp_mode == "ftp":
target_protocol = "ftp"
target_options = {
"host": "ftp.meteo.uni-bonn.de",
"username": "anonymous",
"password": "anonymous@",
}
cache_storage = str(data_path() / ".cache/ftp.meteo.uni-bonn.de")
elif rp_mode == "automount":
target_protocol = "file"
target_options = {}
cache_storage = str(data_path() / ".cache/automount")
else:
raise ValueError("Unknown RP mode!")
rpol_data = fsspec.filesystem(
"filecache",
target_protocol=target_protocol,
target_options=target_options,
cache_storage=cache_storage,
)
class DWDHTMLParser(HTMLParser):
def handle_starttag(self, tag, attrs):
if tag != "a":
return
self.links.append(attrs[0][1])
[docs]
def get_xpol_path(path=None, start_time=None, loc="boxpol", fs=None):
"""Create path of BoXPol/JuXPol radar data files.
Parameter
---------
path : str
Path to root folder of radar data,
defaults to None (/automount/radar/scans or /automount/radar-archiv/scans)
start_time : datetime.datetime
datetime - object to select correct folder
loc : str
"boxpol" or "juxpol" (not case sensitive)
Return
------
radar_path : str
Path to radar data
"""
if start_time is None:
start_time = dt.datetime.today()
loc = "" if loc.lower()[0:2] == "bo" else "_juelich"
if fs is None:
fs = rpol_data
if path is None:
if fs.protocol == "filecache":
_fs = fs.fs
else:
_fs = fs
if _fs.protocol == ("file", "local"):
ins = "-archiv" if start_time < dt.datetime(2015, 1, 1) else ""
path = f"/automount/radar{ins}/scans{loc}"
if not fs.exists(path):
path = os.environ["RADARMET_DATA"]
else:
path = "/radar/rp-data/scans"
radar_path = os.path.join(path, "{0}/{0}-{1:02}/{0}-{1:02d}-{2:02d}")
return radar_path.format(start_time.year, start_time.month, start_time.day)
def get_file_date_regex(filename):
"""Get regex from filename."""
# regex for ""%Y-%m-%d--%H:%M:%S"
reg0 = r"\d{4}.\d{2}.\d{2}..\d{2}.\d{2}.\d{2}"
# regex for "%Y%m%d%H%M%S"
reg1 = r"\d{14}"
match = re.search(reg0, os.path.basename(filename))
return reg1 if match is None else reg0
def get_datetime_from_filename(filename, regex):
"""Get datetime from filename."""
fmt = "%Y%m%d%H%M%S"
match = re.search(regex, os.path.basename(filename))
match = "".join(re.findall(r"[0-9]+", match.group()))
return dt.datetime.strptime(match, fmt)
[docs]
def create_filelist(path_glob, starttime, endtime, fs=None):
"""Create filelist from path_glob and filename dates in a location-agnostic manner."""
if fs is None:
fs = rpol_data
file_names = sorted(fs.glob(path_glob))
regex = get_file_date_regex(file_names[0])
for fname in file_names:
time = get_datetime_from_filename(fname, regex)
if starttime <= time < endtime:
yield fname
[docs]
def load_single_sweep(start_time, scan, reindex_angle=None, fs=None):
stop_time = start_time + dt.timedelta(minutes=5)
path = get_xpol_path(start_time=start_time, fs=fs)
file_path = os.path.join(path, scan)
flist = list(
create_filelist(os.path.join(file_path, "*"), start_time, stop_time, fs=fs)
)
if reindex_angle is None:
reindex_angle = dict(start_angle=0, stop_angle=360, angle_res=1.0, direction=1)
if fs is None:
fs = rpol_data
with fs.open(flist[0], "rb") as f:
swp = xr.open_dataset(
f.name, engine="gamic", group="sweep_0", reindex_angle=reindex_angle
)
swp = swp.set_coords("sweep_mode")
swp = swp.wrl.georef.georeference()
return swp
[docs]
def timeseries_sweep_to_zarr(
start_time,
delta,
scan,
zarr_name,
compressor=None,
filters=None,
dask_scheduler=None,
reindex_angle=None,
fs=None,
):
if dask_scheduler is None:
dask.config.set(scheduler="threads")
else:
dask.config.set(scheduler=dask_scheduler)
if compressor is None:
compressor = numcodecs.Zstd(level=2)
else:
compressor = compressor
stop_time = start_time + delta
if reindex_angle is None:
reindex_angle = dict(start_angle=0, stop_angle=360, angle_res=1.0, direction=1)
if fs is None:
fs = rpol_data
path = get_xpol_path(start_time=start_time, fs=fs)
file_path = os.path.join(path, scan)
flist = list(
create_filelist(os.path.join(file_path, "*"), start_time, stop_time, fs=fs)
)
valid_keys = {"dtype", "scale_factor", "add_offset", "_FillValue"}
for i, f in enumerate(tqdm(flist, desc="Writing consolidated zarr")):
with fs.open(f, "rb") as fh:
local_path = fh.name
ds = xr.open_dataset(
local_path,
engine="gamic",
group="sweep_0",
reindex_angle=reindex_angle,
chunks=-1,
)
ds = ds.set_coords("sweep_mode")
ds = ds.reset_coords(["time", "elevation"])
ds = ds.assign_coords({"time2": ds.time.load().min()})
ds = ds.expand_dims("time2")
ds.time2.encoding = ds.time.encoding.copy()
if i == 0:
encoding = {}
for v in ds.data_vars:
enc = {k: v2 for k, v2 in ds[v].encoding.items() if k in valid_keys}
var_chunks = ds[v].chunks
if var_chunks is None:
zarr_chunks = ds[v].shape
else:
zarr_chunks = tuple(
c[0] if isinstance(c, tuple) else c for c in var_chunks
)
enc["chunks"] = zarr_chunks
enc["compressor"] = compressor
enc["filters"] = filters
encoding[v] = enc
ds.to_zarr(
zarr_name,
mode="w",
consolidated=False,
encoding=encoding,
zarr_version=2,
)
else:
ds.to_zarr(
zarr_name,
mode="a",
append_dim="time2",
consolidated=False,
zarr_version=2,
)
zarr.consolidate_metadata(zarr_name)
# shutil.rmtree(cache_dir, ignore_errors=True)
def preprocess(ds, keep_vars=None):
if keep_vars is not None:
to_drop = [v for v in ds.data_vars if v not in keep_vars]
ds = ds.drop_vars(to_drop, errors="ignore")
# ds = ds.set_coords("sweep_mode")
# ds = ds.assign_coords(time2=ds.time.min())
# dim0 = ds.wrl.util.dim0()
# dim1 = {"azimuth": "elevation", "elevation": "azimuth"}[dim0]
# angle_median = ds[dim1].load().median(dim0, skipna=True)
# angle = xr.where(ds[dim1] != angle_median, angle_median, ds[dim1])
# ds = ds.assign_coords({dim1: angle})
return ds
[docs]
def load_timeseries_sweep(
start_time, delta, scan, reindex_angle=None, keep_vars=None, fs=None
):
if fs is None:
fs = rpol_data
if reindex_angle is None:
reindex_angle = dict(start_angle=0, stop_angle=360, angle_res=1.0, direction=1)
stop_time = start_time + delta
path = get_xpol_path(start_time=start_time, fs=fs)
file_path = os.path.join(path, scan)
flist = list(
create_filelist(os.path.join(file_path, "*"), start_time, stop_time, fs=fs)
)
print(f"loading {len(flist)} files")
preprocess_func = partial(preprocess, keep_vars=keep_vars)
_flist = []
for _i, f in enumerate(tqdm(flist, desc="Claiming files")):
with fs.open(f, "rb") as fh:
_flist.append(fh.name)
sweep_ts = xr.open_mfdataset(
_flist,
concat_dim="time2",
combine="nested",
engine="gamic",
group="sweep_0",
reindex_angle=reindex_angle,
preprocess=preprocess_func,
parallel=True,
)
# sweep_ts = sweep_ts.wrl.georef.georeference()
return sweep_ts
[docs]
def load_volume(flist, reindex_angle=None, fs=None):
if fs is None:
fs = rpol_data
fnames = []
for fpath in flist:
with fs.open(fpath) as f:
fnames.append(f.name)
sweeps = [
xr.open_dataset(
fname, group="sweep_0", engine="gamic", reindex_angle=reindex_angle
)
for fname in fnames
]
sweeps.insert(0, xr.open_dataset(fnames[0], group="/"))
dtree = xr.DataTree(dataset=xd.io.backends.common._assign_root(sweeps), name="root")
dtree = xd.io.backends.common._attach_sweep_groups(dtree, sweeps[1:])
for key in list(dtree.children):
if "sweep" in key:
ds = dtree[key].to_dataset().set_coords("sweep_mode")
ds = ds.assign_coords(elevation=ds.elevation.median())
ds = ds.wrl.georef.georeference()
dtree[key].ds = ds
return dtree
[docs]
def concat_sweep_angle(dtree, rres):
output = []
for key in list(dtree.children):
if "sweep" in key:
ds = dtree[key].to_dataset()
ds = ds.reindex(
range=np.arange(rres / 2, ds.range.max() + rres / 2, rres),
method="nearest",
)
output.append(ds)
return xr.concat(output, "sweep_fixed_angle")
## Plotting
def set_title(da, ax=None):
if ax is None:
ax = plt.gca()
da_time = da.time.min().dt.strftime("%Y-%m-%dT%H:%M:%S")
ax.set_title(f"{da.name} - {da_time.values}")
[docs]
def get_polar_ax(fig=None, ax=None):
if fig is None:
fig = plt.figure(figsize=(10, 10))
if ax is None:
ax = 111
ax = fig.add_subplot(ax, projection="polar")
# set the lable go clockwise and start from the top
ax.set_theta_zero_location("N")
# clockwise
ax.set_theta_direction(-1)
# set ticks proper
ax.set_xticks(np.pi / 180.0 * np.linspace(0, 360, 36, endpoint=False))
return ax
def plot_polar(data, ax, **kwargs):
theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
datap = data.assign_coords(
theta=(["azimuth"], theta, {"standard_name": "azimuth angle"})
)
pm = datap.plot(x="theta", ax=ax, **kwargs)
ax.set_ylabel("")
return pm
[docs]
def get_discrete_cmap(ticks, colors, bad="white", over=None, under=None):
"""Create discrete colormap.
Parameters
----------
ticks : int | sequence
number of ticks or sequence of ticks
colors : colormap | sequence
colormap or sequence of colors
bad : color
over : color
under : color
Returns
-------
matplotlib.colors.ListedColormap
"""
ticks = ticks if isinstance(ticks, int) else len(ticks)
if isinstance(colors, (str, mpl.colors.Colormap)):
cmap = mpl.colormaps.get_cmap(colors)
colors = cmap(np.linspace(0, 1, ticks + 1))
cmap = mpl.colors.ListedColormap(colors[1:-1])
if over is None:
over = colors[-1]
if under is None:
under = colors[0]
cmap.set_under(under)
cmap.set_over(over)
cmap.set_bad(color=bad)
return cmap
[docs]
def get_discrete_norm(ticks):
"""Return discrete boundary norm.
Parameters
----------
ticks : sequence
sequence of ticks
Returns
-------
matplotlib.colors.BoundaryNorm
"""
return mpl.colors.BoundaryNorm(ticks, len(ticks) - 1)
[docs]
def plot_reference_colorbar(ticks, cmap, ax, **kwargs):
"""Plot reference colorbar
Parameters
----------
ticks : sequence
sequence of ticks
cmap : cmap instance
ax : axes instance
Keyword Arguments
-----------------
Additonal Keywords for colorbar
Returns
-------
matplotlib.colorbar instance
"""
vmin = ticks[0] - np.diff(ticks)[0]
vmax = ticks[-1] + np.diff(ticks)[-1]
vmax2 = ticks[-1] + 2 * np.diff(ticks)[-1]
y = [vmin] + list(ticks) + [vmax] + [vmax2]
x = np.arange(2)
data = np.repeat(np.array(y, dtype=float), 2).reshape(-1, 2)
data[-2, :] = np.nan
norm = get_discrete_norm(ticks)
pm = ax.pcolormesh(
x,
y,
data,
cmap=cmap,
norm=norm,
)
cb = plt.colorbar(
pm,
ax=ax,
ticks=ticks,
**kwargs,
)
return cb
colors_prabhakar = np.array(
[
[0.00, 1.00, 1.00],
[0.00, 0.70, 0.93],
[0.00, 0.00, 1.00],
[0.50, 1.00, 0.00],
[0.40, 0.80, 0.00],
[0.27, 0.55, 0.00],
[1.00, 1.00, 0.00],
[0.80, 0.80, 0.00],
[1.00, 0.65, 0.00],
[1.00, 0.27, 0.00],
[0.80, 0.22, 0.00],
[0.55, 0.15, 0.00],
[1.00, 0.00, 1.00],
[0.58, 0.44, 0.86],
]
)
cmap_prabhakar = mpl.colors.ListedColormap(colors_prabhakar)
mpl.colormaps.register(cmap_prabhakar, name="miub")
#: Configuration dictionary for visualizing radar data.
visdict14 = dict(
ZH=dict(
ticks=np.arange(-10, 55, 5),
contours=[0, 5, 10, 15, 20, 25, 30, 35],
cmap=cmap_prabhakar,
name=r"Horizontal Reflectivity (dBz)",
short=r"$\mathrm{\mathsf{Z_{H}}}$",
),
ZDR=dict(
ticks=np.array(
[-1.0, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 2.0, 3.0]
),
contours=np.array([-1, -0.3, 0]),
cmap=cmap_prabhakar,
name=r"Differential Reflectivity (dB)",
short=r"$\mathrm{\mathsf{Z_{DR}}}$",
),
RHO=dict(
ticks=np.array(
[
0.7,
0.8,
0.85,
0.9,
0.92,
0.94,
0.95,
0.96,
0.97,
0.98,
0.99,
0.995,
0.998,
]
),
cmap=cmap_prabhakar,
name=r"Crosscorrelation Coefficient ()",
short=r"$\mathrm{\mathsf{RHO_{HV}}}$",
),
KDP=dict(
ticks=np.array(
[-0.5, -0.1, 0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 2.0, 3.0]
),
cmap=cmap_prabhakar,
name=r"Specific Differental Phase (°/km)",
short=r"$\mathrm{\mathsf{K_{DP}}}$",
),
TEMP=dict(contours=[5, 0, -5, -10, -15, -20, -25, -30, -35, -40]),
PHI=dict(
ticks=np.array([0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90]),
cmap=cmap_prabhakar,
name=r"Differental Phase (°)",
short=r"$\mathrm{\mathsf{\Phi_{DP}}}$",
),
HMC=dict(
norm=mpl.colors.BoundaryNorm(np.arange(-0.5, 11 + 0.6, 1), 12),
ticks=np.arange(0, 11 + 1),
bounds=np.arange(-0.5, 11 + 0.6, 1),
cmap=mpl.colors.ListedColormap(
[
"LightBlue",
"Blue",
"Lime",
"Black",
"Red",
"Yellow",
"Fuchsia",
"LightPink",
"Cyan",
"Gray",
"DarkOrange",
"White",
]
),
name=r"HMC - Zrnic et al. 2001 ",
short=r"$\HMC_{Z}$",
long_name=r"Hydrometeorclass",
labels=[
"Light Rain",
"Moderate Rain",
"Heavy Rain",
"Large Drops",
"Hail",
"Rain/Hail",
"Graupel/Hail",
"Dry Snow",
"Wet Snow",
"H Crystals",
"V Crystals",
"No Rain",
],
labels_short=[
"LR",
"MR",
"HR",
"LD",
"HL",
"RH",
"GH",
"DS",
"WS",
"HC",
"VC",
"NR",
],
),
)
[docs]
def plot_moment(mom, *, ax=111, **kwargs):
sname = mom.attrs.get("standard_name", "")
if "equivalent_reflectivity_factor" in sname:
ticks = visdict14["ZH"]["ticks"]
elif "differential_reflectivity" in sname:
ticks = visdict14["ZDR"]["ticks"]
elif "correlation_coefficient" in sname:
ticks = visdict14["RHO"]["ticks"]
elif "radar_differential_phase" in sname:
ticks = visdict14["PHI"]["ticks"]
elif "radar_specific_differential_phase" in sname:
ticks = visdict14["KDP"]["ticks"]
else:
ticks = 14
cmap = None
ticks = kwargs.pop("ticks", ticks)
cmap = kwargs.pop("cmap", "HomeyerRainbow")
pm = mom.wrl.vis.plot(ax=ax, levels=ticks, cmap=cmap, **kwargs)
if ax == 111:
ax = plt.gca()
if mom["sweep_mode"] != "rhi":
maxrange = mom.range.max().values
ranges = [25e3, 50e3, 75e3, 100e3, 125e3, 150e3]
ranges = [rng for rng in ranges if rng <= maxrange + 1000]
wrl.vis.plot_ppi_crosshair((0, 0, 0), ranges=ranges, ax=ax)
if not ax.get_subplotspec().get_topmost_subplotspec().is_first_col():
ax.set_ylabel("")
if not ax.get_subplotspec().get_topmost_subplotspec().is_last_row():
ax.set_xlabel("")
dim0 = mom.wrl.util.dim0()
secs = "<M8[s]"
t0 = mom.time.load().min(dim0).values.astype(secs)
ax.set_title(f"{mom.name} - {t0}")
return pm
[docs]
def plot_2by2(swp, moments=None, **kwargs):
if moments is None:
moments = ["DBTH", "DBZH", "ZDR", "RHOHV"]
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(13, 10))
for i, mom in enumerate(moments):
plot_moment(swp[mom], ax=axes.flat[i], **kwargs)
plt.tight_layout()
[docs]
def plot_noise_hist(hist):
fig = plt.figure(figsize=(18, 18))
g = hist.histogram_RHO_SNR.plot(
col="noise",
col_wrap=5,
levels=30,
vmin=1,
vmax=None,
cmap="turbo",
norm=LogNorm(),
)
idxmin = hist.stdev.argmin()
for i, ax in enumerate(g.axes.flat):
if not ax.has_data():
continue
h = hist.isel(noise=i).histogram_RHO_SNR
h.idxmax("RHO_bin").plot(ax=ax, c="k")
if i == idxmin:
c = "r"
else:
c = "k"
ax.set_title(
f"NL {np.round(h.noise.values, decimals=1)} \n STD {np.around(h.stdev.values, decimals=5)}",
c=c,
)
g.fig.subplots_adjust(wspace=0.35, hspace=0.8, right=0.8)
fig.tight_layout()
def make_patch_spines_invisible(ax):
ax.set_frame_on(True)
ax.patch.set_visible(False)
for _, sp in ax.spines.items():
sp.set_visible(False)
def set_spine_direction(ax, direction):
if direction in ["right", "left"]:
ax.yaxis.set_ticks_position(direction)
ax.yaxis.set_label_position(direction)
elif direction in ["top", "bottom"]:
ax.xaxis.set_ticks_position(direction)
ax.xaxis.set_label_position(direction)
else:
raise ValueError(f"Unknown Direction: {direction}")
ax.spines[direction].set_visible(True)
[docs]
def create_lineplot(fig, subplot=111, xlabel=None, ylabel=None):
if xlabel is None:
xlabel = "Range Bins"
if ylabel is None:
ylabel = ""
host = fig.add_subplot(subplot)
ax = host
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
return ax
[docs]
def add_axis(host, ylabel=None, pos=1.0):
if ylabel is None:
ylabel = ""
ax = host.twinx()
host.spines["right"].set_visible(False)
ax.spines["right"].set_position(("axes", pos))
make_patch_spines_invisible(ax)
set_spine_direction(ax, "right")
ax.set_ylabel(ylabel)
return ax
def ufunc_wrapper(data, **kwargs):
func = kwargs.pop("func")
return func(data, **kwargs)
# phase processing
[docs]
def filter_data(obj, medwin):
dim0 = obj.wrl.util.dim0()
out = xr.apply_ufunc(
ufunc_wrapper,
obj,
input_core_dims=[[dim0, "range"]],
output_core_dims=[[dim0, "range"]],
output_dtypes=["float"],
dask="parallelized",
vectorize=True,
kwargs=dict(kernel_size=[1, medwin], func=scipy.signal.medfilt2d),
dask_gufunc_kwargs=dict(allow_rechunk=True),
)
return out
[docs]
def gauss_kernel(width, sigma):
dirac = np.zeros(width)
dirac[int(width / 2)] = 1
return scipy.ndimage.gaussian_filter1d(dirac, sigma=sigma)
def convolve(data, kernel=None, mode="same"):
mask = np.isnan(data)
out = np.convolve(np.where(mask, 0, data), kernel, mode=mode) / np.convolve(
~mask, kernel, mode=mode
)
return out
[docs]
def smooth_data(obj, kernel):
out = xr.apply_ufunc(
ufunc_wrapper,
obj,
input_core_dims=[["range"]],
output_core_dims=[["range"]],
output_dtypes=["float"],
dask="parallelized",
vectorize=True,
kwargs=dict(kernel=kernel, func=convolve),
dask_gufunc_kwargs=dict(allow_rechunk=True),
)
return out
def get_peaks(da, height=0):
def process_peaks(arr, height=0):
# Apply find_peaks
arr = arr.copy()
peaks, _ = scipy.signal.find_peaks(arr, height=height)
try:
peak = peaks[0]
except IndexError:
peak = 0
return peak
return xr.apply_ufunc(
process_peaks,
da,
input_core_dims=[["PHIDP_bin"]],
output_core_dims=[[]],
dask="parallelized",
vectorize=True,
kwargs=dict(height=height),
dask_gufunc_kwargs=dict(allow_rechunk=True),
)
[docs]
def phase_offset(phioff, rng=3000.0):
"""Calculate Phase offset.
Parameter
---------
phioff : xarray.DataArray
differential phase array
Keyword Arguments
-----------------
rng : float
range in m to calculate system phase offset
Return
------
start_range : xarray.DataArray
DataArray with start range values
off : xarray.DataArray
DataArray with phase offset values
"""
range_step = np.diff(phioff.range)[0]
nprec = int(rng / range_step)
if nprec % 2:
nprec += 1
# create binary array
phib = xr.where(np.isnan(phioff), 0, 1)
# take nprec range bins and calculate sum
phib_sum = phib.rolling(range=nprec, center=True).sum(skipna=True)
# get start range of first N consecutive precip bins
start_range = phib_sum.idxmax(dim="range") - nprec // 2 * np.diff(phib_sum.range)[0]
# get range of first non-nan value per ray
# start_range = (~np.isnan(phioff)).idxmax(dim='range', skipna=True)
# add range
stop_range = start_range + rng
# get phase values in specified range
off = phioff.where(
(phioff.range >= start_range.load()) & (phioff.range <= stop_range.load()),
drop=True,
)
# calculate nan median over range
off = off.median(dim="range", skipna=True)
return xr.Dataset(
dict(PHIDP_OFFSET=off, start_range=start_range, stop_range=stop_range)
)
[docs]
def phase_zphi(phi, rng=1000.0, **kwargs):
range_step = np.diff(phi.range)[0]
nprec = int(rng / range_step)
if nprec % 2 == 0:
nprec += 1
offset = nprec // 2 * np.diff(phi.range)[0]
offset_idx = nprec // 2
# create binary array
phib = xr.where(np.isnan(phi), 0, 1)
# take nprec range bins and calculate sum
min_periods = kwargs.get("min_periods", nprec)
phib_sum = phib.rolling(range=nprec, center=True, min_periods=min_periods).sum(
skipna=True
)
start_range = phib_sum.idxmax(dim="range") - offset
start_range_idx = phib_sum.argmax(dim="range") - offset_idx
first = phi.where(
(phi.range >= start_range.load()) & (phi.range <= start_range.load() + rng),
).median("range", skipna=True)
stop_range = phib_sum[..., ::-1].idxmax(dim="range") - offset
stop_range_idx = (
len(phib_sum.range) - (phib_sum[..., ::-1].argmax(dim="range") - offset_idx) - 2
)
last = phi.where(
(phi.range >= stop_range.load() - rng) & (phi.range <= stop_range.load()),
).median("range", skipna=True)
return xr.Dataset(
dict(
phib=phib_sum,
offset=offset,
offset_idx=offset_idx,
start_range=start_range,
stop_range=stop_range,
first=first,
first_idx=start_range_idx,
last=last,
last_idx=stop_range_idx,
dphi=last - first,
)
)
# Specific attenutation attributes
spec_att = dict(
units="dB/km",
standard_name="specific_attenuation_h",
long_name="Specific attenuation H",
short_name="AH",
)
[docs]
def zphi(phi, zh, **kwargs):
rng = kwargs.get("rng", 1000.0)
cphase = phase_zphi(phi, rng=rng)
if "alpha" in cphase.coords:
cphase = cphase.reset_coords("alpha", drop=True)
dphi = cphase.dphi.where(cphase.dphi >= 0, 1e-4)
alphax = kwargs.get("alphax", np.arange(0.1, 0.35, 0.01))
alphax = xr.DataArray(alphax, coords={"alpha": alphax}, dims="alpha")
bx = kwargs.get("bx", 0.78)
# need to expand alphax to dphi-shape
fdphi = 10 ** (0.1 * bx * alphax * dphi) - 1
# zhraw = zh.copy()#where((swp.CMAP2==0) & valid & (swp.RHOHV > 0.9))
zhraw = zh.where((zh.range > cphase.start_range) & (zh.range < cphase.stop_range))
zax = zhraw.wrl.trafo.idecibel().fillna(1e-4)
za = zax**bx
# set masked to zero for integration
za_zero = za.fillna(0)
# integrate
iza_x = 0.46 * bx * za_zero.pipe(x_cumtrapz)
iza = iza_x.max("range") - iza_x
iza_fdphi = iza / fdphi
# integral from r1 to r2
iza_first = iza_fdphi[..., cphase.first_idx.compute(), :]
ah = za / (iza_first + iza)
ah.attrs = spec_att
return ah
# to wradlib
def phidp_from_kdp(da):
"""Derive PHIDP from KDP.
Parameter
---------
da : xarray.DataArray
array with specific differential phase data
winlen : int
size of window in range dimension
Return
------
phi : xarray.DataArray
DataArray with differential phase values
"""
dr = da.range.diff("range").median("range").values / 1000.0
print("range res [km]:", dr)
return (
xr.apply_ufunc(
scipy.integrate.cumulative_trapezoid,
da,
input_core_dims=[["range"]],
output_core_dims=[["range"]],
dask="parallelized",
kwargs=dict(dx=dr, initial=0.0, axis=-1),
)
* 2
)
[docs]
def xr_rolling(da, window, window2=None, method="mean", min_periods=2, **kwargs):
"""Apply rolling function `method` to 2D datasets
Parameter
---------
da : xarray.DataArray
array with data to apply rolling function
window : int
size of window in range dimension
Keyword Arguments
-----------------
window2 : int
size of window in azimuth dimension
method : str
function name to apply
min_periods : int
minimum number of valid bins
**kwargs : dict
kwargs to feed to rolling function
Return
------
da_new : xarray.DataArray
DataArray with applied rolling function
"""
prng = window // 2
srng = slice(prng, -prng)
da_new = da.pad(range=prng, mode="reflect", reflect_type="odd")
dim = dict(range=window)
isel = dict(range=srng)
dim0 = da.wrl.util.dim0()
if window2 is not None:
paz = window2 // 2
saz = slice(paz, -paz)
da_new = da_new.pad({dim0: paz}, mode="wrap")
dim.update({dim0: window2})
isel.update({dim0: saz})
rolling = da_new.rolling(dim=dim, center=True, min_periods=min_periods)
da_new = getattr(rolling, method)(**kwargs)
da_new = da_new.isel(**isel)
return da_new
# Hydrometeor Classification
# todo check with wradlib
# def msf_index_indep_xarray(msf_ds, obs):
# def wrap_digitize(data, bins=None):
# return np.digitize(data, bins)
# idp = msf_ds.idp.values
# bins = np.append(idp, idp[-1] + (idp[-1] - idp[-2]))
# idx = (
# xr.apply_ufunc(
# wrap_digitize,
# obs,
# dask="parallelized",
# kwargs=dict(bins=bins),
# output_dtypes=["i4"],
# )
# - 1
# )
# # select bins
# idx = xr.where((idx >= 0) & (idx < bins.shape[0] - 2), idx, 0)
# return msf_ds.isel(idp=idx)
# def trapezoid(msf, obs):
# ones = (obs >= msf[..., 1]) & (obs <= msf[..., 2])
# ((obs < msf[..., 0]) | (obs > msf[..., 3]))
# lower = (obs >= msf[..., 0]) & (obs < msf[..., 1])
# higher = (obs > msf[..., 2]) & (obs <= msf[..., 3])
# obs_lower = obs - msf[..., 0]
# msf_lower = msf[..., 1] - msf[..., 0]
# low = obs_lower / msf_lower
# obs_higher = obs - msf[..., 3]
# msf_higher = msf[..., 2] - msf[..., 3]
# high = obs_higher / msf_higher
# ret = xr.zeros_like(obs) # * np.nan
# # ret = ret.where(zeros, 0).where(ones, 1).where(lower, low).where(higher, high)
# ret = xr.where(ones, 1, ret)
# ret = xr.where(lower, low, ret)
# ret = xr.where(higher, high, ret)
# return ret # .where(ret == np.nan, 0)
# def fuzzify(msf_ds, hmc_ds, msf_obs_mapping):
# fuzz_ds = xr.Dataset()
# for mf, hm in msf_obs_mapping.items():
# obs = hmc_ds[hm]
# msf = msf_ds[mf]
# fuzz_ds = fuzz_ds.assign({mf: trapezoid(msf, obs)})
# return fuzz_ds.transpose("hmc", ...)
# def probability(data, weights):
# """Calculate probability of hmc-class for every data bin.
# Parameters
# ----------
# data : xarray.Dataset
# Dataset containing containing the membership probability values.
# weights : :class:`numpy:numpy.ndarray`
# Array of length (observables) containing the weights for
# each observable.
# Returns
# -------
# out : xarray.DataArray
# Array containing weighted hmc-membership probabilities.
# """
# out = data.to_array(dim="obs")
# w = weights.to_array(dim="obs")
# out = (out * w).sum("obs") / w.sum("obs")
# return out # .transpose("hmc", ...)
# def classify(data, threshold=0.0):
# """Calculate probability of hmc-class for every data bin.
# Parameters
# ----------
# data : np.ndarray
# Array which is of size (hmc-class, data.shape), containing the
# weighted hmc-membership probability values.
# Keyword Arguments
# -----------------
# threshold : float
# Threshold value where probability is considered no precip,
# defaults to 0
# Returns
# -------
# out : xr.DataArray
# DataArray containing containing probability scores.
# No precip is added on the top.
# """
# # handle no precipitation
# nop = xr.where(data.sum("hmc") / len(data.hmc) <= threshold, 1, 0)
# nop = nop.assign_coords({"hmc": "NP"}).expand_dims(dim="hmc", axis=-1)
# return xr.concat([data, nop], dim="hmc")
def calc_gpm_height(data):
# Function for height calculation from bins.....
# GPM DPR ATBD Level 2 2017 Awaka et al S. 21
def _calc_height(bin_range_number):
rangeBinSize = 125.0
binEllipsoid = 176.0
bin_range_number = data["binBBBottom"]
elipsoid_offset = data["ellipsoidBinOffset"]
localZenithAngle = data["localZenithAngle"]
calc_1 = (binEllipsoid - bin_range_number) * rangeBinSize + elipsoid_offset
calc_h = calc_1 * np.cos(np.deg2rad(localZenithAngle))
return calc_h
data = data.assign(dpr_bbb=_calc_height(data["binBBBottom"]))
data = data.assign(dpr_bbt=_calc_height(data["binBBTop"]))
data = data.assign(dpr_cfb=_calc_height(data["binClutterFreeBottom"]))
return data
def read_gpm_file(fname, proj_stereo):
sr_data = wrl.io.open_gpm_dataset(fname, group="NS")
sr_data = sr_data.set_coords(["Longitude", "Latitude"])
sr_data = xr.decode_cf(sr_data)
sr_data.attrs["platform"] = "GPM"
sr_data.attrs = {}
# calculate bb and clutter heights
sr_data = calc_gpm_height(sr_data)
# assign x,y,z coordinates in WGS84
sr_data = sr_data.assign_coords(
x=sr_data.Longitude, y=sr_data.Latitude, z=sr_data.dpr_cfb
) # .load()
# Reproject to RADOLAN
proj_wgs = osr.SpatialReference()
proj_wgs.ImportFromEPSG(4326)
sr_data = sr_data.wrl.georef.reproject(src_crs=proj_wgs, trg_crs=proj_stereo)
# parallax correction
sr_data = sr_data.wrl.georef.correct_parallax(125.0)
return sr_data
def get_radolan_path(t, pr):
return (
f"/automount/radar/dwd/{pr}/{t:%Y}/{t:%Y}-{t:%m}/"
f"{t:%Y}-{t:%m}-{t:%d}/raa01-{pr}_10000-{t:%y%m%d%H%M}-dwd---bi*"
)
def read_radolan(gpm_time):
ts = pd.Timestamp(gpm_time.values)
ts = ts.round(freq="5T")
rx_path = get_radolan_path(ts, "rx")
ry_path = get_radolan_path(ts, "ry")
ry_ds = xr.open_dataset(glob.glob(ry_path)[0], engine="radolan")
ry_ds.RY.attrs["units"] = ry_ds.RY.attrs.pop("unit")
rx_ds = xr.open_dataset(glob.glob(rx_path)[0], engine="radolan")
rx_ds.RX.attrs["units"] = rx_ds.RX.attrs.pop("unit")
n_ds = xr.open_dataset(
"/automount/ftp/radar/wradlib-data/radolan/alt/radolan_geo.nc"
)
rdata = xr.merge([ry_ds, rx_ds, n_ds])
return rdata
def get_gpm_in_radolan(gr_data, sr_data):
from skimage import measure
bingrid = xr.where(gr_data.RY >= 0, 1, 0)
bingrid = bingrid.pad(pad_width=dict(x=1, y=1), constant_values=0)
# Contouren des Bin-grids erstellen
contours = measure.find_contours(
bingrid.values, 0, positive_orientation="high", fully_connected="high"
)
contours.sort(key=lambda x: len(x), reverse=True)
radolan_contour = contours[0]
cx = radolan_contour[:, 1].astype(int)
cy = radolan_contour[:, 0].astype(int)
# pad_coords
xx = gr_data.x.pad(pad_width=dict(x=1), mode="reflect", reflect_type="odd")
yy = gr_data.y.pad(pad_width=dict(y=1), mode="reflect", reflect_type="odd")
xx1 = xx.isel(x=cx)
yy1 = yy.isel(y=cy)
# GPM Coords GeoDataFrame
points_gpm = gpd.GeoDataFrame(
geometry=gpd.points_from_xy(sr_data.x.values.ravel(), sr_data.y.values.ravel())
)
# Radolan Shapely Polygon
xy = np.vstack((xx1.values.ravel(), yy1.values.ravel())).transpose()
radolan_polygon = shapely.Polygon(xy)
# Get point index
gpm_idx = points_gpm.within(radolan_polygon)
# to 1D and select per index
sr_data_1d = sr_data.stack(onedim=("nscan", "nray")).sel(
onedim=gpm_idx.to_xarray().rename(index="onedim")
)
# add xy coords
xy = xr.concat([sr_data_1d.x, sr_data_1d.y], "xyc") # .transpose()
sr_data_1d = sr_data_1d.assign_coords(xy=xy)
# reset Multi-Index
sidx = sr_data_1d.reset_index("onedim").assign_coords(
onedim=np.arange(len(sr_data_1d.onedim))
)
return sidx.transpose("onedim", ...)
def get_radolan_in_gpm(sr_data, gr_data):
dpr_contour_x = dpr_swath_contour(sr_data.x)
dpr_contour_y = dpr_swath_contour(sr_data.y)
# create dpr polygon of dpr xy contours
dpr_xy_poly = np.vstack((dpr_contour_x.ravel(), dpr_contour_y.ravel())).transpose()
gpm_polygon = shapely.Polygon(dpr_xy_poly)
radolan_xy = xr.concat(xr.broadcast(gr_data.x, gr_data.y), dim="xy")
radolan_xy.name = "radolan_coords"
radolan_xy = radolan_xy.stack(onedim=("x", "y")).transpose(..., "xy")
points_radolan = gpd.GeoDataFrame(
geometry=gpd.points_from_xy(radolan_xy[..., 0], radolan_xy[..., 1])
)
radolan_idx = points_radolan.within(gpm_polygon)
rdata_1d = gr_data.stack(onedim=("x", "y")).sel(
onedim=radolan_idx.to_xarray().rename(index="onedim")
)
# reset Multi-Index
ridx = rdata_1d.reset_index("onedim").assign_coords(
onedim=np.arange(len(rdata_1d.onedim))
)
return ridx
def dpr_antenna_weighting(r):
"""
Funktion: explained in
Watters, D., A. Battaglia, K. Mroz, and F. Tridon, 0:
Validation of the GPM Version-5 Surface Rainfall
Products over Great Britain and Ireland.
J. Hydrometeor., 0, https://doi.org/10.1175/JHM-D-18-0051.1
and
Mroz, K., A. Battaglia, T. J. Lang, D. J. Cecil,
S. Tanelli, and F. Tridon, 2017: Hail-detection
algorithm for the GPM Core Observatory satellite sensors.
Journal of Applied Meteorology and Climatology,56 (7), 1939–1957
Input:
------
r ::: Distance to the center point of the footprint in m
Output:
-------
omega ::: Weighting for the point with the distance
"""
r = r / 1000.0
omega = np.exp(-((r / 2.5) ** 2.0) * np.log10(4.0))
return omega
def _weighted_sum(patch, weights, wsum, dim):
is_refl = patch.attrs["units"] == "dBZ" if "units" in patch.attrs else False
if is_refl:
patch = patch.wrl.trafo.idecibel()
patch = (patch * weights).sum(dim, skipna=True) / wsum
if is_refl:
patch = patch.wrl.trafo.decibel()
return patch
def weighted_sum(patch, weights, dim):
patch = patch.copy()
wsum = weights.sum(dim=dim, skipna=True)
for vname, dvar in patch.data_vars.items():
if dim in dvar.dims and dvar.dtype != object:
patch[vname] = _weighted_sum(dvar, weights, wsum, dim)
return patch
def radolan_gpm_matching(ridx, sidx, k=32, weight_gpm=True):
# build tree
tree = KDTree(np.vstack([ridx.x, ridx.y]).T, balanced_tree=False)
# query tree for 32 nearest neighbours
# 32 maximal possible ry pixel in dpr footprint
dists, ix = tree.query(sidx.xy, k=k, workers=-1)
# assign radolan data to dataset
sidx1 = sidx.assign(ridx.sel(onedim=xr.DataArray(ix, dims=("onedim", "bins"))))
# assign distances
sidx1 = sidx1.assign(dists=xr.DataArray(dists, dims=("onedim", "bins")))
# calculate weighting and summation
w = dpr_antenna_weighting(sidx1.dists)
wsum = weighted_sum(sidx1, w, "bins")
w = dpr_antenna_weighting(np.abs(wsum.rhmin - wsum.zp))
wsum = wsum.assign(weights=w)
# select gpm disc's
sr_match = wsum.where(
(wsum.zp >= (wsum.rhmin - wsum.beam_width / 2.0))
& (wsum.zp <= (wsum.rhmin + wsum.beam_width / 2.0))
)
# mean on gpm column
if weight_gpm:
res = weighted_sum(
sr_match.reset_coords(["xp", "yp", "zp"]), sr_match.weights, "nbin"
)
else:
res = sr_match.reset_coords(["xp", "yp", "zp"]).mean(dim="nbin", skipna=True)
return res
def read_gr_sweep(gr_file_path, offset_z=0, offset_phi=0, offset_zdr=0):
"""
Read and organize the GR data for one sweep.The following code reads data
in ODIM H5 format. If your GR data is in some other format
respective adaptions are needed.
# Input:
# ------
gr_file_path ::: path to graund radar sweep data
offset_z ::: Refl. offset in dB
offset_phi ::: System PHI
offset_zdr ::: Diff. Refl offst in dB
# Output:
# ------
ds0_masked ::: Sweep Data with PPB Quality Index
More description in the data itself.
"""
def noise_correction(ds, noise_level):
# noise calculations
snr = (
ds.DBZH
- 20 * np.log10(ds.range * 0.001)
- noise_level
- 0.033 * ds.range / 1000
)
rho = ds.RHOHV * np.sqrt(1.0 + 1.0 / 10.0 ** (snr * 0.1))
return ds.assign({"SNR": snr, "RHOHV": rho})
def xr_rolling(da, window, method="mean"):
da_new = da.pad(range=window // 2, mode="reflect", reflect_type="odd")
da_new = getattr(da_new.rolling(range=window, center=True), method)()
da_new = da_new.isel(range=slice(window // 2, -window // 2))
return da_new
# read one ppi data
# umstellen open_odim_dataset
ds0 = xr.open_dataset(gr_file_path, engine="gamic", group="sweep_0")
# Offsets and Noise Processing
# ------------------------------
# keep uncorrected RHOHV
ds0 = ds0.assign({"URHOHV": ds0.RHOHV})
# add the offsets
ds0["DBZH"] -= offset_z
ds0["PHIDP"] -= offset_phi
ds0["ZDR"] -= offset_zdr
# RHOHV noise correction
# ds0 = ds0.pipe(noise_correction, -22)
# ds0 = ds0.pipe(noise_correction, -22)
# georeference data
ds0 = ds0.wrl.georef.georeference()
## Phase Processing
## -----------------
window = 13
phi_ma = ds0.PHIDP.pipe(xr_rolling, window, method="median")
# mask phi
phi_masked = phi_ma.where((ds0.RHOHV >= 0.95) & (ds0.DBZH >= 0.0))
# not used for moment
# todo: check 2023-11-06
np.isfinite(phi_masked).argmax(dim="range")
phi_masked.sizes["range"] - np.isfinite(
phi_masked.sortby(["range"], ascending=False)
).argmax(dim="range")
# nur bei neuem wradlib
kdp = phi_masked.wrl.dp.kdp_from_phidp(winlen=21, min_periods=3)
phidp = kdp.pipe(phidp_from_kdp)
phi = ds0.PHIDP.where((ds0.RHOHV >= 0.95) & (ds0.DBZH >= 0.0))
phidp2, kdp2 = phi.wrl.dp.phidp_kdp_vulpiani(
winlen=15, th1=-0.5, th2=5, th3=-80, min_periods=3
)
ds0 = ds0.assign({"PHI_SMOOTH": phidp})
ds0 = ds0.assign({"PHI_MASKED": phi_masked})
ds0 = ds0.assign({"KDP1_NEW": kdp})
ds0 = ds0.assign({"PHI1_NEW": phidp})
ds0 = ds0.assign({"KDP2_NEW": kdp2})
ds0 = ds0.assign({"PHI2_NEW": phidp2})
ds0_masked = ds0.where(ds0.RHOHV >= 0.7)
ds0_masked = ds0_masked.assign({"height": ds0.z})
sitecoords = (
ds0_masked.longitude.data,
ds0_masked.latitude.data,
ds0_masked.altitude.data,
)
# BB and Q
nrays = 360 # number of rays
nbins = ds0_masked.range.shape[0] # number of range bins
el = ds0_masked.elevation.data[0] # vertical antenna pointing angle (deg)
bw = 1.0 # half power beam width (deg)
range_res = (
ds0_masked.range.data[1] - ds0_masked.range.data[0]
) # range resolution (meters)
r = np.arange(nbins) * range_res
beamradius = wrl.util.half_power_radius(r, bw)
coord = wrl.georef.sweep_centroids(nrays, range_res, nbins, el)
coords = wrl.georef.spherical_to_proj(
coord[..., 0], coord[..., 1], coord[..., 2], sitecoords
)
lon = coords[..., 0]
lat = coords[..., 1]
alt = coords[..., 2]
polcoords = coords[..., :2]
rlimits = (lon.min(), lat.min(), lon.max(), lat.max())
ds = wrl.io.open_raster(
"/automount/ags/velibor/radcal/scripts/boxpol_processing/bonn_dem.tif"
)
rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(
ds, nodata=-32768.0
)
# Clip the region inside our bounding box
ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
rastercoords = rastercoords[ind[1] : ind[3], ind[0] : ind[2], ...]
rastervalues = rastervalues[ind[1] : ind[3], ind[0] : ind[2]]
# Map rastervalues to polar grid points
polarvalues = wrl.ipol.cart_to_irregular_spline(
rastercoords, rastervalues, polcoords, order=3, prefilter=False
)
PBB = wrl.qual.beam_block_frac(polarvalues, alt, beamradius)
PBB = np.ma.masked_invalid(PBB)
CBB = wrl.qual.cum_beam_block_frac(PBB)
def QBBF(CBB):
"""
Quality Index for PBB
Crisologo 2018
"""
Qbbf = np.zeros_like(CBB) * np.nan
Qbbf[CBB < 0.1] = 1
Qbbf[CBB > 0.5] = 0
Qbbf[(CBB <= 0.5) & (CBB >= 0.1)] = (
1 - (CBB[(CBB <= 0.5) & (CBB >= 0.1)] - 0.1) / 0.4
)
return Qbbf
# ds0_masked = ds0_masked.assign({'QBBF': QBBF(CBB)})
ds0_masked = ds0_masked.assign(
{"QBBF": (["azimuth", "range"], QBBF(CBB), ds0_masked.ZDR.attrs)}
)
return ds0_masked
def get_bbox(gr_data):
"""
Calculation of bounding box around the ground radar
# Input:
# ------
gr_data ::: graund radar sweep data
# Output:
# ------
bbox ::: Coords for bounding box
"""
# number of rays in gr sweep
nray_gr = gr_data.rays.data.shape[0]
# gate length (meters)
dr_gr = gr_data.range.data[1] - gr_data.range.data[0]
# number of gates in gr beam
ngate_gr = gr_data.rays.data.shape[1]
# elevation of sweep (degree)
elev_gr = gr_data.elevation.data[0]
# Longitude of GR
lon0_gr = gr_data.longitude.data
# Latitude of GR
lat0_gr = gr_data.latitude.data
# Altitude of GR (meters)
alt0_gr = gr_data.altitude.data
# spherical to proj
coord = wrl.georef.sweep_centroids(nray_gr, dr_gr, ngate_gr, elev_gr)
coords = wrl.georef.spherical_to_proj(
coord[..., 0],
np.degrees(coord[..., 1]),
coord[..., 2],
(lon0_gr, lat0_gr, alt0_gr),
)
lon = coords[..., 0]
lat = coords[..., 1]
bbox = wrl.zonalstats.get_bbox(lon, lat)
return bbox
def get_grid_polygons(ds):
x = ds.x
for i in reversed(range(x.ndim)):
x = xr.plot.utils._infer_interval_breaks(x, axis=i)
y = ds.y
for i in reversed(range(y.ndim)):
y = xr.plot.utils._infer_interval_breaks(y, axis=i)
coords = np.stack([x, y], axis=-1)
ll = np.dstack([coords[0:-1, 0:-1], ds.z.values[..., None]])
ul = np.dstack([coords[0:-1, 1:], ds.z.values[..., None]])
ur = np.dstack([coords[1:, 1:], ds.z.values[..., None]])
lr = np.dstack([coords[1:, 0:-1], ds.z.values[..., None]])
return np.stack([ll, ul, ur, lr, ll], axis=-2)
def get_sr_cut(swp, sr_data):
# dr_sr = 125.0
# Todo: only keep one method
# # Calculate equivalent earth radius
# wgs84 = wrl.georef.get_default_projection()
# # lon0_gr = gr_data.latitude.values
# re1 = swp.wrl.georef.get_earth_radius(sr=wgs84)
# a = wgs84.GetSemiMajor()
# b = wgs84.GetSemiMinor()
# Set up aeqd-projection gr-centered
# proj_rad = wrl.georef.projstr_to_osr(
# ("+proj=aeqd +lon_0={lon:f} " + "+lat_0={lat:f} +a={a:f} " + "+b={b:f}").format(
# lon=swp.longitude.values, lat=swp.latitude.values, a=a, b=b
# )
# )
# swp.wrl.georef.get_earth_radius(sr=proj_rad)
proj_rad = osr.SpatialReference()
proj_rad.ImportFromWkt(xd.georeference.get_crs(swp).to_wkt())
sr_data = sr_data.wrl.georef.reproject(
coords=dict(x="Longitude", y="Latitude"),
src_crs=wrl.georef.get_default_projection(),
trg_crs=proj_rad,
)
# calculate range and mask srdata
r = np.sqrt(sr_data.x**2 + sr_data.y**2)
sr_cut = sr_data.where(r < swp.gr.max(), drop=True)
return sr_cut # , proj_rad, re1
def calc_sat(sr_cut, swp, dr_sr): # , proj_rad, re1):
sr_cut = sr_cut.wrl.georef.correct_parallax(dr_sr)
proj_rad = osr.SpatialReference()
proj_rad.ImportFromWkt(xd.georeference.get_crs(swp).to_wkt())
re1 = swp.wrl.georef.get_earth_radius()
sr_cut = sr_cut.wrl.georef.xyz_to_spherical(crs=proj_rad)
bw_sr = 0.71
freq = 0
sr_cut = sr_cut.wrl.georef.dist_from_orbit(bw_sr, freq, re1)
rs = sr_cut["pr_dist"]
vol_sr = rs.wrl.qual.pulse_volume(dr_sr, bw_sr)
sr_cut = sr_cut.assign(vol_sr=vol_sr)
bw_gr = 1.0
vol_gr = wrl.qual.pulse_volume(swp.range, swp.range.diff("range").median(), bw_gr)
vol_gr = vol_gr.broadcast_like(swp.DBZH)
freq = 0
alpha = sr_cut["localZenithAngle"].isel(nfreq=freq, missing_dims="ignore")
# Rs
sr_bin_radius = (
0.5 * (1 + np.cos(np.radians(alpha))) * rs * np.tan(np.radians(bw_sr / 2.0))
)
sr_bin_radius.name = "sr_bin_radius"
# Ds
sr_bin_depth = dr_sr / np.cos(np.radians(alpha))
sr_bin_depth.name = "sr_bin_depth"
sr_cut = sr_cut.assign(sr_bin_radius=sr_bin_radius, sr_bin_depth=sr_bin_depth)
return sr_cut
def calculate_transform(sr_cut):
ref_sr = sr_cut.zFactorFinal.isel(nfreq=0, missing_dims="ignore").copy()
# ref_sr = sr_data['refl'].filled(np.nan)
ref_sr_ss = xr.zeros_like(ref_sr) * np.nan
ref_sr_sh = xr.zeros_like(ref_sr) * np.nan
a_s, a_h = (wrl.trafo.KuBandToX.snow, wrl.trafo.KuBandToX.hail)
# mask bb_ratio with bb_mask
ratio = sr_cut["bb_ratio"].where(sr_cut["bb_mask"])
ia = ratio >= 1
ib = ratio <= 0
im = (ratio > 0) & (ratio < 1)
im = xr.where(im, np.round(ratio * 10), 0).astype("int")
# adding ia/ib/im in sequential order to ref_sr_ss/ref_sr/sh
# above melting layer -> no mixing only snow/hail
ref_sr_ss = xr.where(
ia, ref_sr + wrl.util.calculate_polynomial(ref_sr.copy(), a_s[:, 1]), ref_sr_ss
)
ref_sr_sh = xr.where(
ia, ref_sr + wrl.util.calculate_polynomial(ref_sr.copy(), a_h[:, 1]), ref_sr_sh
)
# below the melting layer
ref_sr_ss = xr.where(
ib, ref_sr + wrl.util.calculate_polynomial(ref_sr.copy(), a_s[:, 0]), ref_sr_ss
)
ref_sr_sh = xr.where(
ib, ref_sr + wrl.util.calculate_polynomial(ref_sr.copy(), a_h[:, 0]), ref_sr_sh
)
# ref_sr_ss = xr.where(im, ref_sr + np.nan, ref_sr_ss)
# ref_sr_sh = xr.where(im, ref_sr + np.nan, ref_sr_sh)
ref_sr_ss = xr.where(im, ref_sr + np.nan, ref_sr)
ref_sr_sh = xr.where(im, ref_sr + np.nan, ref_sr)
return ref_sr_ss, ref_sr_sh
[docs]
def volume_matching(swp, sr_data, ref="DBZH", grrefth=15):
dr_sr = 125.0
swp = swp.copy().where(swp[ref] >= grrefth)
gr_poly = get_grid_polygons(swp)
# # Todo: only keep one method
# # Calculate equivalent earth radius
# wgs84 = wrl.georef.get_default_projection()
# # lon0_gr = gr_data.latitude.values
# re1 = swp.wrl.georef.get_earth_radius(sr=wgs84)
# a = wgs84.GetSemiMajor()
# b = wgs84.GetSemiMinor()
# # Set up aeqd-projection gr-centered
# proj_rad = wrl.georef.projstr_to_osr(
# ("+proj=aeqd +lon_0={lon:f} " + "+lat_0={lat:f} +a={a:f} " + "+b={b:f}").format(
# lon=swp.longitude.values, lat=swp.latitude.values, a=a, b=b
# )
# )
# swp.wrl.georef.get_earth_radius(sr=proj_rad)
# sr_data = sr_data.wrl.georef.reproject(
# coords=dict(x="Longitude", y="Latitude"),
# src_crs=wrl.georef.get_default_projection(),
# trg_crs=proj_rad,
# )
# # calculate range and mask srdata
# r = np.sqrt(sr_data.x**2 + sr_data.y**2)
# sr_cut = sr_data.where(r < swp.gr.max(), drop=True)
sr_cut = get_sr_cut(swp, sr_data)
precip_mask = xr.where(sr_cut.flagPrecip > 0, 1, 0)
sr_cut = calc_sat(sr_cut, swp, dr_sr)
sr_cut = sr_cut.wrl.qual.get_bb_ratio()
ref_sr_ss, ref_sr_sh = calculate_transform(sr_cut)
ptype = (sr_cut.typePrecip / 1e7).astype(np.int16)
# Combined conversion
ref_sr_com = xr.zeros_like(ref_sr_ss) * np.nan
ref_sr_com = xr.where(ptype == 1, ref_sr_ss, ref_sr_com)
ref_sr_com = xr.where(ptype == 2, ref_sr_sh, ref_sr_com)
sr_cut = sr_cut.assign(
ref_sr_com=ref_sr_com,
ref_sr_ss=ref_sr_ss,
ref_sr_sh=ref_sr_sh,
precip_mask=precip_mask,
)
## matching below
bw_gr = 1.0
swp = swp.assign(sweep_fixed_angle=swp.sweep_fixed_angle.median(skipna=True))
sr_match = sr_cut.reset_coords(["xp", "yp", "zp"])
sr_match = sr_match.where(sr_cut.precip_mask).where(
(sr_match.elevation >= (swp.sweep_fixed_angle - bw_gr / 2.0))
& (sr_match.elevation <= (swp.sweep_fixed_angle + bw_gr / 2.0))
)
# calculate nbeams
sr_match_mean = sr_match.mean("nbin", skipna=True)
xyzp_sr = xr.concat(
[sr_match_mean.xp, sr_match_mean.yp, sr_match_mean.zp], dim="ncoord"
).transpose(..., "ncoord")
xyzp_sr = (
xyzp_sr.stack(nbeams=("nscan", "nray"))
.transpose("nbeams", "ncoord", ...)
.dropna(
dim="nbeams",
)
)
sr_match2 = sr_match.stack(nbeams=("nscan", "nray")).where(xyzp_sr.nbeams)
# mean along nbin
sr_match_mean = sr_match.mean(dim="nbin", skipna=True)
sr_match2_mean = sr_match2.mean(dim="nbin", skipna=True)
# max along nbin, approximate Rs
# sr_match_max = sr_match.max(dim="nbin", skipna=True)
sr_match2_max = sr_match2.max(dim="nbin", skipna=True)
# sum along nbin
# approximate Ds and Vs
sr_match.sum(dim="nbin", skipna=True)
sr_match2.sum(dim="nbin", skipna=True)
## Create Zonal Data Source of matching SR profiles
proj_rad = osr.SpatialReference()
proj_rad.ImportFromWkt(xd.georeference.get_crs(swp).to_wkt())
xyp_sr = xr.concat([sr_match2_mean.xp, sr_match2_mean.yp], dim="xy").transpose(
..., "xy"
)
zds = wrl.io.VectorSource(xyp_sr.values, trg_crs=proj_rad)
## Create SR beam polygons using sr_bin_radius Rs
tmp_trg_lyr = zds.ds.GetLayer()
trg_poly = []
for i, feat in enumerate(tmp_trg_lyr):
geom = feat.GetGeometryRef().Buffer(sr_match2_max.sr_bin_radius.values[i])
poly = wrl.georef.ogr_to_numpy(geom)
trg_poly.append(poly)
zdp = wrl.zonalstats.ZonalDataPoly(
gr_poly[..., 0:2].reshape(-1, 5, 2),
trg=trg_poly,
crs=proj_rad,
)
zdp.dump_vector("m3d_zonal_poly")
obj3 = wrl.zonalstats.ZonalStatsPoly("m3d_zonal_poly")
# converted for stratiform
ref_sr2gr_m = sr_match2.ref_sr_com.copy()
ref_sr2gr_m = wrl.trafo.idecibel(ref_sr2gr_m)
ref_sr2gr_m = ref_sr2gr_m.mean(dim="nbin", skipna=True)
ref_sr2gr_m = wrl.trafo.decibel(ref_sr2gr_m)
ref_sr2gr_m.name = "mean_sr_reflectivity"
# mean = obj3.mean(wrl.trafo.idecibel(swp.DBZH).values.ravel())
gr_mean = wrl.trafo.idecibel(swp[ref]).values.ravel()
gr_mean = np.array(
[
np.nanmean(gr_mean[obj3.ix[i]])
for i in np.arange(len(obj3.ix))[~obj3.check_empty()]
]
)
obj3.zdata.trg.geo["mean"] = wrl.trafo.decibel(gr_mean)
ref_gr_m = obj3.zdata.trg.geo["mean"].to_xarray()
ref_gr_m.name = "mean_gr_reflectivity"
ref_gr_m = ref_gr_m.rename(index="nbeams")
pbb_gr_i = swp.QBBF.values.ravel()
pbbgr1a = np.array(
[
np.nanmin(pbb_gr_i[obj3.ix[i]])
for i in np.arange(len(obj3.ix))[~obj3.check_empty()]
]
)
qbbf = xr.zeros_like(ref_gr_m)
qbbf.data = pbbgr1a
qbbf.name = "pbb_quality_index"
# qbbf = xr.zeros_like(ref_gr_m)
# qbbf.data = pbbgr1a
#########################
# pia_gr_i = swp.QPIA.values.ravel()
# piagr1a = np.array(
# [
# np.nanmax(pia_gr_i[obj3.ix[i]])
# for i in np.arange(len(obj3.ix))[~obj3.check_empty()]
# ]
# )
# qpia = xr.zeros_like(ref_gr_m)
# qpia.data = piagr1a
# return xr.merge([ref_gr_m, ref_sr2gr_m, qbbf, xyp_sr])
xy = xyp_sr.reset_index("nbeams").reset_coords(drop=True)
ds = xr.merge(
[ref_gr_m, ref_sr2gr_m.reset_index("nbeams").reset_coords(drop=True), qbbf]
)
ds = ds.assign_coords(x=xy.isel(xy=0), y=xy.isel(xy=1))
return ds
# Todo: make height level parameterizable via keyword argument
# Todo: plot gr radar domain and satellite swath outlines
[docs]
def plot_unmatched_refl(sr_data, gr_data, level=-5):
"""
Plotting unmatched SR and GR Reflectivity
# Input:
# ------
gr_data ::: ground radar sweep data
sr_data ::: satellite data
# Output:
# ------
Plot of SR Refl. and GR Refl.
"""
plt.figure(figsize=(12, 4))
ax1 = plt.subplot(1, 2, 1)
pm = sr_data.zFactorCorrected[:, :, level].plot(
x="Longitude",
y="Latitude",
cmap="HomeyerRainbow",
vmin=0,
vmax=40,
add_labels=False,
add_colorbar=False,
)
cbar = plt.colorbar(pm, ax=ax1, ticks=np.arange(0, 42, 2))
cbar.set_label("SR Reflectivity (dBz)", fontsize=12)
cbar.ax.tick_params(labelsize=12)
ax1.set_xlabel("Longitude (°)", fontsize=12)
ax1.set_ylabel("Latitude (°)", fontsize=12)
ax1.tick_params(axis="both", labelsize=12)
ax1.set_title(sr_data.attrs["platform"], fontsize=12, loc="left")
ax1.grid(lw=0.25, color="grey")
ax1.set_xlim(sr_data.Longitude.min(), sr_data.Longitude.max())
ax1.set_ylim(sr_data.Latitude.min(), sr_data.Latitude.max())
ax2 = plt.subplot(1, 2, 2)
swp = gr_data
ll = swp.wrl.georef.spherical_to_proj()
swp = swp.assign_coords(lon=ll[..., 0], lat=ll[..., 1])
pm = swp.DBZH.plot(
x="lon",
y="lat",
cmap="HomeyerRainbow",
vmin=0,
vmax=40,
ax=ax2,
add_labels=False,
add_colorbar=False,
)
cbar = plt.colorbar(pm, ax=ax2, ticks=np.arange(0, 42, 2))
cbar.set_label("GR Reflectivity (dBz)", fontsize=12)
cbar.ax.tick_params(labelsize=12)
ax2.set_xlabel("Longitude (°)", fontsize=12)
ax2.set_ylabel("Latitude (°)", fontsize=12)
ax2.tick_params(axis="both", labelsize=12)
ax2.set_title("Ground Radar", fontsize=12, loc="left")
ax2.grid(lw=0.25, color="grey")
ax2.set_xlim(sr_data.Longitude.min(), sr_data.Longitude.max())
ax2.set_ylim(sr_data.Latitude.min(), sr_data.Latitude.max())
plt.tight_layout()
[docs]
def plot_matched_refl(vmatch):
"""
Plotting unmatched SR and GR Reflectivity
# Input:
# ------
gr_data ::: graund radar sweep data
sr_data ::: satellite data
# Output:
# ------
Plot of SR Refl. and GR Refl.
"""
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.scatter(
vmatch.x,
vmatch.y,
c=vmatch.mean_sr_reflectivity,
s=10,
cmap="HomeyerRainbow",
vmin=0,
vmax=40,
)
cbar = plt.colorbar(ticks=np.arange(0, 42, 2))
cbar.set_label("SR Reflectivity (dBz)", fontsize=12)
cbar.ax.tick_params(labelsize=12)
plt.xlabel("Longitude (°)", fontsize=12)
plt.ylabel("Latitude (°)", fontsize=12)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.title("Satellite Radar", fontsize=12, loc="left")
plt.grid(lw=0.25, color="grey")
plt.subplot(1, 2, 2)
plt.scatter(
vmatch.x,
vmatch.y,
c=vmatch.mean_gr_reflectivity,
s=10,
cmap="HomeyerRainbow",
vmin=0,
vmax=40,
)
cbar = plt.colorbar(ticks=np.arange(0, 42, 2))
cbar.set_label("GR Reflectivity (dBz)", fontsize=12)
cbar.ax.tick_params(labelsize=12)
plt.xlabel("x (km)", fontsize=12)
plt.ylabel("y (km)", fontsize=12)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.title("Ground Radar", fontsize=12, loc="left")
plt.grid(lw=0.25, color="grey")
plt.tight_layout()
def dpr_swath_contour(sr_grid):
"""
Function for dpr swath grid contour
## sr_grid[:, 0] # - complete left column
## sr_grid[-1,1:-1] # - upper edge without left and right corner
## sr_grid[:, -1][::-1] # - complete right column (backward)
## sr_grid[0,0:-1] # - lower edge (without right corner)
Input:
------
sr_grid ::: Space-borne grid x,y
Output:
-------
trg ::: Contour of the space-borne grid
"""
trg = np.r_[
sr_grid[:, 0], sr_grid[-1, 1:-1], sr_grid[:, -1][::-1], sr_grid[0, 0:-1][::-1]
]
return trg
def gpm_scan_time(filename, scanswath="NS"):
from netCDF4 import Dataset
pr_data = Dataset(filename, mode="r")
year = pr_data[scanswath]["ScanTime"].variables["Year"]
month = pr_data[scanswath]["ScanTime"].variables["Month"]
dayofmonth = pr_data[scanswath]["ScanTime"].variables["DayOfMonth"]
# dayofyear = pr_data['NS']['ScanTime'].variables['DayOfYear'][mask]
hour = pr_data[scanswath]["ScanTime"].variables["Hour"]
minute = pr_data[scanswath]["ScanTime"].variables["Minute"]
second = pr_data[scanswath]["ScanTime"].variables["Second"]
# secondofday = pr_data['NS']['ScanTime'].variables['SecondOfDay'][mask]
millisecond = pr_data[scanswath]["ScanTime"].variables["MilliSecond"]
date_array = zip(
year, month, dayofmonth, hour, minute, second, millisecond, strict=True
)
pr_time = np.array(
[dt.datetime(d[0], d[1], d[2], d[3], d[4], d[5], d[6]) for d in date_array]
)
# ttt = gpm_scan_time(gpm_file)
s_n = {"NS": 49, "HS": 24, "MS": 25}
time_array = np.array(s_n[scanswath] * [pr_time]).T
pr_data.close()
return time_array
def get_time_of_gpm2(gpm_time):
"""einfache mittlere zeit"""
iii = int(len(np.array(gpm_time["Year"])) / 2)
gpm_year = np.array(gpm_time["Year"])[iii]
gpm_month = np.array(gpm_time["Month"])[iii]
gpm_day = np.array(gpm_time["DayOfMonth"])[iii]
gpm_hour = np.array(gpm_time["Hour"])[iii]
gpm_min = np.array(gpm_time["Minute"])[iii]
gpm_sek = np.array(gpm_time["Second"])[iii]
gpm_dt = dt.datetime(
gpm_year, gpm_month, gpm_day, gpm_hour, gpm_min, gpm_sek
).strftime("%Y.%m.%d -- %H:%M:%S")
return gpm_dt
def hist_2d(
A,
B,
binsx=35,
binsy=35,
mini=1,
maxi=None,
cmap="turbo",
colsteps=30,
alpha=1,
fsize=15,
colbar=True,
):
"""
# Histogram 2d Quicklooks
# ------------------------
Plotting 2d Histogramm of two varibles
# Input
# -----
A,B ::: Variables
bins1, bins2 ::: x, y bins
mini, maxi ::: min and max
cmap ::: colormap
colsteps ::: number of cmap steps
alpha ::: transperency
fsize ::: fontsize
# Output
# ------
2D Histogramm Plot
::: Hist mode:::
absolute ::: absolute numbers
relative ::: relative numbers
relative_with_y ::: relative numbers of y levels
"""
from matplotlib.colors import LogNorm
# discret cmap
cmap = mpl.colormaps.get_cmap(cmap, colsteps)
# mask array
m = ~np.isnan(A) & ~np.isnan(B)
plt.hist2d(
A[m],
B[m],
bins=(binsx, binsy),
cmap=cmap,
norm=LogNorm(vmin=mini, vmax=maxi),
alpha=alpha,
)
if colbar is True:
cb = plt.colorbar(shrink=1, pad=0.01)
cb.set_label("number of samples", fontsize=fsize)
cb.ax.tick_params(labelsize=fsize)
plt.xticks(fontsize=fsize)
plt.yticks(fontsize=fsize)
[docs]
def set_sweep_mode_as_coords(ds):
ds = ds.set_coords(["sweep_mode"])
return ds
[docs]
def x_cumtrapz(da):
"""xarray wrapper for scipy.cumtrapz"""
dr = da.range.diff("range").median("range").values / 1000.0
print("range res [km]:", dr)
return xr.apply_ufunc(
scipy.integrate.cumulative_trapezoid,
da,
input_core_dims=[["range"]],
output_core_dims=[["range"]],
dask="parallelized",
kwargs=dict(dx=dr, initial=0.0, axis=-1),
)
[docs]
def plot_beamprop(r, theta, sitecoords, er, k):
beamradius = wrl.util.half_power_radius(r, 1.0)
beam_height = wrl.georef.bin_altitude(
r, theta[:, np.newaxis], sitecoords[2], re=er, ke=k
)
arc_distance = wrl.georef.bin_distance(
r, theta[:, np.newaxis], sitecoords[2], re=er, ke=k
)
def height_formatter(x, pos):
x = (x - 6370000) / 1000
fmt_str = f"{x:g}"
return fmt_str
def range_formatter(x, pos):
x = x / 1000.0
fmt_str = f"{x:g}"
return fmt_str
fig = plt.figure(figsize=(10, 6))
cgax, caax, paax = wrl.vis.create_cg(fig=fig, rot=0, scale=1)
# azimuth angle
# fix grid_helper
gh = cgax.get_grid_helper()
gh.grid_finder.grid_locator2._nbins = 40
gh.grid_finder.grid_locator2._steps = [1, 2, 4, 5, 10]
# calculate beam_height and arc_distance for ke=1
# means line of sight
bhe = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=1.0)
ade = wrl.georef.bin_distance(r, 0, sitecoords[2], re=er, ke=1.0)
nn0 = np.zeros_like(r)
# for nice plotting we assume earth_radius = 6370000 m
ecp = nn0 + er
# theta (arc_distance sector angle)
thetap = -np.degrees(ade / er) + 90.0
# zero degree elevation with standard refraction
bh0 = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=4 / 3)
# plot (ecp is earth surface normal null)
(bes,) = paax.plot(thetap, ecp, "-k", linewidth=3, label="Earth Surface NN")
(bc0r,) = paax.plot(thetap, ecp + bh0, "-g", label="0 deg Refraction")
(bc0n,) = paax.plot(thetap, ecp + bhe, "-r", label="0 deg line of sight")
thetap_k = -np.degrees(arc_distance / er) + 90.0
for i in range(len(beam_height) - 1):
(bc,) = paax.plot(
thetap_k[i],
ecp + beam_height[i, :],
linestyle="-",
linewidth=3,
label=f"Beam Center {theta[i]}deg",
)
(b3db,) = paax.plot(
thetap_k[i],
ecp + beam_height[i, :] + beamradius,
linestyle=":",
color="grey",
)
paax.plot(
thetap_k[i],
ecp + beam_height[i, :] - beamradius,
linestyle=":",
color="grey",
)
(bc,) = paax.plot(
thetap_k[i],
ecp + beam_height[i + 1, :],
linestyle="-",
linewidth=3,
label=f"Beam Center {theta[i + 1]}deg",
)
(b3db,) = paax.plot(
thetap_k[i],
ecp + beam_height[i + 1, :] + beamradius,
linestyle=":",
color="grey",
label="+3 dB Beam width",
)
paax.plot(
thetap_k[i],
ecp + beam_height[i + 1, :] - beamradius,
linestyle=":",
color="grey",
label="-3 dB Beam width",
)
# orography
# paax.fill_between(thetap, ecp,
# ecp + polarvalues[angle, :],
# color='0.75')
# shape axes
cgax.set_xlim(0, np.max(ade))
cgax.set_ylim([ecp.min() - 1000, ecp.max() + 1000])
caax.grid(True, axis="x")
cgax.grid(True, axis="y")
cgax.axis["top"].toggle(all=False)
caax.yaxis.set_major_locator(
mpl.ticker.MaxNLocator(steps=[1, 2, 4, 5, 10], nbins=20, prune="both")
)
caax.xaxis.set_major_locator(mpl.ticker.MaxNLocator())
caax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(height_formatter))
caax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(range_formatter))
caax.set_xlabel("Range (km)")
caax.set_ylabel("Altitude (km)")
# Shrink current axis by 20%
box = cgax.get_position()
cgax.set_position([box.x0, box.y0, box.width, box.height * 0.75])
paax.legend(bbox_to_anchor=(0.5, 1.5), loc="upper center", ncol=2)
[docs]
def noise_correction(ds, dbz, rho, noise_level, gas_att=0.033):
"""Calculate SNR, apply correction to RHOHV.
Dataset needs range-coordinate in m.
Parameter
---------
ds : xarray.Dataset
Dataset with radar moments
dbz : str
name of reflectivity moment
rho : str
name of crosscorrelation coefficient moment
noise_level : float
noise level for calculation
Keyword Arguments
-----------------
gas_att : float
gaseous attenuation factor
Return
------
ds : xarray.Dataset
Dataset with added corrected moments
"""
# noise calculations
snrh = (
ds[dbz]
- 20 * np.log10(ds.range * 0.001)
- noise_level
- gas_att * ds.range / 1000
)
# set nan to 0
snrh = snrh.where(snrh >= 0).fillna(0)
attrs = xd.model.sweep_vars_mapping["SNRH"]
attrs.pop("gamic", None)
snrh = snrh.assign_attrs(attrs)
rhohv = ds[rho] * np.sqrt(1.0 + 1.0 / 10.0 ** (snrh * 0.1))
rhohv = rhohv.assign_attrs(ds[rho].attrs)
ds = ds.assign({"SNRH": snrh, rho + "_NC": rhohv})
ds = ds.assign_coords({"NOISE_LEVEL": noise_level})
return ds
def noise_correction2(dbz, rho, noise_level, gas_att=0.033):
"""Calculate SNR, apply to RHOHV.
DataArray needs range-coordinate in m.
Parameter
---------
dbz : xarray.DataArray
reflectivity DataArray
rho : xarray.DataArray
crosscorrelation coefficient DataArray
noise_level : float
noise level for calculation
Keyword Arguments
-----------------
gas_att : float
gaseous attenuation factor
Return
------
ds : xarray.Dataset
Dataset with added corrected moments
"""
# noise calculations
rng_km = dbz.range * 0.001
snrh = dbz - 20 * np.log10(rng_km) - noise_level - gas_att * rng_km
# set nan to 0
snrh = snrh.where(snrh >= 0).fillna(0)
attrs = xd.model.sweep_vars_mapping["SNRH"]
attrs.pop("gamic", None)
snrh = snrh.assign_attrs(attrs)
snrh.name = "SNR"
rho_nc = rho * np.sqrt(1.0 + 1.0 / 10.0 ** (snrh * 0.1))
rho_nc = rho_nc.assign_attrs(rho.attrs)
rho_nc.name = "RHO"
return snrh, rho_nc
def noise_rho_hist(
dbz, rho, noise=(-40, -20, 2), rho_bins=(0.9, 1.1, 0.005), snr_bins=(5.0, 30.0, 0.1)
):
"""Calculate SNR, apply to RHOHV.
DataArray needs range-coordinate in m.
Parameter
---------
dbz : xarray.DataArray
reflectivity DataArray
rho : xarray.DataArray
crosscorrelation coefficient DataArray
Keyword Arguments/RADAR_NG/
-----------------
noise : tuple
tuple to create noise-array
rho_bins : tuple
tuple to create rho-array
snr_bins : tuple
tuple to create snr-array
Return
------
hists : xarray.DataArray
DataArray with Noise Calculations
"""
noise = np.arange(*noise)
rho_bins = np.arange(*rho_bins)
snr_bins = np.arange(*snr_bins)
corr = [noise_correction2(dbz, rho, n) for n in noise]
hist = [
dask.delayed(histogram)(
rho0, snr0, bins=[rho_bins, snr_bins], block_size=rho.time.size
)
for snr0, rho0 in corr
]
with ProgressBar():
hist = dask.compute(*hist)
hists = xr.concat(hist, "noise").to_dataset()
hists = hists.assign_coords({"noise": (["noise"], noise)})
std = hists.idxmax("RHO_bin").std(dim="SNR_bin", skipna=True).histogram_RHO_SNR
hists = hists.assign_coords({"stdev": std})
return hists
# def noise_rho_hist(
# dbz, rho, noise=(-40, -10, 2), rho_bins=(0.9, 1.1, 0.005), snr_bins=(5.0, 30.0, 0.1)
# ):
# """Calculate SNR, apply to RHOHV.
#
# DataArray needs range-coordinate in m.
#
# Parameter
# ---------
# dbz : xarray.DataArray
# reflectivity DataArray
# rho : xarray.DataArray
# crosscorrelation coefficient DataArray
#
# Keyword Arguments
# -----------------
# noise : tuple
# tuple to create noise-array
# rho_bins : tuple
# tuple to create rho-array
# snr_bins : tuple
# tuple to create snr-array
#
# Return
# ------
# hists : xarray.DataArray
# DataArray with Noise Calculations
#
# """
# noise = np.arange(*noise)
# rho_bins = np.arange(*rho_bins)
# snr_bins = np.arange(*snr_bins)
# corr = [noise_correction2(dbz, rho, n) for n in noise]
# hist = [
# dask.delayed(histogram)(
# rho0, snr0, bins=[rho_bins, snr_bins], block_size=rho.time.size
# )
# for snr0, rho0 in corr
# ]
# with ProgressBar():
# hist = dask.compute(*hist)
# hists = xr.concat(hist, "noise").to_dataset()
# hists = hists.assign_coords({"noise": (["noise"], noise)})
# std = hists.idxmax("RHO_bin").std(dim="SNR_bin", skipna=True).histogram_RHO_SNR
# hists = hists.assign_coords({"stdev": std})
# return hists
[docs]
def noise_hist(dbz, rho, noise=(-40, -10, 1)):
hist = noise_rho_hist(dbz, rho, noise=noise)
nl = hist.stdev.idxmin()
noise = (nl - 2.0, nl + 2.1, 0.5)
hist = noise_rho_hist(dbz, rho, noise=noise)
nl = hist.stdev.idxmin()
noise = (nl - 0.5, nl + 0.51, 0.1)
hist = noise_rho_hist(dbz, rho, noise=noise)
return hist
[docs]
def noise_level(hist):
nl = hist.stdev.idxmin()
return nl
[docs]
def estimate_cmap(ds, dbzsum, upper_frequency=0.05, lower_frequency=0.05, name="CMAP"):
"""Retrieve Clutter and Beam blockage.
Dataset needs range-coordinate in m.
Parameter
---------
ds : xarray.Dataset
Dataset with radar moments
dbzsum : xarray.DataArray
DataArray with reflectivity sum
Return
------
ds : xarray.Dataset
Dataset with added clmap moment
"""
# clutter calculation
clmap = wrl.classify.histo_cut(
dbzsum, upper_frequency=upper_frequency, lower_frequency=lower_frequency
).compute()
attrs = xd.model.sweep_vars_mapping["CMAP"]
clmap = clmap.assign_attrs(attrs)
ds = ds.assign({name: clmap})
return ds
[docs]
def get_range_step(ds):
return np.diff(ds.range)[0]
[docs]
def get_bins_from_range(ds, rng):
range_step = get_range_step(ds)
return int(rng / range_step)
def _interpolate_dem(demc, demd, polc):
polarvalues = wrl.ipol.map_coordinates(demc, demd, polc, order=3, prefilter=False)
return polarvalues
[docs]
def interpolate_dem(swp, dem):
# swp = swp.copy()
crs_wkt = dem.coords.get("crs_wkt") or dem.coords.get("spatial_ref", False)
if crs_wkt is not False:
proj_crs = pyproj.CRS.from_cf(crs_wkt.attrs)
else:
raise ValueError("missing spatial ref")
proj_dem = osr.SpatialReference()
proj_dem.ImportFromWkt(proj_crs.to_wkt())
coords = swp.wrl.georef.reproject(trg_crs=proj_dem).reset_coords()
swp = swp.assign_coords(xd=coords.x, yd=coords.y)
polc = xr.concat([coords.x, coords.y], "xy").transpose(..., "xy")
dsign = np.sign(dem.y.diff("y")[0]).astype(int).values
dem = dem.sel(
x=slice(coords.x.min().values, coords.x.max().values),
y=slice(coords.y.min().values, coords.y.max().values, dsign),
)
demc = xr.concat(xr.broadcast(dem.x, dem.y), dim="xy").transpose("y", "x", "xy")
pol = xr.apply_ufunc(
_interpolate_dem,
demc,
dem,
polc,
input_core_dims=[["y", "x", "xy"], ["y", "x"], ["azimuth", "range", "xy"]],
output_core_dims=[["azimuth", "range"]],
)
swp = swp.assign(DEM=pol)
return swp
[docs]
def get_beamblockage(swp, bw=1.0):
bw = 1.0
beamradius = wrl.util.half_power_radius(swp.range.values, bw)
PBB = wrl.qual.beam_block_frac(swp.DEM.values, swp.z.values, beamradius)
PBB = np.ma.masked_invalid(PBB)
CBB = wrl.qual.cum_beam_block_frac(PBB)
swp = swp.assign(PBB=(["azimuth", "range"], PBB), CBB=(["azimuth", "range"], CBB))
return swp
[docs]
def get_qbbf(swp):
"""
Quality Index for PBB
Crisologo 2018
"""
Qbbf = xr.zeros_like(swp.CBB) * np.nan
Qbbf = xr.where(swp.CBB < 0.1, 1, Qbbf)
Qbbf = xr.where(swp.CBB > 0.5, 0, Qbbf)
Qbbf = xr.where(
(swp.CBB <= 0.5) & (swp.CBB >= 0.1), 1 - (swp.CBB - 0.1) / 0.4, Qbbf
)
Qbbf.name = "QBBF"
swp = swp.assign(QBBF=Qbbf)
return swp
[docs]
def save_dwd_station(station_id, start_time, stop_time, path=".", reload=False):
print(f"processing station: {int(station_id):05d}")
outfile = os.path.join(
path, f"{int(station_id):05d}_{start_time:%Y%m%d%H%M}_{stop_time:%Y%m%d%H%M}.nc"
)
# check if ncfile already exists
if os.path.isfile(outfile) and reload is False:
print(f"--- station {int(station_id):05d} already exists")
return
# check if dwd file already exists
fname = f"1minutenwerte_nieder_{station_id:05d}_{start_time:%Y%m}01*_hist.zip"
if not os.path.isfile(fname):
parser = DWDHTMLParser()
opendata_url1 = f"https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/1_minute/precipitation/historical/{start_time.year}/"
http = urllib3.PoolManager(cert_reqs="CERT_REQUIRED", ca_certs=certifi.where())
response = http.request("GET", opendata_url1).data.decode("utf-8")
parser.links = []
parser.feed(response)
filelist = parser.links[1:]
matching = fnmatch.filter(filelist, fname)
for f in matching:
with (
http.request(
"GET", os.path.join(opendata_url1, f), preload_content=False
) as r,
open(f, "wb") as out,
):
shutil.copyfileobj(r, out)
if not matching:
return
columns_to_read = ["STATIONS_ID", "MESS_DATUM_BEGINN", "RS_01"]
df = pd.read_csv(
matching[0], sep=";", parse_dates=["MESS_DATUM_BEGINN"], usecols=columns_to_read
)
df.columns = ["station_id", "time", "precipitation_height"]
rain = df.to_xarray()
rain.precipitation_height.attrs = dict(short_name="precipitation_height", unit="mm")
rain = rain.set_coords(["time", "station_id"])
rain = rain.swap_dims(index="time")
rain = rain.drop("index")
rain["station_id"] = rain["station_id"][0]
# select time range
rain = rain.sel(time=slice(start_time, stop_time))
# do not write, if no values available
if rain.sizes["time"] == 0:
return
# extract time strings
start_time_str = rain.time.dt.strftime("%Y%m%d%H%M")[0].values
stop_time_str = rain.time.dt.strftime("%Y%m%d%H%M")[-1].values
outfile = f"{int(station_id):05d}_{start_time_str}_{stop_time_str}.nc"
print("--- save to ", outfile)
rain.to_netcdf(outfile)
[docs]
def process_regression(uh, rlog):
from scipy import stats
mask = ~np.isnan(rlog) & ~np.isnan(uh)
uh = uh.where(mask, drop=True)
rl = rlog.where(mask, drop=True)
a = np.nan
b = np.nan
slope = np.nan
intercept = np.nan
if len(uh) > 0:
slope, intercept, r_value, p_value, std_err = stats.linregress(uh, rl)
print(f"m:{slope}, n:{intercept}")
b = 1.0 / (10.0 * slope)
a = 10 ** (-b * intercept)
print(f"a:{a}, b:{b}")
return a, b, slope, intercept
[docs]
def read_dwd_stations(fname):
dwd_stations = pd.read_csv(
fname, encoding="iso-8859-1", nrows=1, delimiter=" ", engine="python"
)
names = list(dwd_stations.to_xarray().data_vars)
width = [6, 9, 9, 15, 12, 10, 39]
dwd_stations = pd.read_fwf(
fname, width=width, encoding="iso-8859-1", skiprows=2, header=None
)
dwd_stations.columns = names
dwd = dwd_stations.to_xarray().rename(
Stations_id="station_id", Stationsname="station_name"
)
return dwd
[docs]
def filter_stations(dwd, swp, distance):
x = dwd.geoLaenge
y = dwd.geoBreite
points = np.squeeze(np.dstack((x, y)))
rng = swp.sel(range=distance, method="nearest")
vertices = np.squeeze(np.dstack((rng.xd, rng.yd)))
pth = Path(vertices)
dwdx = dwd.sel(index=pth.contains_points(points))
dwdx = dwdx.set_coords("station_id")
dwdx = dwdx.swap_dims(index="station_id")
return dwdx
[docs]
def filter_radar_bins(swp, dwdx, n_nearest):
ds = swp.stack(bins=["azimuth", "range"])
ds = ds.assign_coords(xyd=xr.concat([ds.xd, ds.yd], "xy").transpose(..., "xy"))
tree = KDTree(ds.xyd)
dist, ix = tree.query(
list(zip(dwdx.geoLaenge, dwdx.geoBreite, strict=True)), k=n_nearest
)
dwdx = dwdx.assign(
bins=(["station_id", "bin"], ix), dist=(["station_id", "bin"], dist)
)
ds_stations = ds.isel(bins=dwdx.bins)
return ds_stations, dwdx
def workflow(sweep_ts, scan, path, cmap2=False):
dim0 = sweep_ts.wrl.util.dim0()
# DEM Processing
# Get DEM-data
if sweep_ts.sweep_mode.values == "azimuth_surveillance":
rasterfile = wrl.util.get_wradlib_data_file("geo/bonn_dem.tif")
dem = xr.open_dataset(rasterfile, engine="rasterio")
# DEM-data on radar grid
sweep_ts = sweep_ts.pipe(interpolate_dem, dem.band_data[0])
# Beamblockage
sweep_ts = sweep_ts.pipe(get_beamblockage)
# QI(BBF)
sweep_ts = sweep_ts.pipe(get_qbbf)
# Clutter and Noise
if cmap2:
start_time = sweep_ts.time.dt.date.min().load().item()
start_time_day = dt.datetime(start_time.year, start_time.month, start_time.day)
delta = dt.timedelta(days=1)
sweep_ts_day = load_timeseries_sweep(start_time_day, delta, scan)
DBTH = sweep_ts_day.DBTH.sum("time2")
sweep_ts = sweep_ts.pipe(
estimate_cmap,
DBTH,
upper_frequency=0.05,
lower_frequency=0.01,
name="CMAP2",
)
noise = (-40, -10, 1)
hist = noise_hist(sweep_ts.DBZH, sweep_ts.RHOHV, noise=noise)
nl = noise_level(hist)
sweep_ts = sweep_ts.pipe(noise_correction, "DBZH", "RHOHV", nl)
# System PHIDP
if "CMAP" not in sweep_ts:
if "time2" in sweep_ts:
CMAP = xr.zeros_like(sweep_ts.DBTH.isel(time2=0), "uint16")
else:
CMAP = xr.zeros_like(sweep_ts.DBTH, "uint16")
# sweep_ts = sweep_ts.assign(CMAP=xr.zeros_like(sweep_ts.DBTH.isel(time2=0), "uint16"))
else:
CMAP = sweep_ts.CMAP.astype("uint16")
if "CMAP2" not in sweep_ts:
if "time2" in sweep_ts:
CMAP2 = xr.zeros_like(sweep_ts.DBTH.isel(time2=0), "uint8")
else:
CMAP2 = xr.zeros_like(sweep_ts.DBTH, "uint8")
# sweep_ts = sweep_ts.assign(CMAP=xr.zeros_like(sweep_ts.DBTH.isel(time2=0), "uint16"))
else:
CMAP2 = sweep_ts.CMAP2
bit = 0
cmp = (CMAP >> bit) & 0x1
# mask phase
no_clutter_bb = (cmp == 0) & (CMAP2 == 0)
phimask = sweep_ts.PHIDP.where(no_clutter_bb & (sweep_ts.RHOHV_NC >= 0.9))
# create binary array
phib = xr.where(np.isnan(phimask), 0, 1)
# Phase V1
rng = 2000.0
N = sweep_ts.pipe(get_bins_from_range, rng)
# take N range bins and calculate sum, for finding N consecutive bins
phib_sum = phib.rolling(range=N, center=True).sum(skipna=True)
# get maximum along range
smax = phib_sum.max(dim="range", skipna=True)
# if smax == N, we've identified the first consecutive range of valid phase values
# derive start_range and stop_range
rstep = get_range_step(phib_sum)
start_range = phib_sum.idxmax(dim="range").where(smax == N) - smax // 2 * rstep
stop_range = start_range + N * rstep
phi_1 = phimask.where(
(phimask.range >= start_range) & (phimask.range <= stop_range)
)
phi_1m = phi_1.median("range", skipna=True)
phi_1m_sort = xr.apply_ufunc(
np.sort,
phi_1m,
input_core_dims=[[dim0]],
output_core_dims=[[dim0]],
dask="parallelized",
)
sysphi_1 = phi_1m_sort.isel({dim0: slice(0, 30)}).median(dim0, skipna=True)
# Phase V2
# get the maximum along the ray
# derive start_range and stop_range
start_range = phib_sum.idxmax(dim="range") - N // 2 * rstep
stop_range = start_range + N * rstep
phi_2 = phimask.where(
(phimask.range >= start_range) & (phimask.range <= stop_range)
)
# phi_2.plot(x="azimuth")
phi_2m = phi_2.median("range", skipna=True)
phi_2m_sort = xr.apply_ufunc(
np.sort,
phi_2m,
input_core_dims=[[dim0]],
output_core_dims=[[dim0]],
dask="parallelized",
)
sysphi_2 = phi_2m_sort.isel({dim0: slice(0, 30)}).median(dim0, skipna=True)
# Phase V3
# cumulative sum, for finding first N valid bins
phib_cumsum = phib.cumsum("range", skipna=True)
# select only values until N is reached
phi_3 = phimask.where(phib_cumsum <= N)
# phi_3.plot(x="azimuth", vmin=75, vmax=95)
phi_3m = phi_3.median("range", skipna=True)
phi_3m_sort = xr.apply_ufunc(
np.sort,
phi_3m,
input_core_dims=[[dim0]],
output_core_dims=[[dim0]],
dask="parallelized",
)
sysphi_3 = phi_3m_sort.isel({dim0: slice(0, 30)}).median(dim0, skipna=True)
# Phase V4
phase_res = 0.1
bins = np.arange(-180, 180, phase_res)
phi_hist = histogram(sweep_ts.PHIDP, dim=("range",), bins=[bins])
phi_hist_smooth = phi_hist.rolling(PHIDP_bin=5, center=True).mean()
pmax = phi_hist_smooth.idxmax("PHIDP_bin").compute()
phi_hist_smooth_cut = phi_hist_smooth.where(
(phi_hist.PHIDP_bin > (pmax - 5)) & (phi_hist.PHIDP_bin < (pmax + 5)), drop=True
)
phi_hist_smooth_cut = phi_hist_smooth_cut / phi_hist_smooth_cut.max()
phi_hist_smooth_cut = phi_hist_smooth_cut.chunk(
PHIDP_bin=phi_hist_smooth_cut.sizes["PHIDP_bin"]
)
sysphi_4 = (
phi_hist_smooth_cut.differentiate("PHIDP_bin")
.rolling(PHIDP_bin=11, center=True)
.mean()
.idxmax("PHIDP_bin")
)
with xr.set_options(keep_attrs=True):
sweep_ts = sweep_ts.assign({"SYSPHI1": sysphi_1})
sweep_ts = sweep_ts.assign({"PHIDP_OC1": sweep_ts.PHIDP - sysphi_1})
sweep_ts = sweep_ts.assign({"SYSPHI2": sysphi_2})
sweep_ts = sweep_ts.assign({"PHIDP_OC2": sweep_ts.PHIDP - sysphi_2})
sweep_ts = sweep_ts.assign({"SYSPHI3": sysphi_3})
sweep_ts = sweep_ts.assign({"PHIDP_OC3": sweep_ts.PHIDP - sysphi_3})
sweep_ts = sweep_ts.assign({"SYSPHI4": sysphi_4})
sweep_ts = sweep_ts.assign({"PHIDP_OC4": sweep_ts.PHIDP - sysphi_4})
# Phase Processing
# CMAP = sweep_ts.CMAP.astype("uint16")
noise = ((CMAP >> 0) & 0x1).astype(bool)
clutter = ((CMAP >> 2) & 0x1).astype(bool)
masked = noise | clutter
valid = ~masked
# CMAP2 = sweep_ts.CMAP2
clutter2 = CMAP2 == 1
bb = CMAP == 2
masked2 = clutter2 | bb
valid2 = ~masked2
val = valid & valid2
sweep_ts_masked = (
sweep_ts.where(val).where(sweep_ts.RHOHV_NC >= 0.95).where(sweep_ts.DBZH >= 0)
)
sweep_ts = sweep_ts.assign(
RHOHV_SMOOTH=sweep_ts_masked.RHOHV_NC.rolling(
range=11, min_periods=2, center=True
).median(skipna=True)
)
medwin = 41
PHIDP_OC = "PHIDP_OC1"
# median filtering 2d
phimed = sweep_ts_masked[PHIDP_OC].copy()
with xr.set_options(keep_attrs=True):
phimed = phimed.pipe(filter_data, medwin)
# gaussian convolution 1d - smoothing
# play with sigma and kwidth
kwidth = 41
sigma = 10
gkern = gauss_kernel(kwidth, sigma)
with xr.set_options(keep_attrs=True):
phiclean = phimed.pipe(smooth_data, gkern)
window = 41
window2 = 3
dim0 = sweep_ts_masked.wrl.util.dim0()
with xr.set_options(keep_attrs=True):
phi_median = (
sweep_ts_masked[PHIDP_OC]
.copy()
.pipe(
xr_rolling,
window,
window2=window2,
method="median",
skipna=True,
min_periods=3,
)
.transpose(..., dim0, "range")
)
# gaussian convolution 1d - smoothing
# play with sigma and kwidth
kwidth = 41
sigma = 10
gkern = gauss_kernel(kwidth, sigma)
phi_clean = phi_median.pipe(smooth_data, gkern)
sweep_ts = sweep_ts.assign(
{
"PHIDP_CORR1": phimed,
"PHIDP_CORR2": phi_median,
"PHIDP_CORR3": phiclean,
"PHIDP_CORR4": phi_clean,
}
)
PHIDP_CORR = "PHIDP_CORR1"
L = 11
kdp_lanczos = sweep_ts[PHIDP_CORR].wrl.dp.kdp_from_phidp(winlen=L, min_periods=2)
phidp_v, kdp_vulpiani = sweep_ts[PHIDP_CORR].wrl.dp.phidp_kdp_vulpiani(winlen=L)
sweep_ts = sweep_ts.assign({"KDP_LANC": kdp_lanczos, "KDP_VULP": kdp_vulpiani})
# Attenuation Correction
uh = "DBTH"
zh = "DBZH"
zdr = "ZDR"
phi_list = [f"PHIDP_CORR{i}" for i in range(1, 5)]
alphax = 0.28
betax = 0.05
with xr.set_options(keep_attrs=True):
for phi in phi_list:
uhcorr = sweep_ts[uh] + alphax * sweep_ts[phi]
uhcorr.name = f"{sweep_ts[uh].name}_{sweep_ts[phi].name.split('_')[-1]}"
zhcorr = sweep_ts[zh] + alphax * sweep_ts[phi]
zhcorr.name = f"{sweep_ts[zh].name}_{sweep_ts[phi].name.split('_')[-1]}"
zdrcorr = sweep_ts[zdr] + betax * sweep_ts[phi]
zdrcorr.name = f"{sweep_ts[zdr].name}_{sweep_ts[phi].name.split('_')[-1]}"
sweep_ts = sweep_ts.assign(
{uhcorr.name: uhcorr, zhcorr.name: zhcorr, zdrcorr.name: zdrcorr}
)
# ZPHI - Method
# CMAP = sweep_ts.CMAP.astype("uint16")
noise = ((CMAP >> 0) & 0x1).astype(bool)
clutter = ((CMAP >> 2) & 0x1).astype(bool)
masked = noise | clutter
valid = ~masked
phi = sweep_ts.PHIDP.where((CMAP2 == 0) & valid & (sweep_ts.RHOHV_NC > 0.9))
rng = 2000.0
zh = sweep_ts.DBZH.where((CMAP2 == 0) & valid & (sweep_ts.RHOHV_NC > 0.9))
ah = zphi(phi, zh, rng=rng, alphax=np.array([0.28])).squeeze()
kdp_ah = ah.fillna(0) / ah.alpha
kdp_ah.attrs = sweep_ts.KDP.attrs
phical = phidp_from_kdp(kdp_ah)
phical = phical.rename("PHICAL")
phical.attrs = sweep_ts.PHIDP.attrs
sweep_ts = sweep_ts.assign(
AH=ah,
KDP_AH=kdp_ah,
PHIDP_CORR5=phical,
)
uh = "DBTH"
zh = "DBZH"
zdr = "ZDR"
phi = "PHIDP_CORR5"
with xr.set_options(keep_attrs=True):
uhcorr = sweep_ts[uh] + sweep_ts.alpha * sweep_ts[phi]
uhcorr.name = f"{sweep_ts[uh].name}_{sweep_ts[phi].name.split('_')[-1]}"
zhcorr = sweep_ts[zh] + sweep_ts.alpha * sweep_ts[phi]
zhcorr.name = f"{sweep_ts[zh].name}_{sweep_ts[phi].name.split('_')[-1]}"
zdrcorr = sweep_ts[zdr] + betax * sweep_ts[phi]
zdrcorr.name = f"{sweep_ts[zdr].name}_{sweep_ts[phi].name.split('_')[-1]}"
sweep_ts = sweep_ts.assign(
{uhcorr.name: uhcorr, zhcorr.name: zhcorr, zdrcorr.name: zdrcorr}
)
# Save data
dtime_min = sweep_ts.time.load().min().dt.strftime("%Y%m%dT%H%M").values
dtime_max = sweep_ts.time.load().max().dt.strftime("%Y%m%dT%H%M").values
fname = os.path.join(path, f"tmp_{dtime_min}_{dtime_max}_workflow.nc")
print(fname)
sweep_ts.to_netcdf(fname)
return sweep_ts
[docs]
def get_temperature_profile(rs_data, hmax=30000.0):
stemp = rs_data["TEMP"]
sheight = rs_data["HGHT"]
# remove nans
idx = np.isfinite(stemp)
stemp = stemp[idx]
sheight = sheight[idx]
stemp_da = xr.DataArray(
data=stemp,
dims=["height"],
coords=dict(height=(["height"], sheight)),
attrs=dict(description="Temperature", units="degC", long_name="Temperature"),
)
ht = np.arange(0.0, hmax)
itemp_da = stemp_da.interp({"height": ht})
itemp_da = itemp_da.bfill(dim="height")
return itemp_da
[docs]
def merge_radar_profile(rds, cds):
cds = cds.interp({"height": rds.z}, method="linear")
rds = rds.assign({"TEMP": cds})
return rds
