Plot moment in polar and projected reference systems#
Import section#
import datetime as dt
import matplotlib.pyplot as plt
import wradlib as wrl
Utility Functions#
import radarpolarimetry as rp
Load Single Sweep#
# choose radar specific parameters
start_time = dt.datetime(2017, 8, 10, 11, 0)
scan = "n_ppi_020deg"
swp = rp.load_single_sweep(start_time, scan)
display(swp)
<xarray.Dataset> Size: 26MB
Dimensions: (azimuth: 360, range: 750)
Coordinates: (12/15)
* azimuth (azimuth) float64 3kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
elevation (azimuth) float64 3kB 2.0 2.0 2.0 2.0 ... 2.0 2.0 2.0 2.0
time (azimuth) datetime64[ns] 3kB ...
* range (range) float32 3kB 100.0 300.0 ... 1.497e+05 1.499e+05
x (azimuth, range) float64 2MB 0.8721 2.616 ... -1.306e+03
y (azimuth, range) float64 2MB 99.93 299.8 ... 1.497e+05
... ...
bins (azimuth, range) float32 1MB 100.0 300.0 ... 1.499e+05
sweep_mode <U20 80B 'azimuth_surveillance'
longitude float64 8B 7.072
latitude float64 8B 50.73
altitude float64 8B 99.5
crs_wkt int64 8B 0
Data variables: (12/18)
DBZH (azimuth, range) float32 1MB ...
DBZV (azimuth, range) float32 1MB ...
KDP (azimuth, range) float32 1MB ...
RHOHV (azimuth, range) float32 1MB ...
CMAP (azimuth, range) float32 1MB ...
DBTH (azimuth, range) float32 1MB ...
... ...
PHIDP (azimuth, range) float32 1MB ...
sweep_number int64 8B ...
prt_mode <U7 28B ...
follow_mode <U7 28B ...
sweep_fixed_angle float64 8B ...
nyquist_velocity float64 8B ...
Attributes:
source: gamic
ant_gain_h: 44.0
ant_gain_v: 44.0
noise_power_h: -0.52988
noise_power_v: -0.339564
rx_loss_h: 0.73
rx_loss_v: 0.73Exercise 1: Simple Plot#
Plot DBZH with minimal effort.
PPI plot#
Using xarray wrl.vis-accessor PPI plots can be created. We use discrete colormaps/colorbars.
Further material on colormaps/colorbars is available for self study in Colormaps_Colorbars notebook.
# create discrete cmap
zh = rp.visdict14["ZH"]
cmap = rp.get_discrete_cmap(zh["ticks"], "HomeyerRainbow")
fig = plt.figure(figsize=(8, 8))
ax1 = fig.add_subplot()
pm = swp.DBTH.wrl.vis.plot(ax=ax1, levels=zh["ticks"], cmap=cmap)
wrl.vis.plot_ppi_crosshair(
(0, 0, 0), ranges=[25e3, 50e3, 75e3, 100e3, 125e3, 150e3], ax=ax1
)
txt = ax1.set_title("Simple PPI Plot")
Define plot_moment function#
We will use the simple plot_moment function to plot moments from now on.
fig, axes = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(12, 8))
rp.plot_moment(swp.DBTH, ax=axes.flat[0])
rp.plot_moment(swp.DBZH, ax=axes.flat[1])
plt.tight_layout()
PPI plot in curvelinear grid#
This is just an example to show wradlib plotting capabilities.
pm = swp.DBZH.wrl.vis.plot(
ax=111, crs="cg", cbar_kwargs=dict(pad=0.1), levels=zh["ticks"], cmap=cmap
)
fig = plt.figure(figsize=(8, 8))
pm = swp.DBZH.wrl.vis.plot(
ax=111,
fig=fig,
crs="cg",
cbar_kwargs=dict(pad=0.075, shrink=0.8),
levels=zh["ticks"],
cmap=cmap,
)
cgax = plt.gca()
caax = cgax.parasites[0]
paax = cgax.parasites[1]
plt.text(1.0, 1.05, "azimuth", transform=caax.transAxes, va="bottom", ha="right")
t = cgax.set_title("Simple Plot in cg grid", y=1.05)
Plot in projected coordinate systems#
Uses cartopy.
import cartopy.crs as ccrs
mercator = ccrs.Mercator(central_longitude=swp.longitude.values)
utm32 = ccrs.UTM(32)
wgs84 = ccrs.PlateCarree(central_longitude=swp.longitude.values)
Mercator#
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=mercator)
pm = swp.DBZH.wrl.vis.plot(
ax=ax, add_labels=True, cbar_kwargs=dict(pad=0.075), levels=zh["ticks"], cmap=cmap
)
t = ax.set_title("Simple Plot in Mercator", y=1.05)
ax.gridlines(draw_labels=True)
<cartopy.mpl.gridliner.Gridliner at 0x7f094d27f230>
Exercise 2: Plot on UTM32 and WGS84#
Create two plots side-by-side.
Exercise 3: Plotting just one sector#
Plot a sector of the data, azimuth between 200 and 250 and range between 10km and 30km.
Plotting Map Features#
import cartopy.feature as cfeature
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=mercator)
pm = swp.DBZH.wrl.vis.plot(
ax=ax, cbar_kwargs=dict(pad=0.1, fraction=0.0476), levels=zh["ticks"], cmap=cmap
)
ax.set_title("radar data on map features")
ax.gridlines(draw_labels=True)
# Create a feature for States/Admin 1 regions at 1:50m from Natural Earth
states_provinces = cfeature.NaturalEarthFeature(
category="cultural",
name="admin_1_states_provinces_lines",
scale="10m",
facecolor="none",
)
countries = cfeature.NaturalEarthFeature(
category="cultural", name="admin_0_countries", scale="10m", facecolor="none"
)
ax.add_feature(states_provinces, edgecolor="black", ls=":")
ax.add_feature(
countries,
edgecolor="black",
)
SOURCE = "Natural Earth"
LICENSE = "public domain"
Solution 1: Simple Plot#
Directly use xarray-machinery.
pm = swp.DBZH.plot()
txt = plt.title("Simple Range-Azimuth Plot")
Solution 2: Plot on UTM32 and WGS84#
fig = plt.figure(figsize=(16, 7))
ax1 = fig.add_subplot(121, projection=utm32)
pm = swp.DBZH.wrl.vis.plot(
ax=ax1, cbar_kwargs=dict(pad=0.1, shrink=0.75), levels=zh["ticks"], cmap=cmap
)
t = ax1.set_title("Simple Plot in UTM zone 32", y=1.05)
ax1.gridlines(draw_labels=True)
ax2 = fig.add_subplot(122, projection=wgs84)
pm = swp.DBZH.wrl.vis.plot(
ax=ax2, cbar_kwargs=dict(pad=0.1, shrink=0.75), levels=zh["ticks"], cmap=cmap
)
t = ax2.set_title("Simple Plot in WGS 84", y=1.05)
ax2.gridlines(draw_labels=True)
<cartopy.mpl.gridliner.Gridliner at 0x7f094d2d65d0>
Solution 3: Plotting just one sector#
For this purpose, we need to give the ranges and azimuths explicitly…
da_sel = swp.DBZH.sel(azimuth=slice(200, 250), range=slice(10e3, 30e3))
pm = da_sel.wrl.vis.plot(add_labels=True, levels=zh["ticks"], cmap=cmap)
txt = plt.title("Sector PPI")
