.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_radar_ppi.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_radar_ppi.py: Radar PPI ============================== Display radar Plan Position Indicators (PPIs) using a combination of *radar_tools*, *geo_tools* and *legs*. - Data consists of a radar Volume scan in the ODIM HDF5 format that are read with *domutils.radar_tools*. - Data is projected on Cartopy geoaxes using *domutils.geo_tools* - Various color mappings are applied with *domutils.legs* Some advanced concepts are demonstrated in this example. - Quality controls are applied to Doppler velocity based on the Depolarization Ratio (DR) and the total quality index. - Matplotlib/cartopy geoaxes are clipped and superposed to show a closeup of a possible meso-cyclone. .. GENERATED FROM PYTHON SOURCE LINES 20-457 .. code-block:: Python def add_feature(ax): """plot geographical and political boundaries in matplotlib axes """ import cartopy.feature as cfeature ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=0.1, edgecolor='0.1',zorder=1) def radar_ax_circ(ax, radar_lat, radar_lon): '''plot azimuth lines and range circles around a given radar ''' import numpy as np import cartopy.crs as ccrs import domutils.geo_tools as geo_tools #cartopy transform for latlons proj_pc = ccrs.PlateCarree() color=(100./256.,100./256.,100./256.) #add a few azimuths lines az_lines = np.arange(0,360.,90.) ranges = np.arange(250.) for this_azimuth in az_lines: lons, lats = geo_tools.lat_lon_range_az(radar_lon, radar_lat, ranges, this_azimuth) ax.plot(lons, lats, transform=proj_pc, c=color, zorder=300, linewidth=.3) #add a few range circles ranges = np.arange(0,250.,100.) azimuths = np.arange(0,361.)#360 degree will be included for full circles for this_range in ranges: lons, lats = geo_tools.lat_lon_range_az(radar_lon, radar_lat, this_range, azimuths) ax.plot(lons, lats, transform=proj_pc, c=color, zorder=300, linewidth=.3) def main(): import os, inspect import copy import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import cartopy.crs as ccrs # imports from domutils import domutils.legs as legs import domutils.geo_tools as geo_tools import domutils.radar_tools as radar_tools #missing values ane the associated color missing = -9999. missing_color = 'grey_160' undetect = -3333. undetect_color = 'grey_180' #file to read currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) parentdir = os.path.dirname(currentdir) #directory where this package lives odim_file = parentdir + '/test_data/odimh5_radar_volume_scans/2019071120_24_ODIMH5_PVOL6S_VOL.qc.casbv.h5' #Read a PPI from a volume scan in the ODIM H5 format #we want the PPI with a nominal elevation of 0.4 degree #we retrieve reflectivity (dbzh) and Doppler velocity (vradh) and the Depolarization ratio (dr) #the reader computes the lat/lon of data points in the PPI for you out_dict = radar_tools.read_h5_vol(odim_file=odim_file, elevations=[0.4], quantities=['dbzh','vradh','dr'], no_data=missing, undetect=undetect, include_quality=True, latlon=True) #radar properties radar_lat = out_dict['radar_lat'] radar_lon = out_dict['radar_lon'] #PPIs # You can see the quantities available in the dict with # print(out_dict['0.4'].keys()) reflectivity = out_dict['0.4']['dbzh'] doppvelocity = out_dict['0.4']['vradh'] dr = out_dict['0.4']['dr'] tot_qi = out_dict['0.4']['quality_qi_total'] latitudes = out_dict['0.4']['latitudes'] longitudes = out_dict['0.4']['longitudes'] # #Quality controls on Doppler velocity #pixel is considered weather is DR < 12 dB weather_target = np.where( dr < -12., 1, 0) #not a weather target when DR is missing weather_target = np.where( np.isclose(dr, missing), 0, weather_target) #not a weather target when DR is undetect weather_target = np.where( np.isclose(dr, undetect), 0, weather_target) #a 3D representation of a boxcar filter see radar_tools API for docs on this function #5x5 search area in polar coordinates remapped_data_5 = radar_tools.median_filter.remap_data(weather_target, window=5, mode='wrap') #if less than 12 points (half the points in a 5x5 window) have good dr, pixel is marked as clutter is_clutter = np.where(np.sum(remapped_data_5, axis=2) <= 12, True, False) #3x3 search area in polar coordinates remapped_data_3 = radar_tools.median_filter.remap_data(weather_target, window=3, mode='wrap') #if less than 9 points (all the points in a 3x3 window) have good dr, pixel is marked as clutter is_clutter = np.where(np.sum(remapped_data_3, axis=2) < 9, True, is_clutter) # doppvelocity_qc = copy.deepcopy(doppvelocity) #DR filtering only applied to points with reflectivities < 35 dB doppvelocity_qc = np.where( np.logical_and(is_clutter , reflectivity < 35.) , undetect, doppvelocity_qc) #add filtering using total quality index doppvelocity_qc = np.where( tot_qi < 0.2 , undetect, doppvelocity_qc) #pixel density of panels in figure ratio = 1. hpix = 800. #number of horizontal pixels E-W vpix = ratio*hpix #number of vertical pixels S-N img_res = (int(hpix),int(vpix)) #cartopy Rotated Pole projection pole_latitude=90. pole_longitude=0. proj_rp = ccrs.RotatedPole(pole_latitude=pole_latitude, pole_longitude=pole_longitude) #plate carree proj_pc = ccrs.PlateCarree() #a domain ~250km around the Montreal area where the radar is lat_0 = 45.7063 delta_lat = 2.18 lon_0 = -73.85 delta_lon = 3.12 map_extent=[lon_0-delta_lon, lon_0+delta_lon, lat_0-delta_lat, lat_0+delta_lat] #a smaller domain for a closeup that will be inlaid in figure lat_0 = 46.4329666 delta_lat = 0.072666 lon_0 = -75.00555 delta_lon = 0.104 map_extent_close=[lon_0-delta_lon, lon_0+delta_lon, lat_0-delta_lat, lat_0+delta_lat] #a circular clipping mask for the closeup axes x = np.linspace(-1., 1, int(hpix)) y = np.linspace(-1., 1, int(vpix)) xx, yy = np.meshgrid(x, y, indexing='ij') clip_alpha = np.where( np.sqrt(xx**2.+yy**2.) < 1., 1., 0. ) #a circular line around the center of the closeup window radius=8. #km azimuths = np.arange(0,361.)#360 degree will be included for full circles closeup_lons, closeup_lats = geo_tools.lat_lon_range_az(lon_0, lat_0, radius, azimuths) #a line 5km long for showing scale in closeup insert azimuth = 90 #meteorological angle distance = np.linspace(0,5.,50)#360 degree will be included for full circles scale_lons, scale_lats = geo_tools.lat_lon_range_az(lon_0-0.07, lat_0-0.04, distance, azimuth) #point density for figure mpl.rcParams['figure.dpi'] = 100. #crank this up for high def images mpl.rcParams.update({'font.size': 25}) mpl.rcParams.update({'font.family':'Latin Modern Roman'}) # dimensions for figure panels and spaces # all sizes are inches for consistency with matplotlib fig_w = 13.5 # Width of figure fig_h = 12.5 # Height of figure rec_w = 5. # Horizontal size of a panel rec_h = ratio * rec_w # Vertical size of a panel sp_w = .5 # horizontal space between panels sp_h = .8 # vertical space between panels fig = plt.figure(figsize=(fig_w,fig_h)) xp = .02 #coords of title (axes normalized coordinates) yp = 1.02 #coords for the closeup that is overlayed x0 = .55 #x-coord of bottom left position of closeup (axes coords) y0 = .55 #y-coord of bottom left position of closeup (axes coords) dx = .4 #x-size of closeup axes (fraction of a "regular" panel) dy = .4 #y-size of closeup axes (fraction of a "regular" panel) #normalize sizes to obtain figure coordinates (0-1 both horizontally and vertically) rec_w = rec_w / fig_w rec_h = rec_h / fig_h sp_w = sp_w / fig_w sp_h = sp_h / fig_h #instantiate objects to handle geographical projection of data proj_inds = geo_tools.ProjInds(src_lon=longitudes, src_lat=latitudes, extent=map_extent, dest_crs=proj_rp, extend_x=False, extend_y=True, image_res=img_res, missing=missing) #for closeup image proj_inds_close = geo_tools.ProjInds(src_lon=longitudes, src_lat=latitudes, extent=map_extent_close, dest_crs=proj_rp, extend_x=False, extend_y=True, image_res=img_res, missing=missing) # #Reflectivity # #axes for this plot pos = [sp_w, sp_h+(sp_h+rec_h), rec_w, rec_h] ax = fig.add_axes(pos, projection=proj_rp) ax.spines['geo'].set_linewidth(0.3) ax.set_extent(proj_inds.rotated_extent, crs=proj_rp) ax.set_aspect('auto') ax.annotate('Qced Reflectivity', size=30, xy=(xp, yp), xycoords='axes fraction') # colormapping object for reflectivity map_reflectivity = legs.PalObj(range_arr=[0.,60], n_col=6, over_high='extend', under_low='white', excep_val=[missing, undetect], excep_col=[missing_color, undetect_color]) #geographical projection of data into axes space proj_data = proj_inds.project_data(reflectivity) #plot data & palette map_reflectivity.plot_data(ax=ax, data=proj_data, zorder=0, palette='right', pal_units='[dBZ]', pal_format='{:3.0f}') #add political boundaries add_feature(ax) #radar circles and azimuths radar_ax_circ(ax, radar_lat, radar_lon) #circle indicating closeup area ax.plot(closeup_lons, closeup_lats, transform=proj_pc, c=(0.,0.,0.), zorder=300, linewidth=.8) #arrow pointing to closeup ax.annotate("", xy=(0.33, 0.67), xytext=(.55, .74), xycoords='axes fraction', arrowprops=dict(arrowstyle="<-")) # #Closeup of reflectivity # pos = [sp_w+x0*rec_w, sp_h+(sp_h+rec_h)+y0*rec_h, dx*rec_w, dy*rec_h] ax2 = fig.add_axes(pos, projection=proj_rp, label='reflectivity overlay') ax2.set_extent(proj_inds_close.rotated_extent, crs=proj_rp) ax2.spines['geo'].set_linewidth(0.0) #no border line ax2.set_facecolor((1.,1.,1.,0.)) #transparent background #geographical projection of data into axes space proj_data = proj_inds_close.project_data(reflectivity) #RGB values for data to plot closeup_rgb = map_reflectivity.to_rgb(proj_data) #get corners of image in data space extent_data_space = ax2.get_extent() ## another way of doing the same thing is to get an object that convers axes coords to data coords ## this method is more powerfull as it will return data coords of any points in axes space #transform_data_to_axes = ax2.transData + ax2.transAxes.inverted() #transform_axes_to_data = transform_data_to_axes.inverted() #pts = ((0.,0.),(1.,1.)) #axes space coords #pt1, pt2 = transform_axes_to_data.transform(pts) #extent_data_space = [pt1[0],pt2[0],pt1[1],pt2[1]] #add alpha channel (transparency) to closeup image rgba = np.concatenate([closeup_rgb/255.,clip_alpha[...,np.newaxis]], axis=2) #plot image ax2.imshow(rgba, interpolation='nearest', origin='upper', extent=extent_data_space, zorder=100) ax2.set_aspect('auto') #circle indicating closeup area circle = ax2.plot(closeup_lons, closeup_lats, transform=proj_pc, c=(0.,0.,0.), zorder=300, linewidth=1.5) #prevent clipping of the circle we just drawn for line in ax2.lines: line.set_clip_on(False) # #Quality Controlled Doppler velocity # #axes for this plot pos = [sp_w+(sp_w+rec_w+1./fig_w), sp_h+(sp_h+rec_h), rec_w, rec_h] ax = fig.add_axes(pos, projection=proj_rp) ax.set_extent(proj_inds.rotated_extent, crs=proj_rp) ax.set_aspect('auto') ax.annotate('Qced Doppler velocity', size=30, xy=(xp, yp), xycoords='axes fraction') #from https://colorbrewer2.org brown_purple=[[ 45, 0, 75], [ 84, 39,136], [128,115,172], [178,171,210], [216,218,235], [247,247,247], [254,224,182], [253,184, 99], [224,130, 20], [179, 88, 6], [127, 59, 8]] range_arr = [-48.,-40.,-30.,-20,-10.,-1.,1.,10.,20.,30.,40.,48.] map_dvel = legs.PalObj(range_arr=range_arr, color_arr = brown_purple, solid='supplied', excep_val=[missing, undetect], excep_col=[missing_color, undetect_color]) #geographical projection of data into axes space proj_data = proj_inds.project_data(doppvelocity_qc) #plot data & palette map_dvel.plot_data(ax=ax,data=proj_data, zorder=0, palette='right', pal_units='[m/s]', pal_format='{:3.0f}') #palette options #add political boundaries add_feature(ax) #radar circles and azimuths radar_ax_circ(ax, radar_lat, radar_lon) #circle indicating closeup area ax.plot(closeup_lons, closeup_lats, transform=proj_pc, c=(0.,0.,0.), zorder=300, linewidth=.8) #arrow pointing to closeup ax.annotate("", xy=(0.33, 0.67), xytext=(.55, .74), xycoords='axes fraction', arrowprops=dict(arrowstyle="<-")) # #Closeup of Doppler velocity # pos = [sp_w+1.*(sp_w+rec_w+1./fig_w)+x0*rec_w, sp_h+(sp_h+rec_h)+y0*rec_h, dx*rec_w, dy*rec_h] ax2 = fig.add_axes(pos, projection=proj_rp, label='overlay') ax2.set_extent(proj_inds_close.rotated_extent, crs=proj_rp) ax2.spines['geo'].set_linewidth(0.0) #no border line ax2.set_facecolor((1.,1.,1.,0.)) #transparent background #geographical projection of data into axes space proj_data = proj_inds_close.project_data(doppvelocity_qc) #RGB values for data to plot closeup_rgb = map_dvel.to_rgb(proj_data) #get corners of image in data space extent_data_space = ax2.get_extent() #add alpha channel (transparency) to closeup image rgba = np.concatenate([closeup_rgb/255.,clip_alpha[...,np.newaxis]], axis=2) #plot image ax2.imshow(rgba, interpolation='nearest', origin='upper', extent=extent_data_space, zorder=100) ax2.set_aspect('auto') #line indicating closeup area circle = ax2.plot(closeup_lons, closeup_lats, transform=proj_pc, c=(0.,0.,0.), zorder=300, linewidth=1.5) for line in ax2.lines: line.set_clip_on(False) #Show scale in inlay ax2.plot(scale_lons, scale_lats, transform=proj_pc, c=(0.,0.,0.), zorder=300, linewidth=.8) ax2.annotate("5 km", size=18, xy=(.16, .25), xycoords='axes fraction', zorder=310) # #DR # #axes for this plot pos = [sp_w, sp_h, rec_w, rec_h] ax = fig.add_axes(pos, projection=proj_rp) ax.set_extent(proj_inds.rotated_extent, crs=proj_rp) ax.set_aspect('auto') ax.annotate('Depolarization ratio', size=30, xy=(xp, yp), xycoords='axes fraction') # Set up colormapping object map_dr = legs.PalObj(range_arr=[-36.,-24.,-12., 0.], color_arr=['purple','blue','orange'], dark_pos =['high', 'high','low'], excep_val=[missing, undetect], excep_col=[missing_color, undetect_color]) #geographical projection of data into axes space proj_data = proj_inds.project_data(dr) #plot data & palette map_dr.plot_data(ax=ax,data=proj_data, zorder=0, palette='right', pal_units='[dB]', pal_format='{:3.0f}') #add political boundaries add_feature(ax) #radar circles and azimuths radar_ax_circ(ax, radar_lat, radar_lon) # #Total quality index # #axes for this plot pos = [sp_w+(sp_w+rec_w+1./fig_w), sp_h, rec_w, rec_h] ax = fig.add_axes(pos, projection=proj_rp) ax.set_extent(proj_inds.rotated_extent, crs=proj_rp) ax.set_aspect('auto') ax.annotate('Total quality index', size=30, xy=(xp, yp), xycoords='axes fraction') # Set up colormapping object pastel = [ [[255,190,187],[230,104, 96]], #pale/dark red [[255,185,255],[147, 78,172]], #pale/dark purple [[255,227,215],[205,144, 73]], #pale/dark brown [[210,235,255],[ 58,134,237]], #pale/dark blue [[223,255,232],[ 61,189, 63]] ] #pale/dark green map_qi = legs.PalObj(range_arr=[0., 1.], dark_pos='high', color_arr=pastel, excep_val=[missing, undetect], excep_col=[missing_color, undetect_color]) #geographical projection of data into axes space proj_data = proj_inds.project_data(tot_qi) #plot data & palette map_qi.plot_data(ax=ax,data=proj_data, zorder=0, palette='right', pal_units='[unitless]', pal_format='{:3.1f}') #palette options #add political boundaries add_feature(ax) #radar circles and azimuths radar_ax_circ(ax, radar_lat, radar_lon) ##uncomment to save figure #plt.savefig('radar_ppi.svg') if __name__ == '__main__': main() .. _sphx_glr_download_auto_examples_plot_radar_ppi.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_radar_ppi.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_radar_ppi.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_radar_ppi.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_