def skyler_processing(filename): #import the needed packages import numpy as np import xarray as xr #to open the desired netCDF and convert it into a dataframe import pandas as pd #to handle the opened dataframe import matplotlib.pyplot as plt #to plot/represent the data visually import numpy as np #for any necessary calculations import pdb #python debugger yyyymmdd = filename[51:59] #UTC time, EDT = UTC - 4 hours hhmmss = filename[60:66] fff = filename[67:70] #fraction of second ele = filename.partition(".")[0][-3:-1] #opens netCDF, ds contains information on the attributes of the netCDF ds = xr.open_dataset(filename) # pdb.set_trace() #now step through next lines, one-by-one by entering "n" at command prompt #save the filename as a string to attach to the plot title and saved image name #for each file, the naming convention follows: LPR2_SN14_20230623_210409_571-0-37-180.netcdf sname = 'skyler-scan_'+str(yyyymmdd)+'_'+str(hhmmss)+'.'+str(fff)+'_'+str(ele)+'.png' tname = 'SKYLER-HU '+str(yyyymmdd)+' '+str(hhmmss)+'.'+str(fff)+' UTC'+' Elevation '+str(ele)+'\u00b0' #converts xarray dataset into pandas dataframe srad = ds.to_dataframe() #Grab the positional data as well as the three categories of interest #For now, we have not used the Elevation Data - later it will be necessary to find the altitude of each measurement #Note that this dataframe has two indices: Gate(distance from the radar) and Radial(looking direction) data = srad[['Azimuth','Elevation','Reflectivity','Velocity','SignalToNoiseRatio']].reset_index() #turn the gates indices from the data into a physical distance #from the header, the gate width is 23.98 m, and the first gate is located 911.4 m in front of the radar data['Distance (m)'] = [911.4 + (23.98*i) for i in data['Gate']] #This finds data points with SNR < 1 and data points with Reflectivity < 0 to replace with nan - this may change as we understand more about the noise floor of Skyler #I pulled this out of the dataframe any work with arrays instead - prevents errors from being thrown sort_data = np.array([data['Reflectivity'],data['Velocity'],data["SignalToNoiseRatio"]]) low_snr = np.where(sort_data[2] < 1) low_refl = np.where(sort_data[0] < 0) #turn the found locations into nan for Reflectivity sort_data[0][low_snr] = np.nan*len(sort_data[0][low_snr]) sort_data[0][low_refl] = np.nan*len(sort_data[0][low_refl]) #turn all locations with nan Reflectivity into nan Velocity no_refl = np.where(np.isnan(sort_data[0])) sort_data[1][no_refl] = np.nan*len(sort_data[1][no_refl]) #drop the old Reflectivity and Velocity data, and insert the new arrays drefl = data.drop(columns = ['Reflectivity','Velocity']) drefl['Reflectivity'] = sort_data[0] drefl['Velocity'] = sort_data[1] snr_1 = drefl.reset_index(drop = True) #add in x y and z coordinates #These are currently in meters ## snr_1['dx'] = [(snr_1['Distance (m)'][i])*np.sin(np.pi*(snr_1['Azimuth'][i])/180)*np.cos(np.pi*(int(filename[-10:-8]))/180) for i in range(len(snr_1['Azimuth']))] ## snr_1['dy'] = [(snr_1['Distance (m)'][i])*np.cos(np.pi*(snr_1['Azimuth'][i])/180)*np.cos(np.pi*(int(filename[-10:-8]))/180) for i in range(len(snr_1['Azimuth']))] #height of the radar is 60.76 m #note the use of int(filename[-10:-8]) - the elevation value! ## snr_1['dz'] = [60.76+snr_1['Distance (m)'][i]*np.sin(np.pi*(int(filename[-10:-8]))/180) for i in range(len(snr_1['Azimuth']))] #add in longitude and latitude coordinates ## snr_1['Latitude'] = [(i/110850) + 37.025833 for i in snr_1['dy']] ## snr_1['Longitude'] = [(snr_1['dx'].loc[j]/(110850*(np.pi*snr_1['Latitude'].loc[j])/180) - 76.34103) for j in range(len(snr_1['dx']))] #Below is only plotting! If we want to interact with data, then we can stop here ##### #Set up Reflectivity data for a contour plot ##### #Grab unique values in both Azimuth and Distance (m) columns u_az = np.unique(snr_1['Azimuth']) u_dist = np.unique(snr_1['Distance (m)']) #Create arrays to work with the velocity and reflectivity data refl_array = [] vel_array = [] #To do a contour plot, I need the data to be in 2D, w/ dimensions of len(Azimuth) and len(Distance (m)) #Right now the data is in 1D w/ dimension len(Azimuth)*len(Distance (m)) for i in u_dist: #using each unique Distance (m) value #Find the Reflectivity and Velocity data at each unique Distance (m) value, there should be no #repeating Azimuth values. x = snr_1['Reflectivity'].loc[np.where(snr_1['Distance (m)'] == i)] y = snr_1['Velocity'].loc[np.where(snr_1['Distance (m)'] == i)] #add the data to the array as a new list refl_array.append(list(x)) vel_array.append(list(y)) #The final array should have dimensions len(Distance (m)) by len(Azimuth) (row by column) #create a dataframe for the reflectivity and velocity that is filled w/ data for each rf2d = pd.DataFrame(refl_array) vel2d = pd.DataFrame(vel_array) #Format the necessary plots, so that the plot represents Earth directional degrees #The full radial plot, where 90 is due East, 0 is due North, and 180 is due South fig,ax = plt.subplots(2,1,figsize = (6,10),subplot_kw = dict(projection = 'polar')) fig.suptitle(tname) #From the filename input to the function ax[0].set_theta_zero_location("N") #Set 0 degrees to be in the North direction ax[0].set_theta_direction(-1) #Set the direction of theta increase to match Earth Directional Degrees ax[0].set_thetamax(np.max(snr_1['Azimuth'])+5) #The choice of 5 here is arbitrary - just to show data edge ax[0].set_thetamin(np.min(snr_1['Azimuth'])-5) ax[1].set_theta_zero_location("N") #Set 0 degrees to be in the North directionz ax[1].set_theta_direction(-1) #Set the direction of theta increase to match Earth Directional Degrees ax[1].set_thetamax(np.max(snr_1['Azimuth'])+5) #The choice of 5 here is arbitrary - just to show data edge ax[1].set_thetamin(np.min(snr_1['Azimuth'])-5) #Printing the label on the radial axis above the mean value of that axis on each plot #Rotate the label to match the angle of the upper image bound label_loc = np.mean(snr_1['Distance (m)']) ax[0].text(0,label_loc,'Meters',rotation = np.min(snr_1['Azimuth'])-5) ax[1].text(0,label_loc,'Meters',rotation = np.min(snr_1['Azimuth'])-5) #the polar projection expects values in radians, so need to convert Azimuth values (given in degrees) rad = np.pi*np.array(u_az)/180 #plot the Reflectivity, Velocity plt.set_cmap('gist_ncar') c_refl = ax[0].contourf(rad,np.array(u_dist),rf2d,levels = np.linspace(0,50,101)) fig.colorbar(c_refl,label = "Reflectivity dBZ",ax = ax[0],ticks = np.linspace(0,50,11)) ax[0].set_title('Reflectivity') c_vel = ax[1].contourf(rad,np.array(u_dist),vel2d,levels = np.linspace(-20,20,41)) fig.colorbar(c_vel,label = "Velocity (m/s)",ax = ax[1],ticks = np.linspace(-20,20,9)) ax[1].set_title('Velocity') plt.savefig(sname) # plt.show()