2D objective analysis of weather radar RHI scan: Fast Barnes Interpolation Example

This example uses the Fast Barnes interpolation methods as presented in Zürcher (2023) available at fast-barnes-py.

The notebook applies the fast Barnes interpolation methods to RHI scan data from the ARM-CSAPR2 radar data from the [TRACER] (https://www.arm.gov/research/campaigns/amf2021tracer) field campaign.

Overview¶

This notebook shows how to download the ARM-CSAPR2 RHI data from ARM data repository and grid radar variables on a two-dimensional Cartesian grid using Python.

Prerequisites¶

Concepts	Importance	Notes
Intro to Xarray	Necessary
Matplotlib Basics	Required	Basic plotting
Py-ART Basics	Required	IO/Visualization
Understanding of NetCDF	Helpful	Familiarity with metadata structure

Time to learn: 10 minutes

Imports¶

import pyart
import act
import glob
import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from pathlib import Path
from copy import deepcopy
from fastbarnes.interpolation import barnes
from datetime import timedelta


## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar data. Py-ART is partly supported
## by the U.S. Department of Energy Office of Science as part of
## the Atmospheric Radiation Measurement (ARM) User Facility.
##
## If you use this software to prepare a publication, please cite:
##
##     JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119

/home/runner/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Setup our Download Query¶

Before downloading our data, we need to make sure we have an ARM Data Account, and ARM Live token. Both of these can be found using this link:

ARM Live Signup

Once you sign up, you will see your token. Copy and replace that where we have arm_username and arm_password below.

arm_username = os.getenv("ARM_USERNAME")
arm_password = os.getenv("ARM_PASSWORD")

datastream = "houcsapr2cfrS2.a1"

start_date = "2022-08-07T18:39:04"
end_date = "2022-08-07T18:39:05"

Grab Data¶

We use the ARM Live API, accessible through ACT!

Download Data API

We are interested in the C-band radar, which is utilizing a cell-tracking algorithm, with the datastream

houcsapr2cfrS2.a1

One of the better cases was from June 2, 2022.

This line will download data into the a new directory with the datastream name (/houcsapr2cfrS2.a1)

aug7_csapr_files = act.discovery.download_arm_data(arm_username,
                                               arm_password,
                                               datastream,
                                               start_date,
                                               end_date)

---------------------------------------------------------------------------
HTTPError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 aug7_csapr_files = act.discovery.download_arm_data(arm_username,
      2                                                arm_password,
      3                                                datastream,
      4                                                start_date,

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/act/discovery/arm.py:122, in download_arm_data(username, token, datastream, startdate, enddate, time, output)
    119 req = Request(query_url, None, headers)
    120 # get url response, read the body of the message,
    121 # and decode from bytes type to utf-8 string
--> 122 response_body = urlopen(req).read().decode('utf-8')
    123 # if the response is an html doc, then there was an error with the user
    124 if response_body[1:14] == '!DOCTYPE html':

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:215, in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    213 else:
    214     opener = _opener
--> 215 return opener.open(url, data, timeout)

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:521, in OpenerDirector.open(self, fullurl, data, timeout)
    519 for processor in self.process_response.get(protocol, []):
    520     meth = getattr(processor, meth_name)
--> 521     response = meth(req, response)
    523 return response

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:630, in HTTPErrorProcessor.http_response(self, request, response)
    627 # According to RFC 2616, "2xx" code indicates that the client's
    628 # request was successfully received, understood, and accepted.
    629 if not (200 <= code < 300):
--> 630     response = self.parent.error(
    631         'http', request, response, code, msg, hdrs)
    633 return response

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:553, in OpenerDirector.error(self, proto, *args)
    551     http_err = 0
    552 args = (dict, proto, meth_name) + args
--> 553 result = self._call_chain(*args)
    554 if result:
    555     return result

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:492, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
    490 for handler in handlers:
    491     func = getattr(handler, meth_name)
--> 492     result = func(*args)
    493     if result is not None:
    494         return result

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:745, in HTTPRedirectHandler.http_error_302(self, req, fp, code, msg, headers)
    742 fp.read()
    743 fp.close()
--> 745 return self.parent.open(new, timeout=req.timeout)

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:521, in OpenerDirector.open(self, fullurl, data, timeout)
    519 for processor in self.process_response.get(protocol, []):
    520     meth = getattr(processor, meth_name)
--> 521     response = meth(req, response)
    523 return response

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:630, in HTTPErrorProcessor.http_response(self, request, response)
    627 # According to RFC 2616, "2xx" code indicates that the client's
    628 # request was successfully received, understood, and accepted.
    629 if not (200 <= code < 300):
--> 630     response = self.parent.error(
    631         'http', request, response, code, msg, hdrs)
    633 return response

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:559, in OpenerDirector.error(self, proto, *args)
    557 if http_err:
    558     args = (dict, 'default', 'http_error_default') + orig_args
--> 559     return self._call_chain(*args)

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:492, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
    490 for handler in handlers:
    491     func = getattr(handler, meth_name)
--> 492     result = func(*args)
    493     if result is not None:
    494         return result

File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/urllib/request.py:639, in HTTPDefaultErrorHandler.http_error_default(self, req, fp, code, msg, hdrs)
    638 def http_error_default(self, req, fp, code, msg, hdrs):
--> 639     raise HTTPError(req.full_url, code, msg, hdrs, fp)

HTTPError: HTTP Error 404: Not Found

Read in and Plot the Data¶

Before following running the next cells, make sure you have created the following directories:

quicklooks

!mkdir quicklooks

radar_file = "houcsapr2cfrS2.a1/houcsapr2cfrS2.a1.20220807.183904.nc"

Plot one of the RHI scans¶

We read in the data corresponding to 7 August 2022 18:39:04 UTC, and plot a basic RadarDisplay which will automatically detect whether the plot is a vertical cross section (RHI or VPT), or a horizontal scan (PPI)

radar = pyart.io.read(radar_file)
display = pyart.graph.RadarDisplay(radar)
display.plot("reflectivity", 0)
plt.savefig(f"quicklooks/{Path(radar_file).stem}.png", dpi=200)
plt.xlim(0,45)
plt.ylim(0,12)
plt.show()
plt.close()

Define a function to grid the RHI data from polar (antenna) coordinates to a two-dimensional Caretsian grid¶

We use numba to vectorize the dist_func function to calculate the distance of each range gate from the radar. This makes our code run faster than simply executing this function for each gate in a for loop.

Next, we use the barnes function from the fastbarnes Python package to interpolate the radar fields such as equivalent reflectivity factor ( $Z_{H}$ ), differential_reflectivity ( $Z_{DR}$ ), and specific_differential_phase ( $K_{DP}$ ) to a uniform range-height Cartesian grid.

def grid_rhi(file,z_res = 100,rng_res = 100,z_limits = (0,15000),rng_limits = (0,109900.0),fields=None):
    
    """
    Input:
    -------
    file (str): Path of the RHI scan file that needs to be gridded.
    z_res (float): Vertical grid spacing of the Cartesian grid.
    rng_res (float): Horizontal grid spacing of the Cartesian grid.
    z_limits (tuple): Lower and upper height limits within which radar data needs to be gridded. 
    rng_limits (tuple): Lower and upper range limits within which radar data needs to be gridded. 
    
    Output:
    -------
    grid_ds (xarray Dataset): Xarray dataset containing radar fields on the Cartesian grid
    """
    
    z_pts = np.arange(z_limits[0],z_limits[1]+z_res,z_res)
    rng_pts = np.arange(rng_limits[0],rng_limits[1]+rng_res,rng_res)
   
    
    rhi = xr.open_dataset(file)
    radar = pyart.io.read(file)
    rhi = rhi.swap_dims({'time':'elevation'})

    lat = float(radar.latitude["data"])
    lon = float(radar.longitude["data"])

    grid_origin = (lat, lon)
    grid_origin_lat, grid_origin_lon = grid_origin

    grid_projection = {"proj": "aeqd", "_include_lon_0_lat_0": True}
    projparams = grid_projection.copy()

    if projparams.pop("_include_lon_0_lat_0", False):
        projparams["lon_0"] = grid_origin_lon
        projparams["lat_0"] = grid_origin_lat

    rg_loc = np.tile(radar.range['data'],len(radar.elevation['data'])).reshape(len(radar.elevation['data']),len(radar.range['data']))

    zg_loc = radar.gate_altitude["data"] - 12 # CSAPR2 antenna altitude = 12 m in this example

     
    # one of [ 'naive', 'radius', 'convolution', 'optimized_convolution' ]
    method = "optimized_convolution"

    # one of [ 0.25, 0.5, 1.0, 2.0, 4.0 ]
    sigma = 1.0

    # applies only to 'convolution' interpolations: one of [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 50 ]
    num_iter = 4
    
    if fields is None:
        fields = ['reflectivity','differential_reflectivity','specific_differential_phase',
                  'copol_correlation_coeff','mean_doppler_velocity','spectral_width']
        
    res_field = np.empty((len(z_pts),len(rng_pts),len(fields)))

    for j in range(len(fields)):

        data = deepcopy(np.array(radar.fields[fields[j]]['data']))
        # data = data.filled(np.nan)

        res_field[:,:,j] = barnes(np.asarray([rg_loc.ravel(),radar.gate_altitude['data'].ravel()]),
                           data.ravel(),
                           100,
                           np.asarray([0,0]),
                           100,
                           (len(z_pts),len(rng_pts)),
                           method=method,
                           num_iter = num_iter,
                           min_weight=0.0002
                          )
    

    data_dict = {}
    for k in range(len(fields)):
        data_dict[fields[k]] = (["z","range"],res_field[:,:,k])

    grid_ds = xr.Dataset(data_vars=data_dict,
                         coords=dict(z=(["z"], z_pts),
                                     range=(["range"],rng_pts),
                                    ),
                         attrs=dict(description="Gridded RHI data."),
                        )
        
    return grid_ds

Apply the gridding function to your radar RHI data file¶

Now, we are ready to grid the RHI data. The grid_rhi function requires the user to specify the vertical and horizontal grid spacing (z_res and rng_res, respectively), as well as the lower and upper limits of the Cartesian grid in the vertical and horizontal dimensions (z_limits and rng_limits, respectively). Custom fields of interest can be specified through the fields parameter. Otherwise, the function grids five radar fields by default i.e., ( $Z_{H}$ , $Z_{DR}$ , $\rho_{hv}$ , $K_{DP}$ , and $\sigma$ ).

grid_ds = grid_rhi(radar_file)

Grid and plot the RHI data¶

# Finally, plot the gridded reflectivity
fig,ax = plt.subplots()
grid_ds['reflectivity'].plot(vmin=0,vmax=70,cmap='HomeyerRainbow',ax=ax)
ax.set_xlim(0,55000)
ax.set_ylim(0,10000)
plt.show()

Summary¶

Within this example, we walked through how to access ARM data from a field campaign in Texas, plot a quick look of the RHI scan data, and grid our RHI data from native (polar) coordinates to a uniform range-height Caretsian grid.

Resources and References¶

ARM Data Discovery
TRACER Field Campaign
CSAPR Radar Data:
- Bharadwaj, N., Collis, S., Hardin, J., Isom, B., Lindenmaier, I., Matthews, A., & Nelson, D. C-Band Scanning ARM Precipitation Radar (CSAPR2CFR). Atmospheric Radiation Measurement (ARM) User Facility. Bharadwaj et al. (2021)
Py-ART:
- Helmus, J.J. & Collis, S.M., (2016). The Python ARM Radar Toolkit (Py-ART), a Library for Working with Weather Radar Data in the Python Programming Language. Journal of Open Research Software. 4(1), p.e25. DOI: Helmus & Collis (2016)
ACT:
- Adam Theisen, Ken Kehoe, Zach Sherman, Bobby Jackson, Alyssa Sockol, Corey Godine, Max Grover, Jason Hemedinger, Jenni Kyrouac, Maxwell Levin, Michael Giansiracusa (2022). The Atmospheric Data Community Toolkit (ACT). Zenodo. DOI: Theisen et al. (2022)

References¶

Zürcher, B. K. (2023). Fast approximate Barnes interpolation: illustrated by Python-Numba implementation fast-barnes-py v1.0. Geoscientific Model Development, 16(6), 1697–1711. 10.5194/gmd-16-1697-2023
Bharadwaj, N., Collis, S., Hardin, J., Isom, B., Lindenmaier, I., Matthews, A., Nelson, D., Feng, Y.-C., Rocque, M., Wendler, T., & Castro, V. (2021). C-Band Scanning ARM Precipitation Radar, 2nd Generation. Atmospheric Radiation Measurement (ARM) Archive, Oak Ridge National Laboratory (ORNL), Oak Ridge, TN (US); ARM Data Center, Oak Ridge National Laboratory (ORNL), Oak Ridge, TN (United States). 10.5439/1467901
Helmus, J. J., & Collis, S. M. (2016). The Python ARM Radar Toolkit (Py-ART), a Library for Working with Weather Radar Data in the Python Programming Language. Journal of Open Research Software, 4(1), 25. 10.5334/jors.119
Theisen, A., Kehoe, K., Sherman, Z., Jackson, B., Sockol, A., Godine, C., Grover, M., Hemedinger, J., Kyrouac, J., Levin, M., & Giansiracusa, M. (2022). ARM-DOE/ACT: ACT Release v1.1.9. Zenodo. 10.5281/ZENODO.6712343