xradar is in an early stage of development, please report any issues here!

Create a Plan Position Indicator (PPI) Plot#

A Plan Position Indicator (PPI) plot is a common plot requested by radar scientists. Let’s show how to create this plot using xradar

Imports#

[1]:
import xarray as xr
import xradar as xd
import cmweather
from open_radar_data import DATASETS
[2]:
import matplotlib.pyplot as plt
import pyproj
import cartopy

Read in some data#

Fetching CfRadial1 radar data file from open-radar-data repository.

[3]:
filename = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc")

Read the data using the cfradial1 engine

[4]:
radar = xd.io.open_cfradial1_datatree(filename, first_dim="auto")
display(radar)
<xarray.DatasetView> Size: 477B
Dimensions:              (sweep: 9)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int32 4B ...
    platform_type        |S32 32B ...
    primary_axis         |S32 32B ...
    status_str           |S1 1B ...
    instrument_type      |S32 32B ...
    time_coverage_start  |S32 32B ...
    time_coverage_end    |S32 32B ...
    latitude             float64 8B ...
    longitude            float64 8B ...
    altitude             float64 8B ...
    sweep_group_name     (sweep) <U7 252B 'sweep_0' 'sweep_1' ... 'sweep_8'
    sweep_fixed_angle    (sweep) float32 36B ...
Attributes: (12/13)
    Conventions:         CF/Radial instrument_parameters radar_parameters rad...
    version:             1.2
    title:               TIMREX
    institution:
    references:
    source:
    ...                  ...
    comment:
    instrument_name:     SPOLRVP8
    site_name:
    scan_name:
    scan_id:             0
    platform_is_mobile:  false

Add georeferencing#

We can use the georeference function, or the accessor to add our georeference information!

Georeference Accessor#

If you prefer the accessor (.xradar.georefence()), this is how you would add georeference information to your radar object.

[5]:
radar = radar.xradar.georeference()
display(radar)
<xarray.DatasetView> Size: 477B
Dimensions:              (sweep: 9)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int32 4B ...
    platform_type        |S32 32B ...
    primary_axis         |S32 32B ...
    status_str           |S1 1B ...
    instrument_type      |S32 32B ...
    time_coverage_start  |S32 32B ...
    time_coverage_end    |S32 32B ...
    latitude             float64 8B ...
    longitude            float64 8B ...
    altitude             float64 8B ...
    sweep_group_name     (sweep) <U7 252B 'sweep_0' 'sweep_1' ... 'sweep_8'
    sweep_fixed_angle    (sweep) float32 36B ...
Attributes: (12/13)
    Conventions:         CF/Radial instrument_parameters radar_parameters rad...
    version:             1.2
    title:               TIMREX
    institution:
    references:
    source:
    ...                  ...
    comment:
    instrument_name:     SPOLRVP8
    site_name:
    scan_name:
    scan_id:             0
    platform_is_mobile:  false

Please observe, that the additional coordinates x, y, z have been added to the dataset. This will also add spatial_ref CRS information on the used Azimuthal Equidistant Projection.

[6]:
radar["sweep_0"]
[6]:
<xarray.DatasetView> Size: 10MB
Dimensions:                    (azimuth: 483, range: 996)
Coordinates:
    time                       (azimuth) datetime64[ns] 4kB 2008-06-04T00:15:...
  * range                      (range) float32 4kB 150.0 300.0 ... 1.494e+05
    elevation                  (azimuth) float32 2kB 0.5164 0.5219 ... 0.5219
  * azimuth                    (azimuth) float32 2kB 0.0 0.75 ... 358.5 359.2
    latitude                   float64 8B 22.53
    longitude                  float64 8B 120.4
    altitude                   float64 8B 45.0
    crs_wkt                    int64 8B 0
    x                          (azimuth, range) float32 2MB 0.0 ... -1.955e+03
    y                          (azimuth, range) float32 2MB 150.0 ... 1.493e+05
    z                          (azimuth, range) float32 2MB 46.0 ... 2.718e+03
Data variables: (12/18)
    sweep_number               int32 4B ...
    sweep_mode                 <U20 80B 'azimuth_surveillance'
    prt_mode                   |S32 32B ...
    follow_mode                |S32 32B ...
    sweep_fixed_angle          float32 4B ...
    pulse_width                (azimuth) timedelta64[ns] 4kB ...
    ...                         ...
    r_calib_index              (azimuth) int8 483B ...
    measured_transmit_power_h  (azimuth) float32 2kB ...
    measured_transmit_power_v  (azimuth) float32 2kB ...
    scan_rate                  (azimuth) float32 2kB ...
    DBZ                        (azimuth, range) float32 2MB ...
    VR                         (azimuth, range) float32 2MB ...

Use the Function#

We can also use the function xd.geoference.get_x_y_z_tree function if you prefer that method.

[7]:
radar = xd.georeference.get_x_y_z_tree(radar)
display(radar["sweep_0"])
<xarray.DatasetView> Size: 10MB
Dimensions:                    (azimuth: 483, range: 996)
Coordinates:
    time                       (azimuth) datetime64[ns] 4kB 2008-06-04T00:15:...
  * range                      (range) float32 4kB 150.0 300.0 ... 1.494e+05
    elevation                  (azimuth) float32 2kB 0.5164 0.5219 ... 0.5219
  * azimuth                    (azimuth) float32 2kB 0.0 0.75 ... 358.5 359.2
    latitude                   float64 8B 22.53
    longitude                  float64 8B 120.4
    altitude                   float64 8B 45.0
    crs_wkt                    int64 8B 0
    x                          (azimuth, range) float32 2MB 0.0 ... -1.955e+03
    y                          (azimuth, range) float32 2MB 150.0 ... 1.493e+05
    z                          (azimuth, range) float32 2MB 46.0 ... 2.718e+03
Data variables: (12/18)
    sweep_number               int32 4B ...
    sweep_mode                 <U20 80B 'azimuth_surveillance'
    prt_mode                   |S32 32B ...
    follow_mode                |S32 32B ...
    sweep_fixed_angle          float32 4B ...
    pulse_width                (azimuth) timedelta64[ns] 4kB ...
    ...                         ...
    r_calib_index              (azimuth) int8 483B ...
    measured_transmit_power_h  (azimuth) float32 2kB ...
    measured_transmit_power_v  (azimuth) float32 2kB ...
    scan_rate                  (azimuth) float32 2kB ...
    DBZ                        (azimuth, range) float32 2MB ...
    VR                         (azimuth, range) float32 2MB ...

Plot our Data#

Plot simple PPI#

Now, let’s create our PPI plot! We just use the newly created 2D-coordinates x and y to create a meshplot.

[8]:
radar["sweep_0"]["DBZ"].plot(x="x", y="y", cmap="ChaseSpectral")
[8]:
<matplotlib.collections.QuadMesh at 0x7f0d4e60c990>
../_images/notebooks_plot-ppi_16_1.png

Plot PPI on map with cartopy#

If you have cartopy installed, you can easily plot on maps. We first have to extract the CRS from the dataset and to wrap it in a cartopy.crs.Projection.

[9]:
proj_crs = xd.georeference.get_crs(radar["sweep_0"].ds)
cart_crs = cartopy.crs.Projection(proj_crs)

Second, we create a matplotlib GeoAxes and a nice map.

[10]:
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=cartopy.crs.PlateCarree())
radar["sweep_0"]["DBZ"].plot(
    x="x",
    y="y",
    cmap="ChaseSpectral",
    transform=cart_crs,
    cbar_kwargs=dict(pad=0.075, shrink=0.75),
)
ax.coastlines()
ax.gridlines(draw_labels=True)
/home/docs/checkouts/readthedocs.org/user_builds/xradar/conda/latest/lib/python3.11/site-packages/cartopy/mpl/geoaxes.py:1781: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
  result = super().pcolormesh(*args, **kwargs)
[10]:
<cartopy.mpl.gridliner.Gridliner at 0x7f0d4c96c450>
../_images/notebooks_plot-ppi_20_2.png