ppi_Pseudo_RHI_azimuth.ipynb

Kernel: Python 3 (ipykernel)

PPI, Pseudo RHI, Py-ART

• @author: Hamid Ali Syed

• @date: Dec 10, 2022

import warnings
warnings.filterwarnings("ignore")
import pyart
import numpy as np
from tqdm.notebook import trange, tqdm
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as feat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

def add_map(ax, b = 0, t=0, l = 0, r = 0):
    gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.3, color='black', alpha=0.3, 
                         linestyle='-', draw_labels=True)
    gl.xlabels_top = t
    gl.xlabels_bottom = b
    gl.ylabels_left = l
    gl.ylabels_right=r
    gl.xlines = True
    gl.ylines = True
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    ax.add_feature(feat.BORDERS, lw = 0.5)
    ax.add_feature(feat.STATES.with_scale("10m"), alpha = 0.5, lw = 0.4, ls=":")

radar = pyart.io.read_cfradial("polar_GOA210516024101.nc")

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return array[idx]

projection=ccrs.LambertAzimuthalEqualArea(radar.longitude['data'][0],radar.latitude['data'][0])
pdisplay = pyart.graph.RadarMapDisplay(radar)
# azi = 270
for i in trange(0, 360, 5):
    azi = find_nearest(radar.azimuth['data'][radar.get_slice(0)], i)
#     print(f"Plotting azimuth {azi}")
    pseudorhi = pyart.util.cross_section_ppi(radar, [azi])
    display = pyart.graph.RadarDisplay(pseudorhi)
    fig = plt.figure(figsize=[14, 6])

    map_panel_axes = [0.05, 0.05, .4, .80]
    x_cut_panel_axes = [0.55, 0.10, .4, .25]
    y_cut_panel_axes = [0.55, 0.50, .4, .25]

    # Panel 1: PPI plot of the second tilt.
    ax1 = fig.add_axes(map_panel_axes, projection=projection)
    pdisplay.plot_ppi_map("REF", cmap = "pyart_NWSRef", projection=projection, resolution="50m", fig=fig, ax=ax1)
    pdisplay.plot_range_rings(np.linspace(50, 250, 3), lw = 0.5, ax=ax1)

    radar_lat = radar.latitude['data']
    radar_lon = radar.longitude['data']

    # Create azimuth lines
    dtor = np.pi/180.0
    max_range = 250
    maxrange_meters = max_range * 1000.
    meters_to_lat = 1. / 111177.
    meters_to_lon =  1. / (111177. * np.cos(radar_lat * dtor))

    # for azi in range(0,360,45): # 45 degree intervals
    azimuth = 90. - azi
    dazimuth = azimuth * dtor
    lon_maxrange = radar_lon + np.cos(dazimuth) * meters_to_lon * maxrange_meters
    lat_maxrange = radar_lat + np.sin(dazimuth) * meters_to_lat * maxrange_meters
    pdisplay.plot_line_geo([radar_lon, lon_maxrange], [radar_lat, lat_maxrange], line_style='-', 
                           lw=2, color='yellow', transform=ccrs.PlateCarree())
    add_map(ax=ax1, b=1, l=1)
    # Panel 2: RHI slice
    ax2 = fig.add_axes(x_cut_panel_axes)
    display.plot_rhi('REF', ax = ax2, cmap = "pyart_NWSRef")
    ax2.set_ylim([0,20])
    ax2.set_xlim([0,200])
    # Panel 3: RHI slice
    ax3 = fig.add_axes(y_cut_panel_axes)
    display.plot_rhi('VELH', ax = ax3, vmin = -30, vmax = 30, cmap = "pyart_NWSVel")
    ax3.set_ylim([0,20])
    ax3.set_xlim([0,200])
    plt.savefig(f"azimuth_ppi_rhi/{i}_RHI_Tauktae_azimuth_{azi}"+".png", bbox_inches="tight")
    plt.close()

  0%|          | 0/72 [00:00<?, ?it/s]

import re
def natural_sort_key(s, _re=re.compile(r'(\d+)')):
    return [int(t) if i & 1 else t.lower() for i, t in enumerate(_re.split(s))]

import glob
from PIL import Image

# filepaths
fp_in = "azimuth_ppi_rhi/*png"
fp_out = 'project.gif'
files = sorted(glob.glob(fp_in))
sorted_files = sorted(files, key=natural_sort_key)
img, *imgs = [Image.open(f) for f in sorted_files]
img.save(fp=fp_out, format='GIF', append_images=imgs,
         save_all=True, duration=300, loop=0)

from IPython.display import Image
Image(url='project.gif')