import pyart from matplotlib import pyplot as plt import numpy as np from netCDF4 import num2date from copy import deepcopy %matplotlib inline print pyart.__version__ excluded_fields = ['velocity', 'spectrum_width', 'differential_phase'] keep_nsweeps = 4 joliet_radar = pyart.io.read('data/KLOT_110622.nexrad', exclude_fields=excluded_fields) joliet_radar = joliet_radar.extract_sweeps(range(keep_nsweeps)) lincoln_radar = pyart.io.read('data/KILX_110437.nexrad', exclude_fields=excluded_fields) lincoln_radar = lincoln_radar.extract_sweeps(range(keep_nsweeps)) davenport_radar = pyart.io.read('data/KDVN_110608.nexrad', exclude_fields=excluded_fields) davenport_radar = davenport_radar.extract_sweeps(range(keep_nsweeps)) for radar in [joliet_radar, lincoln_radar, davenport_radar]: print radar.fields.keys() #set the domain: max_lat = 43. min_lat = 40 min_lon = -90.5 max_lon = -86 # create an instance of the display class display = pyart.graph.RadarMapDisplay(joliet_radar) # generate the figure plt.figure(figsize = [17,4]) # first subplot: Reflectivity plt.subplot(1, 3, 1) display.plot_ppi_map('reflectivity', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon=max_lon, vmin=-8, vmax=64, lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') # second subplot: differential reflectivity plt.subplot(1, 3, 2) display.plot_ppi_map('differential_reflectivity', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon= max_lon, vmin=-2, vmax=6, lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') #third subplot: correlation coef plt.subplot(1, 3, 3) display.plot_ppi_map('cross_correlation_ratio', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon=max_lon, vmin=.5, vmax=1., lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') del display # create an instance of the display class display = pyart.graph.RadarMapDisplay(lincoln_radar) # generate the figure plt.figure(figsize=[17,4]) # first subplot: Reflectivity plt.subplot(1, 3, 1) display.plot_ppi_map('reflectivity', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon=max_lon, vmin=-8, vmax=64, lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') # second subplot: differential reflectivity plt.subplot(1, 3, 2) display.plot_ppi_map('differential_reflectivity', max_lat= max_lat, min_lat=min_lat, min_lon= min_lon, max_lon= max_lon, vmin=-2, vmax=6, lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') # third subplot: correlation coef plt.subplot(1, 3, 3) display.plot_ppi_map('cross_correlation_ratio', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon= max_lon, vmin=.5, vmax=1., lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') del display # create an instance of the display class display = pyart.graph.RadarMapDisplay(davenport_radar) # generate the figure plt.figure(figsize=[17,4]) # first subplot: Reflectivity plt.subplot(1, 3, 1) display.plot_ppi_map('reflectivity', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon= max_lon, vmin=-8, vmax=64, lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') # second subplot: differential reflectivity plt.subplot(1, 3, 2) display.plot_ppi_map('differential_reflectivity', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon= max_lon, vmin=-2, vmax=6, lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') # third subplot: correlation coef plt.subplot(1, 3, 3) display.plot_ppi_map('cross_correlation_ratio', max_lat=max_lat, min_lat=min_lat, min_lon=min_lon, max_lon=max_lon, vmin=.5, vmax=1., lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') del display def smooth_diff_refl(radar): smooth_zdr = radar.fields['differential_reflectivity']['data'].copy() smooth_zdr[np.where(radar.fields['cross_correlation_ratio']['data']<0.5)] = np.ma.masked for i in range(smooth_zdr.shape[0]): smooth_zdr[i,:] = pyart.correct.phase_proc.smooth_and_trim( smooth_zdr[i,:], 8) radar.add_field_like('differential_reflectivity', 'differential_reflectivity_smooth', smooth_zdr) smooth_diff_refl(joliet_radar) smooth_diff_refl(davenport_radar) smooth_diff_refl(lincoln_radar) mesh_mapped_ALL = pyart.map.grid_from_radars( (joliet_radar, davenport_radar, lincoln_radar), grid_shape=(1, 301, 301), grid_limits=((1000, 1000), (-300.*1000., 300.*1000.), (-300.*1000., 300.*1000.)), grid_origin = (41.5, -88.5), fields=['reflectivity', 'differential_reflectivity', 'cross_correlation_ratio', 'differential_reflectivity_smooth'], refl_field='reflectivity', max_refl=100., copy_field_data=False) print mesh_mapped_ALL.fields.keys() print mesh_mapped_ALL.axes.keys() plt.figure(figsize = [15,10]) plt.pcolormesh( mesh_mapped_ALL.axes['x_disp']['data'], mesh_mapped_ALL.axes['y_disp']['data'], mesh_mapped_ALL.fields['reflectivity']['data'][0,:,:], vmin=-8, vmax=64) plt.colorbar() plt.figure(figsize = [15,10]) plt.pcolormesh( mesh_mapped_ALL.axes['x_disp']['data'], mesh_mapped_ALL.axes['y_disp']['data'], mesh_mapped_ALL.fields['ROI']['data'][0,:,:]) plt.colorbar() max_lat = 43. min_lat = 40 min_lon = -90.5 max_lon = -86 plt.figure(figsize = [15,10]) display = pyart.graph.GridMapDisplay(mesh_mapped_ALL) display.plot_basemap(lat_lines=np.arange(min_lat,max_lat,.5), lon_lines=np.arange(min_lon, max_lon, 1), resolution='i') display.plot_grid('reflectivity', vmin=-8, vmax=64) display.plot_colorbar() del display mesh_mapped_ALL_3d = pyart.io.read_grid('data/IL_grid.nc') display = pyart.graph.GridMapDisplay(mesh_mapped_ALL_3d) fig = plt.figure(figsize=[15, 8]) map_panel_axes = [0.05, 0.05, .4, .80] x_cut_panel_axes = [0.55, 0.10, .4, .30] y_cut_panel_axes = [0.55, 0.50, .4, .30] colorbar_panel_axes = [0.05, 0.90, .4, .03] # parameters level = 2 vmin = -8 vmax = 64 lat = 41.7092 lon = -87.9820 # panel 1, basemap and radar reflectivity ax1 = fig.add_axes(map_panel_axes) display.plot_basemap(resolution='i') display.plot_grid('reflectivity', level=level, vmin=vmin, vmax=vmax, cmap=pyart.graph.cm.NWSRef) display.plot_crosshairs(lon=lon, lat=lat) cbax = fig.add_axes(colorbar_panel_axes) display.plot_colorbar(label='Eq. Reflectivity Factor (dBZ)', cax = cbax ) # panel 2, longitude slice. ax2 = fig.add_axes(x_cut_panel_axes) display.plot_longitude_slice('reflectivity', lon=lon, lat=lat, cmap=pyart.graph.cm.NWSRef, vmin=vmin, vmax=vmax) ax2.set_xlabel('Distance from Argonne (km)') # panel 3, latitude slice ax3 = fig.add_axes(y_cut_panel_axes) display.plot_latitude_slice('reflectivity', lon=lon, lat=lat, cmap=pyart.graph.cm.NWSRef, vmin=vmin, vmax=vmax) # add a title slc_height = mesh_mapped_ALL_3d.axes['z_disp']['data'][level] dts = num2date(mesh_mapped_ALL_3d.axes['time']['data'], mesh_mapped_ALL_3d.axes['time']['units']) datestr = dts[0].strftime('%H:%M UTC on %d %b %Y') title = 'Sliced at ' + str(slc_height) + ' meters agl at ' + datestr fig.text(0.5, 0.9, title) fig = plt.figure(figsize=[15, 8]) # parameters vmin = -1 vmax = 6 # panel 1, basemap and radar reflectivity ax1 = fig.add_axes(map_panel_axes) display.plot_basemap(resolution='i') display.plot_grid('differential_reflectivity', level=level, vmin=vmin, vmax=vmax) display.plot_crosshairs(lon=lon, lat=lat) cbax = fig.add_axes(colorbar_panel_axes) display.plot_colorbar(label='Eq. Reflectivity Factor (dBZ)', cax = cbax ) # panel 2, longitude slice. ax2 = fig.add_axes(x_cut_panel_axes) display.plot_longitude_slice('differential_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax) ax2.set_xlabel('Distance from Argonne (km)') # panel 3, latitude slice ax3 = fig.add_axes(y_cut_panel_axes) display.plot_latitude_slice('differential_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax) # add a title slc_height = mesh_mapped_ALL_3d.axes['z_disp']['data'][level] dts = num2date(mesh_mapped_ALL_3d.axes['time']['data'], mesh_mapped_ALL_3d.axes['time']['units']) datestr = dts[0].strftime('%H:%M UTC on %d %b %Y') title = 'Sliced at ' + str(slc_height) + ' meters agl at ' + datestr fig.text(0.5, 0.9, title) fig = plt.figure(figsize=[15, 8]) # parameters vmin = -1 vmax = 6 # panel 1, basemap and radar reflectivity ax1 = fig.add_axes(map_panel_axes) display.plot_basemap(resolution='i') display.plot_grid('differential_reflectivity_smooth', level=level, vmin=vmin, vmax=vmax) display.plot_crosshairs(lon=lon, lat=lat) cbax = fig.add_axes(colorbar_panel_axes) display.plot_colorbar(label='Eq. Reflectivity Factor (dBZ)', cax = cbax ) # panel 2, longitude slice. ax2 = fig.add_axes(x_cut_panel_axes) display.plot_longitude_slice('differential_reflectivity_smooth', lon=lon, lat=lat, vmin=vmin, vmax=vmax) ax2.set_xlabel('Distance from Argonne (km)') # panel 3, latitude slice ax3 = fig.add_axes(y_cut_panel_axes) display.plot_latitude_slice('differential_reflectivity_smooth', lon=lon, lat=lat, vmin=vmin, vmax=vmax) # add a title slc_height = mesh_mapped_ALL_3d.axes['z_disp']['data'][level] dts = num2date(mesh_mapped_ALL_3d.axes['time']['data'], mesh_mapped_ALL_3d.axes['time']['units']) datestr = dts[0].strftime('%H:%M UTC on %d %b %Y') title = 'Sliced at ' + str(slc_height) + ' meters agl at ' + datestr fig.text(0.5, 0.9, title)