<jupyter_start><jupyter_text>Dual Doppler lobe plotterDD lobe plotter. Needs py-ART Grid file since DD lobes are calculated in radar relative coordinates.Based on code created by <NAME>, <NAME>, <NAME>, and myself<jupyter_code>import pyart
import gzip
from matplotlib import pyplot as plt
from matplotlib import rcParams
from scipy import ndimage
import shutil, os
from datetime import timedelta, datetime
import numpy as np
import tempfile
import glob
import re
from copy import deepcopy
from IPython.display import Image, display
import math
%matplotlib inline
import pyproj
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
# Input the range of dates and time wanted for the collection of images
start_year = 2005
start_day = 9
start_month = 12
start_hour = 21
start_minute = 1
start_second = 0

end_year = 2005
end_month = 12
end_day = 9
end_hour = 23
end_minute = 50
end_second = 0

data_path = '/home/rjackson/multidop_grids/'<jupyter_output><empty_output><jupyter_text>This looks for all of the available timeperiods in the data_path directory and pulls outthe file names that match the given time periods above.<jupyter_code># get_radar_times
#     start_year = Start year of animation
#     start_month = Start month of animation
#     start_day = Start day of animation
#     start_hour = Start hour of animation
#     end_year = End year of animation
#     end_month = End month of animation
#     end_day = End day of animation
#     end_minute = End minute of animation
#     minute_interval = Interval in minutes between scans (default is 5)
# This procedure acquires an array of Radar classes between start_time and end_time  
def get_dda_times(start_year, start_month, start_day,
                  start_hour, start_minute, end_year,
                  end_month, end_day, end_hour, 
                  end_minute, minute_interval=5):

    start_time = datetime(start_year,
                      start_month,
                      start_day,
                      start_hour,
                      start_minute,
                      )
    end_time = datetime(end_year,
                      end_month,
                      end_day,
                      end_hour,
                      end_minute,
                      )

    deltatime = end_time - start_time


    if(deltatime.seconds > 0 or deltatime.minute > 0):
        no_days = deltatime.days + 1
    else:
        no_days = deltatime.days
    
    if(start_day != end_day):
        no_days = no_days + 1
        
    days = np.arange(0, no_days, 1)
    print('We are about to load grid files for ' + str(no_days) + ' days')
    

    # Find the list of files for each day
    cur_time = start_time
 
    file_list = []
    time_list = []
    for i in days:
        year_str = "%04d" % cur_time.year
        day_str = "%02d" % cur_time.day
        month_str = "%02d" % cur_time.month
        format_str = (data_path +
                      'cf_compliant_grid' +
                      year_str +
                      month_str +
                      day_str +
                      '*.nc')
    
    
        print('Looking for files with format ' + format_str)
          
        data_list = glob.glob(format_str)
        
        for j in range(0, len(data_list)):
            file_list.append(data_list[j])
        cur_time = cur_time + timedelta(days=1)
    
    
    # Parse all of the dates and time in the interval and add them to the time list
    past_time = []
    for file_name in file_list:
        date_str = file_name[-15:-3]
        year_str = date_str[0:4]
        month_str = date_str[4:6]
        day_str = date_str[6:8]
        hour_str = date_str[8:10]
        minute_str = date_str[10:12]
        second_str = '00'
             
        cur_time = datetime(int(year_str),
                            int(month_str),
                            int(day_str),
                            int(hour_str),
                            int(minute_str),
                            0)
        time_list.append(cur_time)
        
    
    # Sort time list and make sure time are at least xx min apart
    time_list.sort()
    time_list_sorted = deepcopy(time_list)
   
    time_list_final = []
    past_time = []
   
    for times in time_list_sorted:         
        cur_time = times  
        
        if(past_time == []):
            past_time = cur_time
            
        if(cur_time - past_time >= timedelta(minutes=minute_interval)
           and cur_time >= start_time and cur_time <= end_time):
            
            time_list_final.append(cur_time)
            past_time = cur_time
                           
    return time_list_final

# Get a Radar object given a time period in the CPOL dataset
def get_grid_from_dda(time):
    year_str = "%04d" % time.year
    month_str = "%02d" % time.month
    day_str = "%02d" % time.day
    hour_str = "%02d" % time.hour
    minute_str = "%02d" % time.minute
    second_str = "%02d" % time.second
    file_name_str = (data_path +
                    'cf_compliant_grid' +
                     year_str +
                     month_str +
                     day_str +
                     hour_str +
                     minute_str + '.nc')
    
    radar = pyart.io.read_grid(file_name_str)
    return radar
# Plot the radars from given time.

times = get_dda_times(start_year, start_month, start_day,
                      start_hour, start_minute, end_year,
                      end_month, end_day, end_hour, 
                      end_minute, minute_interval=0)
def dms_to_decimal(degrees, minutes, seconds):
    if(degrees > 0):
        return degrees+minutes/60+seconds/3600
    else:
        return degrees-minutes/60-seconds/3600

def get_bca(grid):
    berr_origin = [-12960.1,-23091.1]
    x,y = np.meshgrid(grid.x['data'], grid.y['data'])
    a = np.sqrt(np.multiply(x,x)+np.multiply(y,y))
    b = np.sqrt(pow(x-berr_origin[0],2)+pow(y-berr_origin[1],2))
    c = np.sqrt(berr_origin[0]*berr_origin[0]+berr_origin[1]*berr_origin[1])
    theta_1 = np.arccos(x/a)
    theta_2 = np.arccos((x-berr_origin[1])/b)
    return np.arccos((a*a+b*b-c*c)/(2*a*b))

# Gets beam crossing angle over 2D grid centered over Radar 1. 
# grid_x, grid_y are cartesian coordinates from pyproj.Proj (or basemap)
def get_bca(rad1_lon, rad1_lat, 
            rad2_lon, rad2_lat,
            grid_lon, grid_lat):
    # Beam crossing angle needs cartesian coordinates
    p = ccrs.PlateCarree()
    p = p.as_geocentric()
    rad1 = p.transform_points(ccrs.PlateCarree().as_geodetic(),
                              np.array(rad1_lon), 
                              np.array(rad1_lat))
    rad2 = p.transform_points(ccrs.PlateCarree().as_geodetic(),
                              np.array(rad2_lon), 
                              np.array(rad2_lat))
    grid_lon, grid_lat = np.meshgrid(grid_lon, grid_lat)
    grid = p.transform_points(ccrs.PlateCarree().as_geodetic(),
                              grid_lon,
                              grid_lat,
                              np.zeros(grid_lon.shape))
    # Create grid with Radar 1 in center
    x = grid[:,:,0]-rad1[0,0]
    y = grid[:,:,1]-rad1[0,1]
    rad2 = rad2 - rad1
    a = np.sqrt(np.multiply(x,x)+np.multiply(y,y))
    b = np.sqrt(pow(x-rad2[0,0],2)+pow(y-rad2[0,1],2))
    c = np.sqrt(rad2[0,0]*rad2[0,0]+rad2[0,1]*rad2[0,1])
    theta_1 = np.arccos(x/a)
    theta_2 = np.arccos((x-rad2[0,1])/b)
    return np.arccos((a*a+b*b-c*c)/(2*a*b))<jupyter_output><empty_output><jupyter_text>Grid plotting codeThis code creates plots from all of the Grids developed by multidop This loads the Grid files and creates the animation<jupyter_code># Reflectivity plot callback functions
def plot_dd_lobe(frame_number):

    plt.clf()   
    pyart_grid = get_grid_from_dda(times[frame_number])    
    level = 6    
    CPOLGridDisplay = pyart.graph.GridMapDisplay(pyart_grid)
    CPOLGridDisplay.plot_basemap(min_lat=-12.8, max_lat=-12.0,
                                 min_lon=130.4, max_lon=131.6,
                                 auto_range=False,
                                 lat_lines=[-12.8, -12.6, -12.4, -12.2, -12.0],
                                 lon_lines=[130.4, 130.6, 130.8, 131.0, 131.2, 131.4, 131.6],
                                 resolution='h')
    CPOLGridDisplay.basemap.drawmapscale(131.2, -12.1, 130.6, -12.6, 20)
    CPOLGridDisplay.basemap.fillcontinents()
    lons, lats = pyart_grid.get_point_longitude_latitude(level=level)
     
    # Plot every second arrow
    lons = lons[0::2,0::2]
    lats = lats[0::2,0::2]
    bca = get_bca(pyart_grid)
    bca = bca[0::2, 0::2]
    CSa = CPOLGridDisplay.basemap.contour(lons, lats,
                                          bca, latlon='True',
                                          levels=[math.pi/6, 5*math.pi/6],
                                          linewidths=2)
    
    cpol_latitude = -12.249166
    cpol_longitude = 131.04445
    berr_latitude = -12.456944
    berr_longitude = 130.925
    berr_x, berr_y = CPOLGridDisplay.basemap(berr_longitude, berr_latitude)
    plt.annotate('Berrima', xy=(berr_x, berr_y), 
                 xycoords='data', fontweight='bold',
                 fontsize=20)
    cpol_x, cpol_y = CPOLGridDisplay.basemap(cpol_longitude, cpol_latitude)
    plt.annotate('CPOL', xy=(cpol_x, cpol_y), 
                 xycoords='data', fontweight='bold',
                 fontsize=20)

def scale_bar(ax, length, location=(0.5, 0.05), linewidth=3):
    """
    ax is the axes to draw the scalebar on.
    location is center of the scalebar in axis coordinates ie. 0.5 is the middle of the plot
    length is the length of the scalebar in km.
    linewidth is the thickness of the scalebar.
    """
    #Projection in metres, need to change this to suit your own figure
    utm = ccrs.UTM(18)
    #Get the extent of the plotted area in coordinates in metres
    x0, x1, y0, y1 = ax.get_extent(utm)
    #Turn the specified scalebar location into coordinates in metres
    sbcx, sbcy = x0 + (x1 - x0) * location[0], y0 + (y1 - y0) * location[1]
    #Generate the x coordinate for the ends of the scalebar
    bar_xs = [sbcx - length * 500, sbcx + length * 500]
    #Plot the scalebar
    ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k', linewidth=linewidth)
    #Plot the scalebar label
    ax.text(sbcx, sbcy, str(length) + ' km', transform=utm,
            horizontalalignment='center', verticalalignment='bottom')

def plot_dd_lobes(radar1_loc, radar2_loc, radar1_name, radar2_name):
    ax = plt.axes(projection=ccrs.PlateCarree())
    # Amf locations
    grid_lon = np.arange(radar1_loc[0]-1.5, radar1_loc[0]+1.5, 0.01)
    grid_lat = np.arange(radar1_loc[1]-1.5, radar1_loc[1]+1.5, 0.01)

    i5 = [dms_to_decimal(-97, 35, 37.68), dms_to_decimal(36, 29, 29.4)]
    i4 = [dms_to_decimal(-91, 21, 49.32), dms_to_decimal(36, 34, 44.4)]
    bca = get_bca(radar1_loc[0], radar1_loc[1],
                  radar2_loc[0], radar2_loc[1],
                  grid_lon, grid_lat)
    lon_gridded, lat_gridded = np.meshgrid(grid_lon, grid_lat)
    # Create a feature for States/Admin 1 regions at 1:50m from Natural Earth
    states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')

    SOURCE = 'Natural Earth'
    LICENSE = 'public domain'

    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(states_provinces, edgecolor='gray')
    
    ax.set_xticks(grid_lon[::int(len(grid_lon)/5)], crs=ccrs.PlateCarree())
    ax.set_yticks(grid_lat[::int(len(grid_lon)/5)], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    plt.contour(lon_gridded, lat_gridded, 
                bca, 
                levels=[math.pi/6, 5*math.pi/6],
                linewidths=2,
                transform=ccrs.PlateCarree())
    plt.annotate('ENX', 
                 xy=(ENX[0]+0.02, ENX[1]+0.01),
                 fontweight='bold', 
                 fontsize=8, 
                 transform=ccrs.PlateCarree())
    plt.annotate('KIND', 
                 xy=(KIND[0]+0.02, KIND[1]+0.01),
                 fontweight='bold', 
                 fontsize=8, 
                 transform=ccrs.PlateCarree())
    plt.annotate('Albany', 
                 xy=[Albany_airport[0]+0.02, Albany_airport[1]+0.01],
                 fontweight='bold',
                 fontsize=8, 
                 transform=ccrs.PlateCarree())
    plt.annotate('Delmar', 
                 xy=[Delmar[0]+0.02, Delmar[1]+0.01], 
                 fontweight='bold',
                 fontsize=8, 
                 transform=ccrs.PlateCarree())        
    plt.annotate('VOOR', 
                 xy=[VOOR[0]+0.02, VOOR[1]+0.01], 
                 fontweight='bold',
                 fontsize=8, 
                 transform=ccrs.PlateCarree())        
    plt.annotate('MEDU', 
                 xy=[MEDU[0]+0.02, MEDU[1]+0.01], 
                 fontweight='bold',
                 fontsize=8, 
                 transform=ccrs.PlateCarree())
    plt.plot(ENX[0], ENX[1], marker='d', 
             linewidth=1, color='k')
    plt.plot(KIND[0], KIND[1], marker='d', 
             linewidth=1, color='k')
    plt.plot(Albany_airport[0], Albany_airport[1], 
             linewidth=1, marker='d', color='k')
    plt.plot(Delmar[0], Delmar[1], marker='d', 
             linewidth=1, color='k')
    plt.plot(VOOR[0], VOOR[1], marker='d', 
             linewidth=1, color='k')
    plt.plot(MEDU[0], MEDU[1], marker='d', 
             linewidth=1, color='k')
    scale_bar(ax, 20, location=(0.1, 0.9),)
    ax.coastlines(resolution='10m')
    ax.stock_img()
    plt.xlim((grid_lon[0], grid_lon[-1]))
    plt.ylim((grid_lat[0], grid_lat[-1]))<jupyter_output><empty_output><jupyter_text>Plot DD lobes for XSAPR i4 and i5<jupyter_code>ENX = [dms_to_decimal(-74.0, 3.0, 50.0), dms_to_decimal(42.0, 35.0, 11.0)]
KIND = [dms_to_decimal(-73.0, 42.0, 56.0), dms_to_decimal(42.0, 24.0, 36.0)]

# Amf locations
i5 = [dms_to_decimal(-97, 35, 37.68), dms_to_decimal(36, 29, 29.4)]
i4 = [dms_to_decimal(-91, 21, 49.32), dms_to_decimal(36, 34, 44.4)]

plt.figure(figsize=(8,10))
plot_dd_lobes(i4, i5, 'i4', 'i5')
plt.title('KIND')
radar_file = '/home/rjackson/xsapr/raw_files/sgpxsaprI5.00.20170801.001102.mnt.XSW170731231005.RAWG2TS.maint'
radar = pyart.io.read(radar_file)
print(radar.fields.keys())
max_lat=37.0
min_lat=36.0
max_lon=-97.0
min_lon=-98.3
sweep=3
i5 = [dms_to_decimal(-97, 35, 37.68), dms_to_decimal(36, 29, 29.4)]
i4 = [dms_to_decimal(-97, 21, 49.32), dms_to_decimal(36, 34, 44.4)]
grid_lat = np.arange(min_lat, max_lat, 0.01)
grid_lon = np.arange(min_lon, max_lon, 0.01)
                                 
bca = get_bca(i4[0], i4[1], i5[0], i5[1], grid_lon, grid_lat)
grid_lon, grid_lat = np.meshgrid(grid_lon, grid_lat)
display = pyart.graph.RadarMapDisplayCartopy(radar)
fig = plt.figure(figsize=[10, 8])
display.plot_ppi_map('reflectivity', sweep=sweep, vmin=-8, vmax=64,
                         min_lon=min_lon, max_lon=max_lon, min_lat=min_lat,
                         max_lat=max_lat, resolution='50m',
                         cmap=pyart.graph.cm.NWSRef,
                         lat_lines=np.arange(min_lat, max_lat+0.2, 0.2),
                         lon_lines=np.arange(min_lon, max_lon+0.2, 0.2),
                         projection=ccrs.PlateCarree())
plt.contour(grid_lon, grid_lat, bca, levels=[np.pi/6, 5*np.pi/6], linewidths=2, colors='k')<jupyter_output>dict_keys(['spectrum_width', 'velocity', 'specific_differential_phase', 'normalized_coherent_power', 'differential_phase', 'reflectivity', 'total_power', 'cross_correlation_ratio', 'differential_reflectivity'])
StarcoderdataJupyterStructuredCleanDedup