<jupyter_start><jupyter_text>Look at WG met file (subset of WG L2 data file)  Working toward estimation of bulk fluxesWill need:  - lat/lon (1 Hz and 20 Hz)- SMP21/SGR4 data- wind (u,v,w)- all WXT variables- SST, SSS, sensor depth ("UCTD")1 Hz variables I need:- time_1Hz- latitude_1Hz- longitude_1Hz- UCTD_sea_water_temperature- UCTD_salinity- UCTD_density- UCTD_depth- WXT_rainfall_rate- WXT_atmospheric_pressure- WXT_air_temperature- WXT_relative_humidity- WXT_wind_speed- WXT_wind_direction- SMP21_shortwave_flux- SGR4_longwave_flux- _also need some estimate of surface current_20 Hz variables I need:- time_20Hz- latitude_20Hz- longitude_20Hz- wind_speed- wind_direction- wind_verticalNot sure if I need:  - time_15min- wave_frequency<jupyter_code>import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import matplotlib
import datetime as dt
################
# This allows us to import Tom_tools
import sys
sys.path.append('../../Tom_tools/') # you may need to adjust this path
# sys.path.append('../SWOT_IW_SSH/jtf/Tom_tools/') # you may need to adjust this path
import Tom_tools_v1 as tt
################
# %matplotlib inline
%matplotlib qt 
savefig = True # set to true to save plots as file
plt.rcParams['figure.figsize'] = (5,4)
plt.rcParams['figure.dpi'] = 150
plt.rcParams['savefig.dpi'] = 400
plt.close('all')
plt.rcParams['axes.xmargin'] = 0


__figdir__ = '../plots/WG_timeseries/'
savefig_args = {'bbox_inches':'tight', 'pad_inches':0.2}
plotfiletype='png'
path = '/mnt/e/SMODE_data/pilot/WG/L2/'
WG = 'WHOI43'
file = 'SMODE_PFC_Wavegliders_'+WG+'.nc'
# This needs to be set before any dates are encoded as datetime64
# Because the time increment is so small, we get 'out of range' when starting from the standard datetime64 epoch of 1970
# matplotlib.dates.set_epoch('2000-01-01T00:00:00') 
%%time
ds = xr.open_dataset(path+file, engine = 'netcdf4', decode_times = True) #decode_times = False, 
ds
# Raw met plot from WG:
fig, axs = plt.subplots(5, 1, sharex=True)
fig.autofmt_xdate()
plt.subplot(5,1,1)
h1, = plt.plot(ds.time_1Hz, ds.WXT_air_temperature)
h2, = plt.plot(ds.time_1Hz, ds.UCTD_sea_water_temperature)
plt.legend([h1, h2],['Air temp.','SST'])
plt.ylabel('T [$^\circ$C]')
plt.title('Raw 1 Hz WXT measurements')

plt.subplot(5,1,2)
plt.plot(ds.time_1Hz, ds.WXT_relative_humidity)
plt.ylabel('[%]')
plt.legend(['Rel. Humidity'])

plt.subplot(5,1,3)
plt.plot(ds.time_15min, ds.wave_significant_height)
plt.ylabel('[m]')
plt.legend(['Sig. wave height'],loc='upper right')

plt.subplot(5,1,4)
plt.plot(ds.time_20Hz, ds.wind_speed)
plt.plot(ds.time_1Hz, ds.WXT_wind_speed)
plt.ylabel('[m/s]')
plt.legend(['Gill Wind speed','WXT wind speed'],loc='upper right')

plt.subplot(5,1,5)
plt.plot(ds.Workhorse_time, ds.Workhorse_altitude)
plt.ylabel('[m]')
plt.legend(['Workhorse altitude (from IMU/GPS)'],loc='upper right')

if savefig:
    plt.savefig(__figdir__+WG+'_raw_met' + '.' +plotfiletype,**savefig_args)

# Raw met plot from WG:
fig, axs = plt.subplots(4, 1, sharex=True)
fig.autofmt_xdate()
plt.subplot(4,1,1)
plt.plot(ds.time_1Hz, ds.WXT_atmospheric_pressure)
plt.ylabel('[mbar]')
plt.legend(['Atm. pressure'])

plt.subplot(4,1,2)
plt.plot(ds.time_1Hz, ds.SGR4_longwave_flux)
plt.ylabel('W/m^2')
plt.legend(['Longwave radiation'])

plt.subplot(4,1,3)
plt.plot(ds.time_1Hz, ds.SMP21_shortwave_flux)
plt.ylabel('W/m^2')
plt.legend(['Shortwave radiation'])

plt.subplot(4,1,4)
plt.plot(ds.time_1Hz, ds.WXT_rainfall_rate)
plt.ylabel('[mm/hr]')
plt.legend(['Precip. rate'])

if savefig:
    plt.savefig(__figdir__+WG+'_raw_met2' + '.' +plotfiletype,**savefig_args)

# Raw met plot from WG:
fig, axs = plt.subplots(5, 1, sharex=True,figsize=(5,7))
fig.autofmt_xdate()
plt.subplot(5,1,1)
h1, = plt.plot(ds.time_1Hz, ds.WXT_air_temperature)
h2, = plt.plot(ds.time_1Hz, ds.UCTD_sea_water_temperature)
plt.legend([h1, h2],['Air temp.','SST'])
plt.ylabel('T [$^\circ$C]')
plt.title('Wave Glider surface measurements')

plt.subplot(5,1,2)
plt.plot(ds.time_1Hz, tt.run_avg1d(ds.SGR4_longwave_flux,15*60))
plt.ylabel('W/m^2')
plt.legend(['Longwave radiation'])
plt.setp(plt.gca().get_xticklabels(), visible=False)

plt.subplot(5,1,3)
plt.plot(ds.time_1Hz, tt.run_avg1d(ds.SMP21_shortwave_flux,15*60))
plt.ylabel('W/m^2')
plt.legend(['Shortwave radiation'])

plt.subplot(5,1,4)
plt.plot(ds.time_1Hz, tt.run_avg1d(ds.WXT_wind_speed,15*60))
plt.ylabel('[m/s]')
plt.legend(['wind speed'],loc='upper right')

plt.subplot(5,1,5)
plt.plot(ds.time_15min, ds.wave_significant_height)
plt.ylabel('[m]')
plt.legend(['Sig. wave height'],loc='upper right')

if savefig:
    plt.savefig(__figdir__+WG+'_raw_met3' + '.' +plotfiletype,**savefig_args)

# There is a limited number of NaNs in WXT wind speed
fig, axs = plt.subplots(1, 1)
fig.autofmt_xdate()

plt.plot(ds.time_1Hz, np.isnan(ds.WXT_wind_speed.values))
plt.title('NaNs in WXT wind speed')
# Let's interpolate them
ff = np.where(np.isnan(ds.WXT_wind_speed.values)==0)
t = ds.time_1Hz[ff]
wspd = ds.WXT_wind_speed[ff]
WXT_wspd_interp = np.interp(ds.time_1Hz, t, wspd)
'''
fig, axs = plt.subplots(1, 1)
fig.autofmt_xdate()

# plt.plot(t, np.isnan(wspd))
plt.plot(np.isnan(WXT_wspd_interp))
plt.title('NaNs in WXT wind speed after trying to remove nans')
''' # this checks out
# There is also a limited number of NaNs in Gill wind speed
fig, axs = plt.subplots(1, 1)
fig.autofmt_xdate()

plt.plot(ds.time_20Hz, np.isnan(ds.wind_speed.values))
plt.title('NaNs in Gill wind speed')
# Let's interpolate them
ff = np.where(np.isnan(ds.wind_speed.values)==0)
t = ds.time_20Hz[ff]
wspd = ds.wind_speed[ff]
wspd_interp = np.interp(ds.time_20Hz, t, wspd)
fig, axs = plt.subplots(1, 1)
M=151
tt.spectrum_band_avg(wspd_interp,1/20,M,winstr=None,plotflag=True,ebarflag=None)
tt.spectrum_band_avg(WXT_wspd_interp,1,M,winstr=None,plotflag=True,ebarflag=False)
plt.title('Gill, WXT wind speed')
plt.xlabel('Hz')
plt.ylabel('Spectral density [m$^2$/s$^2$/Hz]')
plt.legend(['Gill Wind speed','WXT wind speed'],loc='upper right')
if savefig:
    plt.savefig(__figdir__+WG+'_wsp_spectra' + '.' +plotfiletype,**savefig_args)<jupyter_output><empty_output><jupyter_text>There appears to be some residual wave contamination in the motion corected velocities, but it shouldn't be a big issue when we average to 1 min or longer<jupyter_code>fig, axs = plt.subplots(1, 1, sharex=True)
fig.autofmt_xdate()
plt.plot(ds.time_20Hz, tt.run_avg1d(ds.wind_speed,15*60*20))
plt.plot(ds.time_1Hz, tt.run_avg1d(ds.WXT_wind_speed,15*60))
plt.ylabel('[m/s]')
plt.legend(['Gill Wind speed ('+str(ds.INST_Gill.height_m)+'m)','WXT wind speed ('+str(ds.INST_WXT.height_above_waterline_m)+'m)'],loc='upper right')
plt.title('15-min average Gill R-3 and WXT wind speed')
plt.axis('tight')
if savefig:
    plt.savefig(__figdir__+WG+'_wspd_smoothed' + '.' +plotfiletype,**savefig_args)

str(ds.INST_Gill.height_m)
ds.INST_WXT.height_above_waterline_m
fig, axs = plt.subplots(1, 1)
fig.autofmt_xdate()

plt.plot(ds.Workhorse_time, ds.Workhorse_altitude)
plt.title('Workhorse altitude\n Something bad happened with the VectorNav?')
plt.xlim(18935.65223174665, 18936.366246499798)

if savefig:
    plt.savefig(__figdir__+WG+'_Workhorse_altitude' + '.' +plotfiletype,**savefig_args)

xlim = plt.xlim()
xlim<jupyter_output><empty_output><jupyter_text>Make an averaged data set (say, 1 min)It would probably be smarter to automate this to loojp through all ds.time_1Hz variables (for example).  That might look something like this:  ```var_list = []  for var in ds.data_vars.variables:      var_list.append(var)      print(var)  ```<jupyter_code># compare different time vectors to decide how to line up variables from different time bases
# ds.time_20Hz
# ds.time_1Hz
# ds.Workhorse_time

# Try the simple thing
# xy, x_ind, y_ind = np.intersect1d(ds.time_1Hz, ds.time_20Hz, assume_unique=True, return_indices=True)
# That doesn't work (no surprise) because the times don't line up exactly (possibly numerical precision issue)

# I think they are already nearly lined up
ds.time_1Hz[0:10:1]-ds.time_20Hz[0:200:20]
# so, try 0,20,40... in 20 Hz array
tfoo=[]
for n in range(len(ds.time_1Hz)):
    tfoo.append(ds.time_20Hz[0:200:20]
len(ds.time_1Hz)
nsec=60*15 # 15 minute average
# ds.time_20Hz variables
# wdir is positive clockwise from North
# need to make U, V
Gill_U=ds.wind_speed*np.cos(ds.wind_direction/np.pi/2)
Gill_V=ds.wind_speed*np.sin(ds.wind_direction/np.pi/2)
Gill_wspd_low = tt.run_avg1d(ds.wind_speed, nsec*20)
Gill_U_low = tt.run_avg1d(Gill_U,nsec*20)
Gill_V_low = tt.run_avg1d(Gill_V,nsec*20)
pitch_low = tt.run_avg1d(ds.pitch,nsec*20)
roll_low = tt.run_avg1d(ds.roll,nsec*20)

# ds.time_1Hz variables
WXT_U=ds.WXT_wind_speed*np.cos(ds.WXT_wind_direction/np.pi/2)
WXT_V=ds.WXT_wind_speed*np.sin(ds.WXT_wind_direction/np.pi/2)
WXT_wspd_low = tt.run_avg1d(ds.WXT_wind_speed,nsec)
WXT_atmp_low = tt.run_avg1d(ds.WXT_air_temperature,nsec)
WXT_rh_low = tt.run_avg1d(ds.WXT_relative_humidity,nsec)
swr_low = tt.run_avg1d(ds.SMP21_shortwave_flux,nsec)
lwr_low = tt.run_avg1d(ds.SGR4_longwave_flux,nsec)
lat_low = tt.run_avg1d(ds.latitude_1Hz,nsec)
lon_low = tt.run_avg1d(ds.longitude_1Hz,nsec)

# Workhorse_time variables
U_low=tt.run_avg2d(ds.Workhorse_vel_east,nsec,1)
V_low=tt.run_avg2d(ds.Workhorse_vel_north,nsec,1)
W_low=tt.run_avg2d(ds.Workhorse_vel_up,nsec,1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15475297395777515, 'fraction_numerical': 0.03600900225056264, 'mean_word_length': 6.290625, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '586857', 'license': 'MIT', 'n_tokens_mistral': 4238, 'n_tokens_neox': 3934, 'n_words': 728, 'path': 'code/WG_L2b_look.ipynb', 'repo_name': 'jtomfarrar/S-MODE_analysis', 'size': '9331'}