In [1]:
import matplotlib.pyplot as plt
from IPython.display import Image, display
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import gdal
import osr
import numpy as np
import netCDF4
import re
import pyart
import wradlib as wrl
import matplotlib
from matplotlib import cm
from IPython.display import Image, display
from mpl_toolkits.basemap import Basemap
get_ipython().magic(
'install_ext https://raw.github.com/cpcloud/ipython-\
autotime/master/autotime.py')
get_ipython().magic('load_ext autotime')
get_ipython().magic(
'install_ext https://raw.github.com/cpcloud/ipython-\
autotime/master/autotime.py')
get_ipython().magic('load_ext autotime')
In [2]:
"""
pyart.retrieve.beam_block_calc
=======================================
Calculates partial beam block(PBB) and cumulative beam block(CBB)
by using wradlib's beamblock and geotiff functions. PBB and CBB
are then used to created flags when a certain beam block fraction
is passed. Empty radar object is created using Py-ART and then
is filled with beam block data.
.. autosummary::
:toctreeL generated/
:template: dev_template.rst
beam_block
empty_radar_beam_block
beam_block_flag
_arrays_to_dict
_flags_to_dict
"""
import pyart
import numpy as np
import wradlib as wrl
def beam_block(radar, tif_name,
beam_width=1.0):
"""
Beam Block Calculation
Parameters
----------
radar : Radar
Radar object used.
tif_name : string
Name of geotiff file to use for the
calculation
Other Parameters
----------------
beam_width : float
Radar's beam width for calculation.
Default value is 1.0.
Returns
-------
pbb : array
Array of partial beam block fractions for each
gate in each ray.
cbb: array
Array of cumulative beam block fractions for
each gate in each ray.
References
----------
Bech, J., B. Codina, J. Lorente, and D. Bebbington,
2003: The sensitivity of single polarization weather
radar beam blockage correction to variability in the
vertical refractivity gradient. J. Atmos. Oceanic
Technol., 20, 845–855
Heistermann, M., Jacobi, S., and Pfaff, T., 2013:
Technical Note: An open source library for processing
weather radar data (wradlib), Hydrol. Earth Syst.
Sci., 17, 863-871, doi:10.5194/hess-17-863-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), p.e25.
DOI: http://doi.org/10.5334/jors.119
"""
radar.fields.clear()
_range = radar.range['data']
#beam_width = np.float(radar.instrument_parameters[
#'radar_beam_width_v']['data'])
beamradius = wrl.util.half_power_radius(_range, beam_width)
rasterfile = tif_name
data_raster = wrl.io.open_raster(rasterfile)
proj_raster = wrl.georef.wkt_to_osr(data_raster.GetProjection())
rastercoords, rastervalues = wrl.io.read_raster_data(rasterfile)
sitecoords = (np.float(radar.longitude['data']),
np.float(radar.latitude['data']),
np.float(radar.altitude['data']))
nrays = radar.nrays
nbins = radar.ngates
elev = radar.elevation['data']
azimuths = radar.azimuth['data']
rg, azg = np.meshgrid(_range, azimuths)
rg, eleg = np.meshgrid(_range, elev)
lon, lat, alt = wrl.georef.polar2lonlatalt_n(rg, azg,
eleg, sitecoords)
#lon = radar.gate_longitude['data']
#lat = radar.gate_latitude['data']
#alt = radar.gate_altitude['data']
x_pol, y_pol, = wrl.georef.reproject(
lon, lat, projection_target=proj_raster)
polcoords = np.dstack((x_pol, y_pol))
rlimits = (x_pol.min(), y_pol.min(),
x_pol.max(), y_pol.max())
# Clip the region inside our bounding box
ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
rastercoords = rastercoords[ind[1]:ind[3], ind[0]:ind[2], ...]
rastervalues = rastervalues[ind[1]:ind[3], ind[0]:ind[2]]
# Map rastervalues to polar grid points
polarvalues = wrl.ipol.cart2irregular_spline(
rastercoords, rastervalues, polcoords)
#bh = wrl.qual.beam_height_ft_doviak(
# _range, elev, degrees=True, re=6371000)
pbb = wrl.qual.beam_block_frac(polarvalues, alt, beamradius)
pbb = np.ma.masked_invalid(pbb)
maxindex = np.nanargmax(pbb, axis=1)
cbb = np.copy(pbb)
# Iterate over all beams
for ii, index in enumerate(maxindex):
premax = 0.
for jj in range(index):
# Only iterate to max index to make this faster
if pbb[ii, jj] > premax:
cbb[ii, jj] = pbb[ii, jj]
premax = cbb[ii, jj]
else:
cbb[ii, jj] = premax
# beyond max index, everything is max anyway
cbb[ii, index:] = pbb[ii, index]
return pbb, cbb
def beam_block_flag(pbb_all, cbb_all, pbb_threshold,
cbb_threshold):
""" Takes PBB and CBB arrays created from the
beam_block function and user chosen thresholds
to create and array of 1s and 0s, 1 is a flagged gate
where the fraction value is past the threshold. """
pbb_flags = np.empty_like(pbb_all)
pbb_flags[pbb_all > pbb_threshold] = True
pbb_flags[pbb_all < pbb_threshold] = False
cbb_flags = np.empty_like(cbb_all)
cbb_flags[cbb_all > cbb_threshold] = True
cbb_flags[cbb_all < cbb_threshold] = False
return pbb_flags, cbb_flags
def empty_radar_beam_block_rhi(ngates, rays_per_sweep, nsweeps,
lon, lat, alt, range_start,
gate_space, azimuth, elev_start,
elev_end):
""" Creates a radar object with no fields based on
user inputed dimensions. The empty radar is to then
be used to add PBB, CBB and the flags for both. """
radar = pyart.testing.make_empty_rhi_radar(
ngates, rays_per_sweep, nsweeps)
nrays = rays_per_sweep * nsweeps
radar.longitude['data'] = np.array([lon])
radar.latitude['data'] = np.array([lat])
radar.altitude['data'] = np.array([alt])
# radar.azimuth['data'] = np.linspace(0, 360, rays_per_sweep)
radar.range['data'] = np.linspace(
range_start, (ngates - 1)*gate_space + range_start, ngates)
radar.elevation['data'] = np.linspace(elev_start, elev_end, nrays)
radar.azimuth['data'] = np.array([azimuth] * nrays)
radar.fixed_angle['data'] = np.array([azimuth])
radar.metadata['instrument_name'] = 'beam_block_radar_object'
return radar
def empty_radar_beam_block_ppi(ngates, rays_per_sweep, nsweeps,
lon, lat, alt, range_start,
gate_space, elevations):
""" Creates a radar object with no fields based on
user inputed dimensions. The empty radar is to then
be used to add PBB, CBB and the flags for both. """
radar = pyart.testing.make_empty_ppi_radar(
ngates, rays_per_sweep, nsweeps)
radar.longitude['data'] = np.array([lon])
radar.latitude['data'] = np.array([lat])
radar.altitude['data'] = np.array([alt])
radar.azimuth['data'] = np.linspace(0, 360, rays_per_sweep)
radar.range['data'] = np.linspace(
range_start, (ngates - 1)*gate_space + range_start, ngates)
radar.fixed_angle['data'] = elevations
radar.elevation['data'] = np.array([elevations] * rays_per_sweep)
radar.metadata['instrument_name'] = 'beam_block_radar_object'
return radar
def _arrays_to_dict(pbb_all, cbb_all):
""" Function that takes the PBB and CBB arrays
and turns them into dictionaries to be used and added
to the pyart radar object. """
pbb_dict = {}
pbb_dict['coordinates'] = 'elevation, azimuth, range'
pbb_dict['units'] = 'unitless'
pbb_dict['data'] = pbb_all
pbb_dict['standard_name'] = 'partial_beam_block'
pbb_dict['long_name'] = 'Partial Beam Block Fraction'
pbb_dict['comment'] = 'Partial beam block fraction due to terrain'
cbb_dict = {}
cbb_dict['coordinates'] = 'elevation, azimuth, range'
cbb_dict['units'] = 'unitless'
cbb_dict['data'] = cbb_all
cbb_dict['standard_name'] = 'cumulative_beam_block'
cbb_dict['long_name'] = 'Cumulative Beam Block Fraction'
cbb_dict['comment'] = 'Cumulative beam block fraction due to terrain'
return pbb_dict, cbb_dict
def _flags_to_dict(pbb_flags, cbb_flags):
""" Function that takes the PBB_flag and CBB_flag
arrays and turns them into dictionaries to be used
and added to the pyart radar object. """
pbb_flag_dict = {}
pbb_flag_dict['units'] = 'unitless'
pbb_flag_dict['data'] = pbb_flags
pbb_flag_dict['standard_name'] = 'partial_beam_block_flag'
pbb_flag_dict['long_name'] = 'Partial Beam Block Flag'
pbb_flag_dict['comment'] = 'Partial beam block fraction flag, ' \
'1 for flagged values, 0 for non-flagged.'
cbb_flag_dict = {}
cbb_flag_dict['units'] = 'unitless'
cbb_flag_dict['data'] = cbb_flags
cbb_flag_dict['standard_name'] = 'cumulative_beam_block_flag'
cbb_flag_dict['long_name'] = 'Cumulative Beam Block Flag'
cbb_flag_dict['comment'] = 'Cumulative beam block fraction flag, ' \
'1 for flagged values, 0 for non-flagged.'
return pbb_flag_dict, cbb_flag_dict
In [3]:
def empty_radar_beam_block(filename):
""" Creates a radar object with no fields based on
user inputed dimensions. The empty radar is to then
be used to add PBB, CBB and the flags for both. """
radar = pyart.io.read(filename)
ngates = radar.ngates
nrays = radar.nrays
nsweeps = radar.nsweeps
empty_radar = pyart.testing.make_empty_rhi_radar(
ngates, 1, nsweeps)
empty_radar.nrays = nrays
empty_radar.time = radar.time
empty_radar.longitude = radar.longitude
empty_radar.latitude = radar.latitude
empty_radar.altitude = radar.altitude
empty_radar.range = radar.range
empty_radar.elevation = radar.elevation
empty_radar.azimuth = radar.azimuth
empty_radar.fixed_angle = radar.fixed_angle
empty_radar.metadata = radar.metadata
ray_start = []
ray_end = []
for i in range(len(radar.sweep_start_ray_index['data'])):
index_start = radar.sweep_start_ray_index['data'][i]
index_end = radar.sweep_end_ray_index['data'][i]
ray_start.append(index_start)
ray_end.append(index_end)
empty_radar.sweep_start_ray_index['data'] = np.array(ray_start)
empty_radar.sweep_end_ray_index['data'] = np.array(ray_end)
return empty_radar
In [4]:
# Examples
# --------
# >>> import pyart
# >>> radar = pyart.io.read('radar_file.nc')
# >>> gatefilter = pyart.correct.GateFilter(radar)
# >>> beam_block = radar.fields['partial_beam_block_flag']['data']
# >>> gatefilter.exclude_beam_block(beam_block)
def exclude_beam_block(self, beam_block, exclude_masked=True, op='or'):
"""
Exclude gates where a beam block is equal to True.
Parameters
----------
beam_block : numpy array
Boolean numpy array with same shape as a field array.
exclude_masked : bool, optional
True to filter masked values in the specified mask if it is
a masked array, False to include any masked values.
op : {'and', 'or', 'new'}
Operation to perform when merging the existing set of excluded
gates with the excluded gates from the current operation.
'and' will perform a logical AND operation, 'or' a logical OR,
and 'new' will replace the existing excluded gates with the one
generated here. 'or', the default for exclude methods, is
typically desired when building up a set of conditions for
excluding gates where the desired effect is to exclude gates which
meet any of the conditions. 'and', the default for include
methods, is typically desired when building up a set of conditions
where the desired effect is to include gates which meet any of the
conditions. Note that the 'and' method MAY results in including
gates which have previously been excluded because they were masked
or invalid.
"""
fdata = next(iter(self._radar.fields.values()))['data']
if beam_block.shape != fdata.shape:
raise ValueError("beam_block array must be the same size as a field.")
marked = np.array(beam_block, dtype='bool')
return self._merge(marked, op, exclude_masked)
In [5]:
azimuth = 147
radar = empty_radar_beam_block_rhi(983, 361, 1, -28.0257,
39.0916, 40.0, 0, 100,
azimuth, 0, 90)
In [6]:
tif_name = '/home/zsherman/beam_block/data/dtm_gra.tif'
In [7]:
pbb_all, cbb_all = beam_block(
radar, tif_name, 1.0)
In [8]:
print(pbb_all.max())
print(cbb_all.max())
In [9]:
pbb_flags, cbb_flags = beam_block_flag(
pbb_all, cbb_all, 0.2, 0.2)
In [10]:
pbb_dict, cbb_dict = _arrays_to_dict(
pbb_all, cbb_all)
pbb_flag_dict, cbb_flag_dict = _flags_to_dict(
pbb_flags, cbb_flags)
In [11]:
radar.add_field('partial_beam_block',
pbb_dict, replace_existing=True)
radar.add_field('cumulative_beam_block',
cbb_dict, replace_existing=True)
radar.add_field('partial_beam_block_flag',
pbb_flag_dict, replace_existing=True)
radar.add_field('cumulative_beam_block_flag',
cbb_flag_dict, replace_existing=True)
In [12]:
pyart.io.write_cfradial(
'/home/zsherman/beam_block/data/radar_object_rhi.nc', radar)
In [13]:
bb_radar = pyart.io.read(
'/home/zsherman/beam_block/data/radar_object_rhi.nc')
In [14]:
np.argwhere(bb_radar.fields[
'partial_beam_block_flag']['data'] == 1).shape
Out[14]:
In [15]:
bb_radar.fields[
'partial_beam_block_flag']['data']
Out[15]:
In [16]:
np.argwhere(bb_radar.fields['partial_beam_block']['data'] > 0.25).shape
Out[16]:
In [17]:
PBB = bb_radar.fields['partial_beam_block']['data']
CBB = bb_radar.fields['cumulative_beam_block']['data']
r = bb_radar.range['data']
th = bb_radar.elevation['data']
az = bb_radar.azimuth['data']
fig = plt.figure(figsize=(10, 7))
angle=520
ax, dem = wrl.vis.plot_rhi(CBB, r=r,
th=th,
cmap=plt.cm.PuRd)
ax.set_xlim(0, 12000)
ax.set_ylim(0, 1000)
ax.plot(0, 0, 'ro', )
ax.grid(True)
ax.annotate(' ARM ENA Site', (0, 0))
ax.set_title('Partial Beam Block 2.0 Degrees')
ax.set_xlabel("Meters")
ax.set_ylabel("Meters")
ax.set_axis_bgcolor('#E0E0E0')
plt.colorbar(dem, ax=ax)
plt.show()
In [18]:
azimuth = 147
radar = empty_radar_beam_block_rhi(1600, 540, 1, -28.0257,
39.0916, 40.0, 0, 10,
azimuth, 0.8, 90)
tif_name = '/home/zsherman/beam_block/data/dtm_gra.tif'
In [19]:
radar.latitude
Out[19]:
In [20]:
beam_width = 1.0
radar.fields.clear()
_range = radar.range['data']
#beam_width = np.float(radar.instrument_parameters[
#'radar_beam_width_v']['data'])
beamradius = wrl.util.half_power_radius(_range, beam_width)
rasterfile = tif_name
data_raster = wrl.io.open_raster(rasterfile)
proj_raster = wrl.georef.wkt_to_osr(data_raster.GetProjection())
rastercoords, rastervalues = wrl.io.read_raster_data(rasterfile)
sitecoords = (np.float(radar.longitude['data']),
np.float(radar.latitude['data']),
np.float(radar.altitude['data']))
nrays = radar.nrays
nbins = radar.ngates
elev = radar.elevation['data']
azimuths = radar.azimuth['data']
rg, azg = np.meshgrid(_range, azimuths)
rg, eleg = np.meshgrid(_range, elev)
lon, lat, alt = wrl.georef.polar2lonlatalt_n(rg, azg,
eleg, sitecoords)
#lon = radar.gate_longitude['data']
#lat = radar.gate_latitude['data']
#alt = radar.gate_altitude['data']
x_pol, y_pol = wrl.georef.reproject(
lon, lat, projection_target=proj_raster)
polcoords = np.dstack((x_pol, y_pol))
rlimits = (x_pol.min(), y_pol.min(),
x_pol.max(), y_pol.max())
# Clip the region inside our bounding box
ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
rastercoords = rastercoords[ind[1]:ind[3], ind[0]:ind[2], ...]
rastervalues = rastervalues[ind[1]:ind[3], ind[0]:ind[2]]
# Map rastervalues to polar grid points
polarvalues = wrl.ipol.cart2irregular_spline(
rastercoords, rastervalues, polcoords)
pbb = wrl.qual.beam_block_frac(polarvalues, alt, beamradius)
pbb = np.ma.masked_invalid(pbb)
maxindex = np.nanargmax(pbb, axis=1)
cbb = np.copy(pbb)
# Iterate over all beams
for ii, index in enumerate(maxindex):
premax = 0.
for jj in range(index):
# Only iterate to max index to make this faster
if pbb[ii, jj] > premax:
cbb[ii, jj] = pbb[ii, jj]
premax = pbb[ii, jj]
else:
cbb[ii, jj] = premax
# beyond max index, everything is max anyway
cbb[ii, index:] = pbb[ii, index]
In [21]:
r = radar.range['data']
th = radar.elevation['data']
az = radar.azimuth['data']
fig = plt.figure(figsize=(10, 7))
ax, dem = wrl.vis.plot_rhi(cbb, r=r,
th=th, vmin=0, vmax=1,
cmap=plt.cm.PuRd)
ax.fill_between(r, 0.0,
polarvalues[0, :]-40,
color='0.45')
ax.set_xlim(0, 11000)
ax.set_ylim(0, 1000)
ax.plot(0, 0, 'ro', markersize=15)
ax.grid(True)
ax.annotate(' ARM ENA Site', (0, 50))
ax.set_title(
'Cumulative Beam Block Azimuth ' + str(azimuth) +
' Degrees \n Elevation 0.8 degrees to 90 degrees')
ax.set_xlabel("Meters From Radar")
ax.set_ylabel("Meters Above Radar")
ax.set_axis_bgcolor('#E0E0E0')
plt.colorbar(dem, ax=ax)
plt.savefig('/home/zsherman/beam_block/images/rhi_az147_gate10m.png',
bbox_inches='tight')
plt.show()
In [22]:
mdv = pyart.io.read(
'/home/zsherman/Downloads/example_mdv_rhi.mdv')
In [23]:
mdv.nsweeps
Out[23]:
In [24]:
mdv.azimuth['data'].shape
Out[24]:
In [25]:
mdv.elevation['data'].shape
Out[25]:
In [26]:
mdv.elevation['data']
Out[26]:
In [27]:
mdv.range
Out[27]:
In [28]:
float(mdv.instrument_parameters['radar_beam_width_v']['data'])
Out[28]:
In [29]:
sitecoords = (np.float(mdv.longitude['data']),
np.float(mdv.latitude['data']),
np.float(mdv.altitude['data']))
In [30]:
coord = wrl.georef.sweep_centroids(mdv.nrays, 100,
mdv.ngates, mdv.elevation['data'][:, np.newaxis])
In [31]:
lona, lata, alta = wrl.georef.polar2lonlatalt_n(
mdv.get_gate_x_y_z(0)[0], mdv.get_gate_x_y_z(0)[1],
mdv.get_gate_x_y_z(0)[2], sitecoords)
In [32]:
lata
Out[32]:
In [33]:
lona
Out[33]:
In [34]:
alta
Out[34]: