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:114, in _error_wrapper(func, args, kwargs, retries)
113 try:
--> 114 return await func(*args, **kwargs)
115 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:1338, in AbstractFileSystem.open(self, path, mode, block_size, cache_options, compression, **kwargs)
1336 else:
1337 ac = kwargs.pop("autocommit", not self._intrans)
-> 1338 f = self._open(
1339 path,
1340 mode=mode,
1341 block_size=block_size,
1342 autocommit=ac,
1343 cache_options=cache_options,
1344 **kwargs,
1345 )
1346 if compression is not None:
1347 from fsspec.compression import compr
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:720, in S3FileSystem._open(self, path, mode, block_size, acl, version_id, fill_cache, cache_type, autocommit, size, requester_pays, cache_options, **kwargs)
717 if cache_type is None:
718 cache_type = self.default_cache_type
--> 720 return S3File(
721 self,
722 path,
723 mode,
724 block_size=block_size,
725 acl=acl,
726 version_id=version_id,
727 fill_cache=fill_cache,
728 s3_additional_kwargs=kw,
729 cache_type=cache_type,
730 autocommit=autocommit,
731 requester_pays=requester_pays,
732 cache_options=cache_options,
733 size=size,
734 )
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:2257, 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)
2255 self.details = s3.info(path)
2256 self.version_id = self.details.get("VersionId")
-> 2257 super().__init__(
2258 s3,
2259 path,
2260 mode,
2261 block_size,
2262 autocommit=autocommit,
2263 cache_type=cache_type,
2264 cache_options=cache_options,
2265 size=size,
2266 )
2267 self.s3 = self.fs # compatibility
2269 # 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:1912, in AbstractBufferedFile.__init__(self, fs, path, mode, block_size, autocommit, cache_type, cache_options, size, **kwargs)
1910 self.size = size
1911 else:
-> 1912 self.size = self.details["size"]
1913 self.cache = caches[cache_type](
1914 self.blocksize, self._fetch_range, self.size, **cache_options
1915 )
1916 else:
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/fsspec/spec.py:1925, in AbstractBufferedFile.details(self)
1922 @property
1923 def details(self):
1924 if self._details is None:
-> 1925 self._details = self.fs.info(self.path)
1926 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:1445, in S3FileSystem._info(self, path, bucket, key, refresh, version_id)
1443 if key:
1444 try:
-> 1445 out = await self._call_s3(
1446 "head_object",
1447 self.kwargs,
1448 Bucket=bucket,
1449 Key=key,
1450 **version_id_kw(version_id),
1451 **self.req_kw,
1452 )
1453 return {
1454 "ETag": out.get("ETag", ""),
1455 "LastModified": out.get("LastModified", ""),
(...) 1461 "ContentType": out.get("ContentType"),
1462 }
1463 except FileNotFoundError:
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:371, in S3FileSystem._call_s3(self, method, *akwarglist, **kwargs)
369 logger.debug("CALL: %s - %s - %s", method.__name__, akwarglist, kw2)
370 additional_kwargs = self._get_s3_method_kwargs(method, *akwarglist, **kwargs)
--> 371 return await _error_wrapper(
372 method, kwargs=additional_kwargs, retries=self.retries
373 )
File ~/micromamba/envs/radar-cookbook-dev/lib/python3.12/site-packages/s3fs/core.py:146, in _error_wrapper(func, args, kwargs, retries)
144 err = e
145 err = translate_boto_error(err)
--> 146 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