In [839]:
import gdal
import osr
import numpy as np
import netCDF4
import re
import pyart
import wradlib as wrl
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation
from mpl_toolkits.basemap import Basemap
from matplotlib import cm
from mpl_toolkits.basemap import Basemap
from IPython.display import Image, display
get_ipython().magic(
'install_ext https://raw.github.com/cpcloud/ipython-\
autotime/master/autotime.py')
get_ipython().magic('load_ext autotime')
In [840]:
"""
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_ppi
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']
beamradius = wrl.util.half_power_radius(_range, beam_width)
rasterfile = tif_name
data_raster = wrl.io.open_raster(rasterfile)
rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(
data_raster, nodata=None)
sitecoords = (np.float(radar.longitude['data']),
np.float(radar.latitude['data']),
np.float(radar.altitude['data']))
pbb_arrays = []
cbb_arrays = []
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] + 1
#nrays = index_end - index_start + 1
#nbins = radar.ngates
#range_res = (np.max(_range) - np.min(_range)) / (nbins - 1)
#range_res = _range[1] - _range[0]
#range_res = wrl.georef._get_range_resolution(_range)
#elev = radar.fixed_angle['data'][i]
elevs = radar.elevation['data'][index_start:index_end]
azimuths = radar.azimuth['data'][index_start:index_end]
rg, azg = np.meshgrid(_range, azimuths)
rg, eleg = np.meshgrid(_range, elevs)
lon, lat, alt = wrl.georef.polar2lonlatalt_n(
rg, azg, eleg, sitecoords)
#lon = radar.gate_longitude['data'][index_start:index_end +1]
#lat = radar.gate_latitude['data'][index_start:index_end + 1]
#alt = radar.gate_altitude['data'][index_start:index_end + 1]
x_pol, y_pol = wrl.georef.reproject(
lon, lat, projection_target=proj)
polcoords = np.dstack((x_pol, y_pol))
rlimits = (x_pol.min(), y_pol.min(), x_pol.max(), y_pol.max())
ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
rastercoords = rastercoords[0:ind[3], ind[0]:ind[2], ...]
rastervalues = rastervalues[0: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)
pbb_arrays.append(pbb)
cbb = wrl.qual.cum_beam_block_frac(pbb)
cbb_arrays.append(cbb)
pbb_all = np.ma.concatenate(pbb_arrays)
cbb_all = np.ma.concatenate(cbb_arrays)
return pbb_all, cbb_all
def beam_block_flag(pbb_all, cbb_all, no_block_thresh=0.01,
complete_block_thresh=0.95):
""" 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 > complete_block_thresh] = 2
pbb_flags[(pbb_all > no_block_thresh) & (pbb_all < complete_block_thresh)] = 1
pbb_flags[pbb_all < no_block_thresh] = 0
cbb_flags = np.empty_like(cbb_all)
cbb_flags[cbb_all > complete_block_thresh] = 2
cbb_flags[(cbb_all > no_block_thresh) & (cbb_all < complete_block_thresh)] = 1
cbb_flags[cbb_all < no_block_thresh] = 0
return pbb_flags, cbb_flags
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_flags'
pbb_flag_dict['long_name'] = 'Partial Beam Block Flags'
pbb_flag_dict['notes'] = '0: No Blockage, 1: Partial Blockage,' \
' 2: Complete Blockage'
cbb_flag_dict = {}
cbb_flag_dict['units'] = 'unitless'
cbb_flag_dict['data'] = cbb_flags
cbb_flag_dict['standard_name'] = 'cumulative_beam_block_flags'
cbb_flag_dict['long_name'] = 'Cumulative Beam Block Flags'
cbb_flag_dict['notes'] = '0: No Blockage, 1: Partial Blockage,' \
' 2: Complete Blockage'
return pbb_flag_dict, cbb_flag_dict
In [363]:
# 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 [859]:
# sitecoords = (-28.0257, 39.0916, 40.0)
elevations = np.array([0.5, 1.0, 2.0, 3.0, 4.0, 10.0, 11.0, 15.0,
20.0, 30.0, 40.0])
radar = empty_radar_beam_block_ppi(1800, 1090, 1, -28.0257,
39.0916, 40.0, 0, 10,
elevations=2.0)
tif_name = '/home/zsherman/beam_block/data/dtm_gra.tif'
In [84]:
elevations = np.array([0.5, 1.0, 2.0, 3.0, 4.0, 10.0, 11.0, 15.0,
20.0, 30.0, 40.0])
radar = empty_radar_beam_block_ppi(1832, 720, 1, -122.4959,
48.1946, 195.0, 2125, 250,
elevations=0.4833)
tif_name = '/home/zsherman/seattle.tif'
In [741]:
#pyart.io.write_cfradial(
# '/home/zsherman/sample_radar.nc', radar)
In [860]:
pbb_all, cbb_all = beam_block(
radar, tif_name, 1.0)
In [861]:
print(pbb_all.max())
print(cbb_all.max())
In [744]:
pbb_all.shape
Out[744]:
In [706]:
wrl.util.get_raster_origin(rastercoords)
Out[706]:
In [698]:
pract = {}
pract['notes'] = '0: No Blockage, 1: Partial Blockage, 2: Complete Blockage'
In [699]:
for pair_str in pract['notes'].split(','):
print(pair_str)
pract.update({pair_str.split(':')[1]:int(pair_str.split(':')[0])})
In [700]:
pract
Out[700]:
In [511]:
tif_name = '/home/zsherman/beam_block/data/dtm_gra.tif'
sitecoords = (-28.0257, 39.0916, 40.0)
nrays = 1200
nbins = 1000
el = 2.0
bw = 1.0
range_res = 100.0
r = np.arange(nbins) * range_res
beamradius = wrl.util.half_power_radius(r, bw)
rasterfile = tif_name
data_raster = wrl.io.open_raster(rasterfile)
rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(
data_raster)
print('rastercoords.shape = ' ,rastercoords.shape)
print('rastervalues.shape = ', rastervalues.shape)
coord = wrl.georef.sweep_centroids(nrays, range_res, nbins, el)
lon, lat, alt = wrl.georef.polar2lonlatalt_n(
coord[..., 0], np.degrees(coord[..., 1]),
coord[..., 2], sitecoords)
x_pol, y_pol = wrl.georef.reproject(
lon, lat, projection_target=proj)
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)
print('ind = ', ind)
rastercoords = rastercoords[ind[1]:ind[3], ind[0]:ind[2], ...]
rastervalues = rastervalues[ind[1]:ind[3], ind[0]:ind[2]]
print('rastercoords values = ', rastercoords)
print('rastervalues values = ', rastervalues)
polarvalues = wrl.ipol.cart2irregular_spline(rastercoords, rastervalues,
polcoords, order=3,
prefilter=False)
In [504]:
ind
Out[504]:
In [517]:
rastercoords.shape
Out[517]:
In [518]:
rastervalues.shape
Out[518]:
In [506]:
ind
Out[506]:
In [862]:
pbb_flags, cbb_flags = beam_block_flag(
pbb_all, cbb_all)
In [863]:
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 [864]:
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_flags',
pbb_flag_dict, replace_existing=True)
radar.add_field('cumulative_beam_block_flags',
cbb_flag_dict, replace_existing=True)
In [865]:
pyart.io.write_cfradial(
'/home/zsherman/flag_radar.nc', radar)
In [866]:
bb_radar = pyart.io.read(
'/home/zsherman/flag_radar.nc')
In [825]:
bb_radar.fields['cumulative_beam_block_flags']['data'].max()
Out[825]:
In [749]:
bb_radar = pyart.io.read(
'/home/zsherman/beam_block/data/radar_object_pbb_cbb_all_sweeps')
In [850]:
bb_radar.fields['partial_beam_block']['data']
Out[850]:
In [868]:
first_sweep = bb_radar.extract_sweeps([0])
In [869]:
first_sweep.azimuth['data'] = np.linspace(0, 360, 1090)
In [853]:
PBB = first_sweep.fields['partial_beam_block']['data']
In [649]:
mask = np.argwhere(PBB > 0.0)[:, 0]
mask
Out[649]:
In [633]:
loc = np.argwhere(PBB > 0.01)
type(loc)
Out[633]:
In [92]:
from collections import OrderedDict
list_two = list(OrderedDict.fromkeys(mask))
In [93]:
list_one = np.linspace(0, 360, 990)
In [94]:
list_three = list_one[list_two]
In [95]:
list_three_vals = [x for x in list_one if x not in list_three]
In [96]:
print(str(np.round(list_three_vals, decimals=3)))
In [97]:
# az_strings = (first_sweep.azimuth['data'][list_three_vals])
az_strings = list_three_vals
In [98]:
new_string = np.array(np.round(az_strings, decimals=3))
In [99]:
string = (', '.join([str(x) for x in new_string]))
In [100]:
spl = string.split(",")
new_fig_text = (','"\n".join([",".join(spl[i:i+15]) for i in range(0, len(spl), 15)]))
In [101]:
PBB = first_sweep.fields['partial_beam_block']['data']
CBB = first_sweep.fields['cumulative_beam_block']['data']
r = first_sweep.range['data']
az = first_sweep.azimuth['data']
el = first_sweep.azimuth['data']
azimuths = first_sweep.azimuth['data']
fig = plt.figure(figsize=(14, 9))
ax, dem = wrl.vis.plot_ppi(PBB, r=r,
az=az,
cmap=plt.cm.PuRd)
ax.set_xlim(-8000, 10000)
ax.set_ylim(-12000, 5000)
ax.plot(0, 0, 'ro', )
for az_string in az_strings:
ax.plot(
[0, np.sin(np.radians(az_string))*1e5],
[0, np.cos(np.radians(az_string))*1e5],
"r-", linewidth=0.2)
ax.grid(True)
ax.annotate(' ARM ENA Site', (r[0], az[1]))
ticks = (ax.get_xticks()/1000).astype(np.int)
ax.set_xticklabels(ticks)
ticks = (ax.get_yticks()/1000).astype(np.int)
ax.set_yticklabels(ticks)
ax.set_title('Partial Beam Block 0.5 Degrees Elevation')
ax.set_xlabel("Kilometers")
ax.set_ylabel("Kilometers")
ax.set_axis_bgcolor('#E0E0E0')
plt.colorbar(dem, ax=ax)
plt.figtext(.85, .04, ('Azimuths No Blockage:\n ' + new_fig_text))
plt.savefig(
'/home/zsherman/beam_block/images/graciaso_0.5elev_azimuths_no_blockage.png',
bbox_inches='tight')
plt.show()
In [559]:
def convertXY(xy_source, inproj, outproj):
# function to convert coordinates
shape = xy_source[0, :, :].shape
size = xy_source[0, :, :].size
# the ct object takes and returns pairs of x,y, not 2d grids
# so the the grid needs to be reshaped (flattened) and back.
ct = osr.CoordinateTransformation(inproj, outproj)
xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T))
xx = xy_target[:, 0].reshape(shape)
yy = xy_target[:, 1].reshape(shape)
return xx, yy
In [870]:
ds = gdal.Open("/home/zsherman/beam_block/data/dtm_gra.tif")
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
xres = gt[1]
yres = gt[5]
# get the edge coordinates and add half the resolution
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
ymax = gt[3] - yres * 0.5
ds = None
# create a grid of xy coordinates in the original projection
xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
In [871]:
llcrnrlat=39
urcrnrlon=-27.9
llcrnrlon=-28.1
urcrnrlat=39.1
lats = np.arange(llcrnrlat, urcrnrlat, .04)
lons = np.arange(llcrnrlon, urcrnrlon, .04)
fig = plt.figure(figsize=(12, 7))
m = Basemap(projection='lcc', lon_0=-28, lat_0=39, resolution='c',
llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat,
urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat)
inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)
# Get the target projection from the basemap object
outproj = osr.SpatialReference()
outproj.ImportFromProj4(m.proj4string)
# Convert from source projection to basemap projection
xx, yy = convertXY(xy_source, inproj, outproj)
elev = data.T
the_sea = data.T < 0.1
masked_data = np.ma.masked_where(the_sea, elev)
In [874]:
title = 'Cumulative Beam Blockage Over Terrain 2.0 Degree Elevation'
fig = plt.figure(figsize=(14, 7))
my_display = pyart.graph.RadarMapDisplay(first_sweep)
my_display.plot_ppi_map('cumulative_beam_block', basemap=m,
lat_lines=lats, lon_lines=lons,
vmin=0, vmax=1, cmap=plt.cm.PuRd,
title=title)
my_display.plot_point(-28.0257, 39.0916,
symbol='ro')
plt.annotate('ARM ENA Site', m(-28.0217, 39.0916), color='k')
CS = m.contour(xx, yy, masked_data,
levels=np.arange(0, 400, 100),
linewidths=1, cmap=plt.cm.terrain)
plt.clabel(CS, inline=1, fontsize=10, fmt='%1.0f')
plt.savefig(
'/home/zsherman/beam_block/images/cbb_graciosa_contour.png',
bbox_inches='tight')
plt.show()
In [890]:
ticks = ([0.33, 1, 1.67])
ticklabs = (['No\nBlockage', 'Partial\nBlockage', 'Complete\nBlockage'])
title = 'Partial Beam Blockage Flag 2.0 Degree Elevation'
fig = plt.figure(figsize=(14, 7))
my_display = pyart.graph.RadarMapDisplay(first_sweep)
lab_colors=['grey','yellow', 'red']
cmap = matplotlib.colors.ListedColormap(lab_colors)
my_display.plot_ppi_map('partial_beam_block_flags', basemap=m,
lat_lines=lats, lon_lines=lons,
cmap=cmap, title=title, ticklabs=ticklabs,
ticks=ticks)
my_display.plot_point(-28.0257, 39.0916,
symbol='ro')
plt.annotate('ARM ENA Site', m(-28.0217, 39.0916), color='k')
CS = m.contour(xx, yy, masked_data,
levels=np.arange(0, 400, 100),
linewidths=1, cmap=plt.cm.terrain)
plt.clabel(CS, inline=1, fontsize=10, fmt='%1.0f')
plt.savefig(
'/home/zsherman/beam_block/images/pbb_flag_2degree_elev.png',
bbox_inches='tight')
plt.show()
In [889]:
title = 'Cumulative Beam Blockage Flag 2.0 Degree Elevation'
fig = plt.figure(figsize=(14, 7))
my_display = pyart.graph.RadarMapDisplay(first_sweep)
lab_colors=['grey','yellow', 'red']
cmap = matplotlib.colors.ListedColormap(lab_colors)
my_display.plot_ppi_map('cumulative_beam_block_flags', basemap=m,
lat_lines=lats, lon_lines=lons,
cmap=cmap, title=title, ticklabs=ticklabs,
ticks=ticks)
my_display.plot_point(-28.0257, 39.0916,
symbol='ro')
plt.annotate('ARM ENA Site', m(-28.0217, 39.0916), color='k')
CS = m.contour(xx, yy, masked_data,
levels=np.arange(0, 400, 100),
linewidths=1, cmap=plt.cm.terrain)
plt.clabel(CS, inline=1, fontsize=10, fmt='%1.0f')
plt.savefig(
'/home/zsherman/beam_block/images/cbb_flag_2degree_elev.png',
bbox_inches='tight')
plt.show()
In [558]:
tif_name = '/home/zsherman/beam_block/data/dtm_gra.tif
elevations = np.array([0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0])
radars = []
for elevation in elevations:
radar = empty_radar_beam_block_ppi(
1800, 1090, 1, -28.0257, 39.0916, 40.0, 0, 10,
elevations=elevation)
pbb_all, cbb_all = beam_block(
radar, tif_name, 1.0)
pbb_dict, cbb_dict = _arrays_to_dict(
pbb_all, cbb_all)
radar.add_field('partial_beam_block',
pbb_dict, replace_existing=True)
radar.add_field('cumulative_beam_block',
cbb_dict, replace_existing=True)
radar.azimuth['data'] = np.linspace(0, 360, 1090)
radars += [radar]
In [560]:
llcrnrlat=39
urcrnrlon=-27.9
llcrnrlon=-28.1
urcrnrlat=39.1
lats = np.arange(llcrnrlat, urcrnrlat, .02)
lons = np.arange(llcrnrlon, urcrnrlon, .04)
In [561]:
m = Basemap(projection='lcc', lon_0=-28, lat_0=39, resolution='c',
llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat,
urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat)
In [578]:
def animate(nframe):
plt.clf()
title = 'Partial Beam Block ' + str(radars[nframe].fixed_angle['data']) + ' Degrees Elevation'
ax.set_axis_bgcolor('#E0E0E0')
display = pyart.graph.RadarMapDisplay(radars[nframe])
display.plot_ppi_map('partial_beam_block', resolution='l',
vmin=0, vmax=1, basemap=m, lat_lines=lats,
lon_lines=lons, cmap=plt.cm.PuRd, title=title)
display.plot_point(-28.0257, 39.0916,
symbol='ro')
plt.annotate('ARM ENA Site', m(-28.0217, 39.0916), color='k')
CS = m.contour(xx, yy, masked_data,
levels=np.arange(0, 400, 100),
linewidths=1, cmap=plt.cm.terrain)
plt.clabel(CS, inline=1, fontsize=10, fmt='%1.0f')
fig = plt.figure(figsize=(15, 7))
anim = animation.FuncAnimation(fig, animate, frames=len(radars))
anim.save('/home/zsherman/pbb_graciosa_animation.gif',
writer='imagemagick', fps=1)
plt.close()
In [579]:
def animate(nframe):
plt.clf()
title = 'Cumulative Beam Block ' + str(radars[nframe].fixed_angle['data']) + ' Degrees Elevation'
ax.set_axis_bgcolor('#E0E0E0')
display = pyart.graph.RadarMapDisplay(radars[nframe])
display.plot_ppi_map('cumulative_beam_block', resolution='l',
vmin=0, vmax=1, basemap=m, lat_lines=lats,
lon_lines=lons, cmap=plt.cm.PuRd, title=title)
display.plot_point(-28.0257, 39.0916,
symbol='ro')
plt.annotate('ARM ENA Site', m(-28.0217, 39.0916), color='k')
CS = m.contour(xx, yy, masked_data,
levels=np.arange(0, 400, 100),
linewidths=1, cmap=plt.cm.terrain)
plt.clabel(CS, inline=1, fontsize=10, fmt='%1.0f')
fig = plt.figure(figsize=(15, 7))
anim = animation.FuncAnimation(fig, animate, frames=len(radars))
anim.save('/home/zsherman/cbb_graciosa_animation.gif',
writer='imagemagick', fps=1)
plt.close()
In [ ]: