import matplotlib.pyplot as plt
import xarray as xr
from tqdm import tqdm_notebook
import cmocean.cm as cmo
import pandas as pd
import numpy as np
import glidertools as gt
import my_functions as my
/Users/marcel/opt/anaconda3/envs/duplessis2021_JGR/lib/python3.9/site-packages/xarray/backends/cfgrib_.py:27: UserWarning: Failed to load cfgrib - most likely there is a problem accessing the ecCodes library. Try `import cfgrib` to get the full error message warnings.warn(
dat_saz = xr.open_dataset('../data/sg542_saz_L2.nc')
dat_pfz = xr.open_dataset('../data/slocum_pfz_L2.nc')
dat_miz = xr.open_dataset('../data/sg643_miz_L2.nc')
interp_res = '6H'
offset = '3H'
dt=21600
dat_saz_interp = dat_saz.resample(time=interp_res, loffset=offset).mean()
dat_pfz_interp = dat_pfz.resample(time=interp_res, loffset=offset).mean()
dat_miz_interp = dat_miz.resample(time=interp_res, loffset=offset).mean()
dat_saz_interp=dat_saz_interp.interpolate_na(dim='time')
dat_pfz_interp=dat_pfz_interp.interpolate_na(dim='time')
dat_miz_interp=dat_miz_interp.interpolate_na(dim='time')
for i, dat in enumerate([dat_saz_interp, dat_pfz_interp, dat_miz_interp]):
dat['mld_03'] = (('time'), my.calc_mld(dat.density, dat.depth, den_lim=0.03))
for dat in [dat_saz_interp, dat_pfz_interp, dat_miz_interp]:
xt=[]
xs=[]
for i, m in enumerate(np.round(dat['mld_03'].astype(int))):
xt += dat.temp.sel(depth=slice(5, m-5)).mean(dim='depth').values[i],
xs += dat.salt.sel(depth=slice(5, m-5)).mean(dim='depth').values[i],
dat['ml_s'] = (('time'), xs)
dat['ml_t'] = (('time'), xt)
dat['ml_s_smooth'] = (('time'), gt.cleaning.rolling_window(dat['ml_s'], func=np.nanmean, window=6).data)
dat['ml_t_smooth'] = (('time'), gt.cleaning.rolling_window(dat['ml_t'], func=np.nanmean, window=6).data)
era5 = xr.open_mfdataset('../data/ERA5_flux_data/*.nc')
era5 = era5.sel(latitude=slice(-40, -65), longitude=slice(-5, 10))
var_name=['slhf', 'sshf', 'str', 'ssr', 'tp', 'e']
era5=my.convert_era5_to_Wm2(ds=era5, var_name=var_name)
era5_interp = era5.resample(time=interp_res, loffset=offset).mean()
era5['emp'] = (('time', 'latitude', 'longitude'), (era5['e']-era5['tp']).data)
era5['qnet'] = (('time', 'latitude', 'longitude'), (era5['slhf']+era5['sshf']+era5['str' ]+era5['ssr' ]).data)
from tqdm import tqdm
var_name = [
'slhf',
'sshf',
'str' ,
'ssr' ,
'tp' ,
'e' ,
# 'emp',
'u10',
'v10',
# 'qnet'
]
for var in tqdm(var_name):
dat_saz_interp[var] = (('time'), my.interp_glider_era5(dat_saz_interp, era5_interp, var=var))
dat_pfz_interp[var] = (('time'), my.interp_glider_era5(dat_pfz_interp, era5_interp, var=var))
dat_miz_interp[var] = (('time'), my.interp_glider_era5(dat_miz_interp, era5_interp, var=var))
for i, dat in enumerate([dat_saz_interp, dat_pfz_interp, dat_miz_interp]):
dat['emp'] = (('time'), (dat['e']-dat['tp']).data)
dat['qnet'] = (('time'), (dat['slhf']+dat['sshf']+dat['str' ]+dat['ssr' ]).data)
100%|████████████████████████████████████████████| 8/8 [17:09<00:00, 128.63s/it]
def getWindDir(uwind, vwind):
data_dir = []
for i in range(len(uwind)):
data_dir += math.atan2(vwind[i], uwind[i])/math.pi*180,
data_dir = np.array(data_dir)
return data_dir
def setWindDirZero360(wind_dir):
ind1 = np.nonzero((wind_dir > 0) & (wind_dir < 90))
ind2 = np.nonzero((wind_dir > 90) & (wind_dir < 190))
ind3 = np.nonzero((wind_dir > -180) & (wind_dir < -90))
ind4 = np.nonzero((wind_dir > -90) & (wind_dir < 0))
wind_dir[ind1] = 90 - wind_dir[ind1]
wind_dir[ind2] = 450 - wind_dir[ind2]
wind_dir[ind3] = abs(wind_dir[ind3] - 90)
wind_dir[ind4] = abs(wind_dir[ind4]) + 90
return wind_dir
import windstress
import math
airT = [14, 0, -5]
for i, dat in enumerate([dat_saz_interp, dat_pfz_interp, dat_miz_interp]):
dat['taux'] = (('time'), windstress.stress(dat.u10, z=10., drag='largepond', rho_air=1.22, Ta=airT[i]))
dat['tauy'] = (('time'), windstress.stress(dat.v10, z=10., drag='largepond', rho_air=1.22, Ta=airT[i]))
dat['wind_dir'] = (('time'), getWindDir(dat.u10, dat.v10).data)
dat['wind_dir'] = (('time'), setWindDirZero360(dat['wind_dir'].values).data)
dat_saz_interp.to_netcdf('../data/dat_saz_6H.nc')
dat_pfz_interp.to_netcdf('../data/dat_pfz_6H.nc')
dat_miz_interp.to_netcdf('../data/dat_miz_6H.nc')