Using wradlib to obtain PBB and CBB, testing the terrain display of GeoTIFF file.
Code has been adapted by code from a notebook by Scott Collis: https://github.com/EVS-ATMOS/ARM-Azores/blob/master/notebooks/2.%20Terrain%20of%20ENA.ipynb
Code for projection conversion by Kai Muehlbauer and Wradlib notebook by Kai: http://wradlib.org/wradlib-docs/latest/notebooks/beamblockage/wradlib_beamblock.html
And Code by Nick Guy from PyRadarMet: https://github.com/nguy/PyRadarMet
In [1]:
import gdal
import osr
import numpy as np
import netCDF4
import re
import pyart
import matplotlib
import wradlib as wrl
import mayavi.mlab as mylab
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.basemap import Basemap
%matplotlib inline
In [2]:
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 [3]:
ds = gdal.Open("/home/zsherman/seattle.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 [4]:
print(xy_source.shape)
In [5]:
proj
Out[5]:
In [6]:
# Currently consumes too much memory.
"""
# Create the figure and basemap object
llcrnrlat=46.5
urcrnrlon=-120.5
llcrnrlon=-124.1
urcrnrlat=49.0
lats = np.arange(llcrnrlat, urcrnrlat, .5)
lons = np.arange(llcrnrlon, urcrnrlon, .5)
m = Basemap(
projection='lcc', lon_0=-122.4959, lat_0=48.1946,
resolution='l', llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat,
urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat)
# Create the projection objects for the convertion
# original (Albers)
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)
fig = plt.figure(figsize=(12, 6))
elev = data.T
the_sea = data.T < 0.1
masked_data = np.ma.masked_where(the_sea, elev)
# plot the data (first layer)
im1 = m.pcolormesh(xx, yy,
masked_data, cmap=plt.cm.terrain)
plt.colorbar(label = 'Elevation above MSL (m)')
# labels = [left,right,top,bottom]
m.drawparallels(lats, labels=[True,False,False,False])
m.drawmeridians(lons, labels=[False,False,False,True])
# annotate
m.drawcountries()
m.drawcoastlines(linewidth=.5)
plt.gca().set_axis_bgcolor('aqua')
"""
Out[6]:
In [7]:
radar = pyart.io.read('/home/zsherman/beam_block/data/seattle_radar.nc')
In [8]:
radar.altitude
Out[8]:
In [9]:
radar.fixed_angle
Out[9]:
In [10]:
sitecoords = (-122.4959, 48.1946, 195.0)
In [11]:
nrays = 360
nbins = 1832
el = 0.5
bw = 1
range_res = 250
In [12]:
r = np.linspace(2125, (nbins - 1)*range_res + 2125, nbins)
#r = np.arange(nbins) * range_res
beamradius = wrl.util.half_power_radius(r, bw)
In [13]:
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)
print("lon, lat, alt:", lon.shape, lat.shape, alt.shape)
In [14]:
rasterfile = '/home/zsherman/seattle.tif'
In [15]:
data_raster = wrl.io.open_raster(rasterfile)
In [16]:
proj_raster = wrl.georef.wkt_to_osr(data_raster.GetProjection())
In [17]:
x_pol, y_pol = wrl.georef.reproject(lon, lat, projection_target=proj_raster)
In [18]:
x_pol
Out[18]:
In [19]:
polcoords = np.dstack((x_pol, y_pol))
print("Polcoords:", polcoords.shape)
In [20]:
wrl.georef._get_range_resolution(r)
Out[20]:
In [21]:
rlimits = (x_pol.min(), y_pol.min(), x_pol.max(), y_pol.max())
print("Radar bounding box:\n\t%.2f\n%.2f %.2f\n\t%.2f" %
(y_pol.max(), x_pol.min(), x_pol.max(), y_pol.min()))
In [22]:
rastercoords, rastervalues = wrl.io.read_raster_data(rasterfile)
In [23]:
# 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]]
In [24]:
# Map rastervalues to polar grid points
polarvalues = wrl.ipol.cart2irregular_spline(
rastercoords, rastervalues, polcoords,
order=3, prefilter=False)
In [25]:
im2 = plt.pcolormesh(x_pol, y_pol,
polarvalues, cmap=plt.cm.terrain)
In [26]:
polarvalues.shape
Out[26]:
In [27]:
beam_height = wrl.qual.beam_height_ft(r, el)
In [28]:
elev = polarvalues
the_sea = polarvalues < 0.1
masked_pol_data = np.ma.masked_where(the_sea, elev)
In [29]:
PBB = wrl.qual.beam_block_frac(polarvalues, alt, beamradius)
PBB = np.ma.masked_invalid(PBB)
print(PBB.shape)
In [30]:
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 [31]:
# Also currently consumes too much memory.
"""
ds = gdal.Open('/home/zsherman/.tif')
data = ds.ReadAsArray()
the_sea = data < 0.1
masked_data = np.ma.masked_where(the_sea, data)
mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))
mlab.surf(masked_data, warp_scale=0.2)
mlap.points3d(x_pol, )
mlab.show()
"""
Out[31]:
In [32]:
fig = plt.figure(figsize=(10, 7))
PBB_sea = polarvalues < 0.1
PBB_masked_data = np.ma.masked_where(PBB_sea, PBB)
ax, dem = wrl.vis.plot_ppi(PBB_masked_data, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd)
ax.plot(0, 0, 'ro', )
ax.grid(True)
ax.annotate(' ARM ENA Site', (0, 0))
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')
ax.set_xlabel("Kilometers")
ax.set_ylabel("Kilometers")
ax.set_axis_bgcolor('#E0E0E0')
ax.set_xlim(-50000, 50000)
ax.set_ylim(-50000, 50000)
plt.colorbar(dem, ax=ax)
plt.savefig('/home/zsherman/beam_block/images/PBB_0_5deg.png', bbox_inches='tight')
plt.show()
In [33]:
fig = plt.figure(figsize=(10, 7))
#CBB_sea = polarvalues < 0.1
#CBB_masked_data = np.ma.masked_where(CBB_sea, CBB)
ax2, dem = wrl.vis.plot_ppi(CBB, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd, vmin=0, vmax=1)
#wrl.vis.plot_ppi_crosshair(dem, 10)
#ax2.set_xlim(-5000, 8000)
#ax2.set_ylim(-10000, 2000)
ax2.plot(0, 0, 'ro', )
ax2.grid(True)
ax2.annotate(' ARM ENA Site', (0, 0))
ticks = (ax2.get_xticks()/1000).astype(np.int)
ax2.set_xticklabels(ticks)
ticks = (ax2.get_yticks()/1000).astype(np.int)
ax2.set_yticklabels(ticks)
ax2.set_title('Cumulative Beam Block 0.5 Degrees')
ax2.set_xlabel("Kilometers")
ax2.set_ylabel("Kilometers")
ax2.set_axis_bgcolor('#E0E0E0')
plt.colorbar(dem, ax=ax2)
plt.savefig('/home/zsherman/beam_block/images/CBB_0_5deg.png', bbox_inches='tight')
plt.show()
In [34]:
fig = plt.figure(figsize=(10, 7))
CBB_sea = polarvalues < 0.1
CBB_masked_data = np.ma.masked_where(CBB_sea, CBB)
ax2, dem = wrl.vis.plot_ppi(CBB, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd, vmin=0, vmax=1)
ax2.set_xlim(-12000, 16000)
ax2.set_ylim(-8000, 11000)
ax2.plot(0, 0, 'ro', )
ax2.grid(True)
ax2.annotate(' ARM ENA Site', (0, 0))
ticks = (ax2.get_xticks()/1000).astype(np.int)
ax2.set_xticklabels(ticks)
ticks = (ax2.get_yticks()/1000).astype(np.int)
ax2.set_yticklabels(ticks)
ax2.set_title('Cumulative Beam Block 0.5 Degrees')
ax2.set_xlabel("Kilometers")
ax2.set_ylabel("Kilometers")
plt.colorbar(dem, ax=ax2)
plt.savefig('/home/zsherman/beam_block/images/CBB_zoomed_0_5deg.png', bbox_inches='tight')
plt.show()
In [35]:
# just a little helper function to style x and y axes of our maps
def annotate_map(ax, cm=None, title=""):
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_xlabel("Kilometers")
ax.set_ylabel("Kilometers")
if not cm is None:
plt.colorbar(cm, ax=ax)
if not title=="":
ax.set_title(title)
ax.grid()
polar_sea = polarvalues < 0.1
polar_masked_data = np.ma.masked_where(polar_sea, polarvalues)
fig = plt.figure(figsize=(10, 8))
# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)
# azimuth angle
angle = 50
# Plot terrain (on ax1)
ax1, dem = wrl.vis.plot_ppi(polarvalues,
ax=ax1, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.terrain, vmin=0.)
ax1.plot([0, np.sin(np.radians(angle))*1e5],
[0, np.cos(np.radians(angle))*1e5], "r-")
ax1.plot(0, 0, 'ro')
# ax1.annotate(' NEXRAD KATX Site',
#(sitecoords[0], sitecoords[1]))
ax1.set_xlim(-50000, 50000)
ax1.set_ylim(-50000, 50000)
annotate_map(ax1, dem, 'Terrain within 50 km range')
# Plot CBB (on ax2)
ax2, cbb = wrl.vis.plot_ppi(CBB, ax=ax2, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd, vmin=0, vmax=1)
ax2.plot(0, 0, 'ro', )
#ax2.annotate(' NEXRAD KATX Site', (sitecoords[0], sitecoords[1]))
ax2.set_xlim(-50000, 50000)
ax2.set_ylim(-50000, 50000)
annotate_map(ax2, cbb, 'CBB fraction 0.5 Degree Elevation')
# Plot single ray terrain profile on ax3
bc, = ax3.plot(r / 1000., alt[angle, :], '-b',
linewidth=3, label='Beam Center')
b3db, = ax3.plot(r / 1000., (alt[angle, :] + beamradius), ':b',
linewidth=1.5, label='3 dB Beam width')
ax3.plot(r / 1000., (alt[angle, :] - beamradius), ':b')
ax3.fill_between(r / 1000., 0.,
polarvalues[angle, :],
color='0.75')
ax3.set_xlim(0., np.max(r / 1000.) + 0.1)
ax3.set_ylim(0., 3000)
ax3.set_xlabel('Range (km)')
ax3.set_ylabel('Altitude (m)')
ax3.grid()
axb = ax3.twinx()
bbf, = axb.plot(r / 1000., CBB[angle, :], '-k',
label='BBF')
axb.set_ylabel('Beam-blockage fraction')
axb.set_ylim(0., 1.)
axb.set_xlim(0., 50)
legend = ax3.legend((bc, b3db, bbf),
('Beam Center', '3 dB Beam width', 'BBF'),
loc='upper left', fontsize=10)
plt.savefig('/home/zsherman/beam_block/images/3plot_seattle_50km_cbb.png',
bbox_inches='tight')
In [36]:
# just a little helper function to style x and y axes of our maps
def annotate_map(ax, cm=None, title=""):
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_xlabel("Kilometers")
ax.set_ylabel("Kilometers")
if not cm is None:
plt.colorbar(cm, ax=ax)
if not title=="":
ax.set_title(title)
ax.grid()
polar_sea = polarvalues < 0.1
polar_masked_data = np.ma.masked_where(polar_sea, polarvalues)
fig = plt.figure(figsize=(10, 8))
# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)
# azimuth angle
angle = 50
# Plot terrain (on ax1)
ax1, dem = wrl.vis.plot_ppi(polarvalues,
ax=ax1, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.terrain, vmin=0.)
ax1.plot([0,np.sin(np.radians(angle))*1e5],
[0,np.cos(np.radians(angle))*1e5],"r-")
ax1.plot(0, 0, 'ro')
ax1.set_xlim(-100000, 100000)
ax1.set_ylim(-100000, 100000)
annotate_map(ax1, dem, 'Terrain within 100 km range')
# Plot CBB (on ax2)
ax2, cbb = wrl.vis.plot_ppi(CBB, ax=ax2, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd, vmin=0, vmax=1)
ax2.plot(0, 0, 'ro', )
ax2.set_xlim(-100000, 100000)
ax2.set_ylim(-100000, 100000)
annotate_map(ax2, cbb, 'CBB fraction 0.5 Degree Elevation')
# Plot single ray terrain profile on ax3
bc, = ax3.plot(r / 1000., alt[angle, :], '-b',
linewidth=3, label='Beam Center')
b3db, = ax3.plot(r / 1000., (alt[angle, :] + beamradius), ':b',
linewidth=1.5, label='3 dB Beam width')
ax3.plot(r / 1000., (alt[angle, :] - beamradius), ':b')
ax3.fill_between(r / 1000., 0.,
polarvalues[angle, :],
color='0.75')
ax3.set_xlim(0., np.max(r / 1000.) + 0.1)
ax3.set_ylim(0., 3000)
ax3.set_xlabel('Range (km)')
ax3.set_ylabel('Altitude (m)')
ax3.grid()
axb = ax3.twinx()
bbf, = axb.plot(r / 1000., CBB[angle, :], '-k',
label='BBF')
axb.set_ylabel('Beam-blockage fraction')
axb.set_ylim(0., 1.)
axb.set_xlim(0., 100)
legend = ax3.legend((bc, b3db, bbf),
('Beam Center', '3 dB Beam width', 'BBF'),
loc='upper left', fontsize=10)
plt.savefig('/home/zsherman/beam_block/images/3plot_seattle_100km_cbb.png',
bbox_inches='tight')
In [37]:
fig = plt.figure(figsize=(10, 8))
# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)
# azimuth angle
angle = 50
# Plot terrain (on ax1)
ax1, dem = wrl.vis.plot_ppi(polarvalues,
ax=ax1, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.terrain, vmin=0.)
ax1.plot([0,np.sin(np.radians(angle))*1e5],
[0,np.cos(np.radians(angle))*1e5],"r-")
ax1.plot(0, 0, 'ro')
# ax1.annotate(' NEXRAD KATX Site',
#(sitecoords[0], sitecoords[1]))
ax1.set_xlim(-50000, 50000)
ax1.set_ylim(-50000, 50000)
annotate_map(ax1, dem, 'Terrain within 50 km range')
# Plot CBB (on ax2)
ax2, cbb = wrl.vis.plot_ppi(PBB, ax=ax2, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd, vmin=0, vmax=1)
ax2.plot(0, 0, 'ro')
#ax2.annotate(' NEXRAD KATX Site', (sitecoords[0], sitecoords[1]))
ax2.set_xlim(-50000, 50000)
ax2.set_ylim(-50000, 50000)
annotate_map(ax2, cbb, 'PBB fraction 0.5 Degree Elevation')
# Plot single ray terrain profile on ax3
bc, = ax3.plot(r / 1000., alt[angle, :], '-b',
linewidth=3, label='Beam Center')
b3db, = ax3.plot(r / 1000., (alt[angle, :] + beamradius), ':b',
linewidth=1.5, label='3 dB Beam width')
ax3.plot(r / 1000., (alt[angle, :] - beamradius), ':b')
ax3.fill_between(r / 1000., 0.,
polarvalues[angle, :],
color='0.75')
ax3.set_xlim(0., np.max(r / 1000.) + 0.1)
ax3.set_ylim(0., 3000)
ax3.set_xlabel('Range (km)')
ax3.set_ylabel('Altitude (m)')
ax3.grid()
axb = ax3.twinx()
bbf, = axb.plot(r / 1000., CBB[angle, :], '-k',
label='BBF')
axb.set_ylabel('Beam-blockage fraction')
axb.set_ylim(0., 1.)
axb.set_xlim(0., 50)
legend = ax3.legend((bc, b3db, bbf),
('Beam Center', '3 dB Beam width', 'BBF'),
loc='upper left', fontsize=10)
plt.savefig('/home/zsherman/beam_block/images/3plot_seattle_50km_pbb.png',
bbox_inches='tight')
In [38]:
fig = plt.figure(figsize=(10, 8))
# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)
# azimuth angle
angle = 50
# Plot terrain (on ax1)
ax1, dem = wrl.vis.plot_ppi(polarvalues,
ax=ax1, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.terrain, vmin=0.)
ax1.plot([0,np.sin(np.radians(angle))*1e5],
[0,np.cos(np.radians(angle))*1e5],"r-")
ax1.plot(0, 0, 'ro')
ax1.set_xlim(-100000, 100000)
ax1.set_ylim(-100000, 100000)
annotate_map(ax1, dem, 'Terrain within 100 km range')
# Plot CBB (on ax2)
ax2, cbb = wrl.vis.plot_ppi(PBB, ax=ax2, r=r,
az=np.degrees(coord[:, 0, 1]),
cmap=plt.cm.PuRd, vmin=0, vmax=1)
ax2.plot(0, 0, 'ro', )
ax2.set_xlim(-100000, 100000)
ax2.set_ylim(-100000, 100000)
annotate_map(ax2, cbb, 'PBB fraction 0.5 Degree Elevation')
# Plot single ray terrain profile on ax3
bc, = ax3.plot(r / 1000., alt[angle, :], '-b',
linewidth=3, label='Beam Center')
b3db, = ax3.plot(r / 1000., (alt[angle, :] + beamradius), ':b',
linewidth=1.5, label='3 dB Beam width')
ax3.plot(r / 1000., (alt[angle, :] - beamradius), ':b')
ax3.fill_between(r / 1000., 0.,
polarvalues[angle, :],
color='0.75')
ax3.set_xlim(0., np.max(r / 1000.) + 0.1)
ax3.set_ylim(0., 3000)
ax3.set_xlabel('Range (km)')
ax3.set_ylabel('Altitude (m)')
ax3.grid()
axb = ax3.twinx()
bbf, = axb.plot(r / 1000., CBB[angle, :], '-k',
label='BBF')
axb.set_ylabel('Beam-blockage fraction')
axb.set_ylim(0., 1.)
axb.set_xlim(0., 100)
legend = ax3.legend((bc, b3db, bbf),
('Beam Center', '3 dB Beam width', 'BBF'),
loc='upper left', fontsize=10)
plt.savefig('/home/zsherman/beam_block/images/3plot_seattle_100km_pbb.png',
bbox_inches='tight')