# -*- coding: utf-8 -*- """ Created on Mon May 10 09:51:55 2021 @author: ppantina """ # -*- coding: utf-8 -*- """ Created on Mon Dec 21 11:35:09 2020 @author: ppantina """ # -*- coding: utf-8 -*- """ Created on Mon Nov 9 10:38:04 2020 @author: ppantina """ import glob from netCDF4 import MFDataset import h5py import numpy as np from scipy import interpolate def L1B_sub_add_horiz_wind (radName, flightLines, L1B_path, h_filepath): ##Isolate data files L1B_files = np.sort(glob.glob(L1B_path )) h_files = np.sort(glob.glob(h_filepath)) ##Bug to get MFDataset to read the file array h_files2 = [] for i in range(len(h_files)): h_files2.append(h_files[i]) ##Read in datasets hrrr = MFDataset(h_files2) L1B = h5py.File(L1B_files[0], 'r') ##Save some vars for later L1B_time = L1B['Time'] ['Data'] ['TimeUTC'] [:] L1B_latitude = L1B['Navigation']['Data'] ['Latitude'] [:] L1B_longitude= L1B['Navigation']['Data'] ['Longitude'][:] L1B_height = L1B['Navigation']['Data'] ['Height'][:] L1B_range = L1B['Products'] ['Information']['Range'] [:] ##Create storage array for ultimate u/v wind. L1B_uwind = np.zeros((len(L1B_time), len(L1B_range)), dtype = np.float32) * np.nan L1B_vwind = np.zeros((len(L1B_time), len(L1B_range)), dtype = np.float32) * np.nan ##Convert hrrr ISO times to Linux times hrrr_time = hrrr['time'][:] ##Begin calculation for ff in range(len(flightLines)): ##for all flightlines: print ('...Flight Leg', ff, 'of', len(flightLines)-1) ##Find start/stop index of the flightline in the larger time array. start = np.argmin(np.abs(L1B_time - flightLines[ff][0] )) stop = np.argmin(np.abs(L1B_time - flightLines[ff][-1])) ##Select the times/lat/lon at which to calculate u/v winds. ##We will select a few per flight leg, then interpolate between these. ##Dont expect the windfield to vary much here. i_time = np.hstack((L1B_time [start:stop:60], L1B_time [stop])) i_lat = np.hstack((L1B_latitude [start:stop:60], L1B_latitude [stop])) i_lon = np.hstack((L1B_longitude[start:stop:60], L1B_longitude[stop])) ##Interpolate WIND and HEIGHT to specific points in time/lat/lon ##using the interpn function pres_levs = hrrr['hgt'].shape[1] ##number of hrrr levels num_times = round(len(i_time)) ##number of timesteps to interpolate interp_coarse_uwind = np.zeros((num_times,pres_levs)) ##storage array interp_coarse_vwind = np.zeros((num_times,pres_levs)) ##storage array interp_coarse_height = np.zeros((num_times,pres_levs)) ##storage array ##For all times/pressure levels, interpolate WIND and HEIGHT and remove LAT/LON dimensions ##The resulting arrays will be coarse functions of time and height. for tt in range(num_times): #for all times print ('......Interpolating segment ', tt, 'of', num_times-1) for pp in range(pres_levs): #for all levs interp_coarse_uwind [tt, pp] = interpolate.interpn((hrrr_time, hrrr['lon'][:,0], hrrr['lat'][0,:]), hrrr['uwnd'][:,pp,:,:], [i_time[tt], i_lon[tt], i_lat[tt]]) interp_coarse_vwind [tt, pp] = interpolate.interpn((hrrr_time, hrrr['lon'][:,0], hrrr['lat'][0,:]), hrrr['vwnd'][:,pp,:,:], [i_time[tt], i_lon[tt], i_lat[tt]]) interp_coarse_height[tt, pp] = interpolate.interpn((hrrr_time, hrrr['lon'][:,0], hrrr['lat'][0,:]), hrrr['hgt' ][:,pp,:,:], [i_time[tt], i_lon[tt], i_lat[tt]]) ##Declare higher resolution time/height arrays. if 'HIWRAP' in radName: new_time = np.arange(i_time[0], i_time[-1], 0.5) else: new_time = np.arange(i_time[0], i_time[-1], 0.25) new_height = np.arange(0, 20000, 10) ##For each coarse time, interpolate vertically to the higher-res height array ##The resulting arrays will be coarse functions of time and fine functions of height. interp_hires_heights_uwind = np.zeros((num_times,len(new_height))) interp_hires_heights_vwind = np.zeros((num_times,len(new_height))) for i, j in enumerate(i_time): interp_hires_heights_uwind[i,:] = np.interp(new_height, interp_coarse_height[i,:], interp_coarse_uwind[i,:]) interp_hires_heights_vwind[i,:] = np.interp(new_height, interp_coarse_height[i,:], interp_coarse_vwind[i,:]) ##For each (new) height, interpolate in time to the higher-res time array. ##The resulting arrays will be fine functions of time and height. interp_hires_uwind = np.zeros((len(new_time),len(new_height))) interp_hires_vwind = np.zeros((len(new_time),len(new_height))) for i, j in enumerate(new_height): interp_hires_uwind[:,i] = np.interp(new_time, i_time, interp_hires_heights_uwind[:,i]) interp_hires_vwind[:,i] = np.interp(new_time, i_time, interp_hires_heights_vwind[:,i]) ##Calculate the height of each range gate, for all times, since we are only provided with "range" ##This is on the same grid as the original L1B data. L1B_altitude = np.zeros((len(L1B_height), len(L1B_range))) for i, j in enumerate(L1B_altitude): #for all times. ##The height is just aircraft height - range L1B_altitude[i]= L1B_height[i] - L1B_range L1B_altitude[i]= L1B_height[i] - L1B_range ##Get rid of negative altitudes. L1B_altitude[L1B_altitude<0 ] = np.nan ##Now, create a U/V wind array with same dimensions as the L1B data. for i, j in enumerate(new_time): #for all times ##Find the index where our flight line time lies within the larger L1B time. if i%1000 == 0: print ('.........Creating wind array', i, 'of', len(new_time)) idx = np.where(L1B_time == new_time[i])[0][0] ##For each altitude at time [i], find the index at which ##it is closest to a height in the hires height array for k, l in enumerate(L1B_altitude[idx]): if np.isfinite(l): idx2 = np.argmin(np.abs(l - new_height)) ##Select the winds at this point. L1B_uwind[idx,k] = interp_hires_uwind[i,idx2] L1B_vwind[idx,k] = interp_hires_vwind[i,idx2] return(L1B_uwind, L1B_vwind)