An echo top is the radar indicated top of an area of precipitation. This notebook demonstrates how to calculate the echo top height (ETH) in a NEXRAD PPI volume scan to determine the maximum elevation angle at which a certain reflectivity threshold is exceeded.
This example uses the echo top height (ETH) calculation code written by Valentin Louf, available at this github repository.
Overview¶
The notebook applies the modified ETH algorithm proposed by Lakshmanan et al. (2013) to a NEXRAD PPI volume scan.
The modified algorithm comprises these steps:
Find the maximum elevation angle () where reflectivity () exceeds the echo-top reflectivity threshold. If is not the highest elevation scan in the virtual volume, obtain the reflectivity value () at the next higher elevation angle (). Then, the echo-top height is given by the height of the radar beam at an elevation angle:
where is the threshold value (e.g., 0 dBZ, 18 dBZ) used to compute the echo top.
If is the highest elevation scan available, set , where is the half-power beamwidth. This condition is met far away from the radar if higher-elevation scans have shorter ranges than a base “surveillance” scan and very close to the radar if the highest-elevation scan does not sample the top of the cloud. Under these circumstances, is set to be the top of the beam containing dBZ ; that is, the traditional echo-top algorithm is followed when there are no data available from a higher-elevation scan.
Prerequisites¶
| Concepts | Importance | Notes |
|---|---|---|
| Matplotlib Basics | Required | Basic plotting |
| Intro to Cartopy | Required | Geospatial plotting |
| Py-ART Basics | Required | IO/Visualization |
| Numba | Helpful | Familiarity with vectorization/optimization of Python functions |
Time to learn: 25 minutes
Imports¶
import pyart
import numpy as np
from echotop import cloud_top_height, grid_cloud_top
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
## 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 as part of the Atmospheric
## Radiation Measurement (ARM) Climate Research Facility, an Office of
## Science 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
Define our function to compute ETH on a uniform x-y grid¶
The input file can be a Next Generation Weather Radar (NEXRAD) archive file from Amazon Web Services. We will remotely access this file and use a reflectivity threshold eth_thld to define our echo tops.
def compute_eth(infile, eth_thld=0):
"""
Compute the Echo Top Height on a grid
Parameters
==========
infile (str): Filename of NEXRAD Level 2 Archive file. The files hosted by
at the NOAA National Climate Data Center [1]_ as well as on the
UCAR THREDDS Data Server [2]_ have been tested. Other NEXRAD
Level 2 Archive files may or may not work. Message type 1 file
and message type 31 files are supported.
eth_thld (float): Reflectivity threshold for which we want to compute
the echo top height.
Returns:
========
cth_grid: ndarray <x, y>
Echo top height on a grid of dimension (x, y).
References
----------
.. [1] http://www.ncdc.noaa.gov/
.. [2] http://thredds.ucar.edu/thredds/catalog.html
"""
# Reading NEXRAD L2 data stored on AWS cloud
radar = pyart.io.read_nexrad_archive(infile)
r = radar.range['data']
azimuth = radar.azimuth['data']
elevation = radar.elevation['data']
refl = np.array(radar.fields['reflectivity']['data'])
st_sweep = radar.sweep_start_ray_index['data']
ed_sweep = radar.sweep_end_ray_index['data']
# Compute ETH. The 'echotop' package uses @jit decorator to optimize
# the 'cloud_top_height' function
cth = cloud_top_height(r, azimuth, elevation, st_sweep, ed_sweep, refl, eth_thld=eth_thld)
# Grid data
th = 450 - azimuth[slice(st_sweep[0], ed_sweep[0] + 1)]
th[th < 0] += 360
R, A = np.meshgrid(r, th)
x = R * np.cos(np.pi * A / 180)
y = R * np.sin(np.pi * A / 180)
xgrid = np.arange(-MAX_RANGE, MAX_RANGE + RANGE_STEP / 2, RANGE_STEP).astype(np.int32)
[X, Y] = np.meshgrid(xgrid, xgrid)
cth_grid = grid_cloud_top(
cth, x, y, X, Y, nnearest=24, maxdist=2500
) # nearest=24 should be enough to sample out to 2500m on a 1000m grid
cth_grid = np.ma.masked_invalid(cth_grid).astype(np.int32).filled(FILLVALUE)
return cth_gridRead and plot reflectivity and velocity fields for a sample file¶
aws_nexrad_level2_file = (
"s3://noaa-nexrad-level2/2022/03/22/KHGX/KHGX20220322_120125_V06"
)
radar = pyart.io.read_nexrad_archive(aws_nexrad_level2_file)
fig = plt.figure(figsize=(12, 4))
display = pyart.graph.RadarMapDisplay(radar)
ax = plt.subplot(121, projection=ccrs.PlateCarree())
display.plot_ppi_map(
"reflectivity",
sweep=0,
ax=ax,
colorbar_label="Equivalent Relectivity ($Z_{e}$) \n (dBZ)",
vmin=-20,
vmax=60,
)
ax = plt.subplot(122, projection=ccrs.PlateCarree())
display.plot_ppi_map(
"velocity",
sweep=1,
ax=ax,
colorbar_label="Radial Velocity ($V_{r}$) \n (m/s)",
vmin=-70,
vmax=70,
)---------------------------------------------------------------------------
ClientError Traceback (most recent call last)
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:115, in _error_wrapper(func, args, kwargs, retries)
114 try:
--> 115 return await func(*args, **kwargs)
116 except S3_RETRYABLE_ERRORS as e:
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/aiobotocore/context.py:36, in with_current_context.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
35 await resolve_awaitable(hook())
---> 36 return await func(*args, **kwargs)
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/aiobotocore/client.py:424, in AioBaseClient._make_api_call(self, operation_name, api_params)
423 error_class = self.exceptions.from_code(error_code)
--> 424 raise error_class(parsed_response, operation_name)
425 else:
ClientError: An error occurred (403) when calling the HeadObject operation: Forbidden
The above exception was the direct cause of the following exception:
PermissionError Traceback (most recent call last)
Cell In[3], line 5
1 aws_nexrad_level2_file = (
2 "s3://noaa-nexrad-level2/2022/03/22/KHGX/KHGX20220322_120125_V06"
3 )
----> 5 radar = pyart.io.read_nexrad_archive(aws_nexrad_level2_file)
8 fig = plt.figure(figsize=(12, 4))
9 display = pyart.graph.RadarMapDisplay(radar)
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/pyart/io/nexrad_archive.py:126, in read_nexrad_archive(filename, field_names, additional_metadata, file_field_names, exclude_fields, include_fields, delay_field_loading, station, scans, linear_interp, storage_options, **kwargs)
115 filemetadata = FileMetadata(
116 "nexrad_archive",
117 field_names,
(...) 121 include_fields,
122 )
124 # open the file and retrieve scan information
125 nfile = NEXRADLevel2File(
--> 126 prepare_for_read(filename, storage_options=storage_options)
127 )
128 scan_info = nfile.scan_info(scans)
130 # time
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/pyart/io/common.py:42, in prepare_for_read(filename, storage_options)
39 return filename
41 # look for compressed data by examining the first few bytes
---> 42 fh = fsspec.open(filename, mode="rb", compression="infer", **storage_options).open()
43 magic = fh.read(3)
44 fh.close()
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/core.py:147, in OpenFile.open(self)
140 def open(self):
141 """Materialise this as a real open file without context
142
143 The OpenFile object should be explicitly closed to avoid enclosed file
144 instances persisting. You must, therefore, keep a reference to the OpenFile
145 during the life of the file-like it generates.
146 """
--> 147 return self.__enter__()
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/core.py:105, in OpenFile.__enter__(self)
102 mode = self.mode.replace("t", "").replace("b", "") + "b"
104 try:
--> 105 f = self.fs.open(self.path, mode=mode)
106 except FileNotFoundError as e:
107 if has_magic(self.path):
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/spec.py:1349, in AbstractFileSystem.open(self, path, mode, block_size, cache_options, compression, **kwargs)
1347 else:
1348 ac = kwargs.pop("autocommit", not self._intrans)
-> 1349 f = self._open(
1350 path,
1351 mode=mode,
1352 block_size=block_size,
1353 autocommit=ac,
1354 cache_options=cache_options,
1355 **kwargs,
1356 )
1357 if compression is not None:
1358 from fsspec.compression import compr
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:733, in S3FileSystem._open(self, path, mode, block_size, acl, version_id, fill_cache, cache_type, autocommit, size, requester_pays, cache_options, **kwargs)
730 if cache_type is None:
731 cache_type = self.default_cache_type
--> 733 return S3File(
734 self,
735 path,
736 mode,
737 block_size=block_size,
738 acl=acl,
739 version_id=version_id,
740 fill_cache=fill_cache,
741 s3_additional_kwargs=kw,
742 cache_type=cache_type,
743 autocommit=autocommit,
744 requester_pays=requester_pays,
745 cache_options=cache_options,
746 size=size,
747 )
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:2295, in S3File.__init__(self, s3, path, mode, block_size, acl, version_id, fill_cache, s3_additional_kwargs, autocommit, cache_type, requester_pays, cache_options, size)
2293 self.details = s3.info(path)
2294 self.version_id = self.details.get("VersionId")
-> 2295 super().__init__(
2296 s3,
2297 path,
2298 mode,
2299 block_size,
2300 autocommit=autocommit,
2301 cache_type=cache_type,
2302 cache_options=cache_options,
2303 size=size,
2304 )
2305 self.s3 = self.fs # compatibility
2307 # when not using autocommit we want to have transactional state to manage
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/spec.py:1923, in AbstractBufferedFile.__init__(self, fs, path, mode, block_size, autocommit, cache_type, cache_options, size, **kwargs)
1921 self.size = size
1922 else:
-> 1923 self.size = self.details["size"]
1924 self.cache = caches[cache_type](
1925 self.blocksize, self._fetch_range, self.size, **cache_options
1926 )
1927 else:
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/spec.py:1936, in AbstractBufferedFile.details(self)
1933 @property
1934 def details(self):
1935 if self._details is None:
-> 1936 self._details = self.fs.info(self.path)
1937 return self._details
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/asyn.py:118, in sync_wrapper.<locals>.wrapper(*args, **kwargs)
115 @functools.wraps(func)
116 def wrapper(*args, **kwargs):
117 self = obj or args[0]
--> 118 return sync(self.loop, func, *args, **kwargs)
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/asyn.py:103, in sync(loop, func, timeout, *args, **kwargs)
101 raise FSTimeoutError from return_result
102 elif isinstance(return_result, BaseException):
--> 103 raise return_result
104 else:
105 return return_result
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/asyn.py:56, in _runner(event, coro, result, timeout)
54 coro = asyncio.wait_for(coro, timeout=timeout)
55 try:
---> 56 result[0] = await coro
57 except Exception as ex:
58 result[0] = ex
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:1471, in S3FileSystem._info(self, path, bucket, key, refresh, version_id)
1469 if key:
1470 try:
-> 1471 out = await self._call_s3(
1472 "head_object",
1473 self.kwargs,
1474 Bucket=bucket,
1475 Key=key,
1476 **version_id_kw(version_id),
1477 **self.req_kw,
1478 )
1479 return {
1480 "ETag": out.get("ETag", ""),
1481 "LastModified": out.get("LastModified", ""),
(...) 1487 "ContentType": out.get("ContentType"),
1488 }
1489 except FileNotFoundError:
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:384, in S3FileSystem._call_s3(self, method, *akwarglist, **kwargs)
382 logger.debug("CALL: %s - %s - %s", method.__name__, akwarglist, kw2)
383 additional_kwargs = self._get_s3_method_kwargs(method, *akwarglist, **kwargs)
--> 384 return await _error_wrapper(
385 method, kwargs=additional_kwargs, retries=self.retries
386 )
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:147, in _error_wrapper(func, args, kwargs, retries)
145 err = e
146 err = translate_boto_error(err)
--> 147 raise err
PermissionError: ForbiddenDefine some global constants and compute the ETH on a horizontally uniform grid¶
These constants are required by the compute_eth function to grid the data from polar coordinates to a horizontally uniform grid. These three constants are defined as:
FILLVALUE (Value to replace missing data)RANGE_STEP (Uniform horizontal grid spacing in x and y dimensions)MAX_RANGE (Maximum range up to which ETH to be calculated from gridded data)
FILLVALUE: int = -9999
RANGE_STEP: int = 1000
MAX_RANGE: float = 250e3
cth_grid = compute_eth(aws_nexrad_level2_file, eth_thld=20)Plot the gridded echo/cloud top height using matplotlib¶
p = plt.pcolormesh(cth_grid,
vmin=0,
vmax=15000,
cmap='ChaseSpectral')
plt.colorbar(mappable=p)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¶
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)
Echo-top height algorithm:
Lakshmanan, V., Hondl, K., Potvin, C. K., & Preignitz, D. (2013). An improved method for estimating radar echo-top height. Weather and Forecasting, 28(2), 481-488. DOI: Lakshmanan et al. (2013)
- 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
- Lakshmanan, V., Hondl, K., Potvin, C. K., & Preignitz, D. (2013). An Improved Method for Estimating Radar Echo-Top Height. Weather and Forecasting, 28(2), 481–488. 10.1175/waf-d-12-00084.1