"""
example of running AirSeaFluxCode with
1. R/V data (data_all.csv) or
2. one day era5 hourly data (era5_r360x180.nc)
compute fluxes
output NetCDF4
and statistics in stats.txt
@author: sbiri
"""
#%% import packages
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import pandas as pd
from AirSeaFluxCode import AirSeaFluxCode
import time
from tabulate import tabulate
#%%
def reject_outliers(data, m=2):
x = np.copy(data)
x = np.where(np.abs(x - np.nanmean(x)) < m*np.nanstd(x),
x, np.nan)
return x
def toy_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
"""
Parameters
----------
inF : str
input filename either data_all.csv or era5_r360x180.nc
outF : str
output filename
gustIn : float
gustiness option e.g. [1, 1.2, 800]
cskinIn : int
cool skin option input 0 or 1
tolIn : float
tolerance input option e.g. ['all', 0.01, 0.01, 5e-05, 0.01, 1, 1]
meth : str
parametrisation method option
Returns
-------
res : float
AirSeaFluxCode output
lon : float
longitude from input netCDF file
lat : float
latitude from input netCDF file
"""
if (inF == "data_all.csv"):
#%% load data_all
inDt = pd.read_csv("data_all.csv")
date = np.asarray(inDt["Date"])
lon = np.asarray(inDt["Longitude"])
lat = np.asarray(inDt["Latitude"])
spd = np.asarray(inDt["Wind speed"])
t = np.asarray(inDt["Air temperature"])
sst = np.asarray(inDt["SST"])
rh = np.asarray(inDt["RH"])
p = np.asarray(inDt["P"])
sw = np.asarray(inDt["Rs"])
hu = np.asarray(inDt["zu"])
ht = np.asarray(inDt["zt"])
hin = np.array([hu, ht, ht])
del hu, ht, inDt
#%% run AirSeaFluxCode
temp = AirSeaFluxCode(spd, t, sst, lat=lat, hum=['rh', rh], P=p,
hin=hin, Rs=sw, tol=tolIn, gust=gustIn,
cskin=cskinIn, meth=meth, L="ecmwf", n=30)
res = temp.loc[:,"tau":"Rnl"]
res = res.to_numpy().T
del temp
#%% delete variables
del spd, t, sst, rh, p, sw, hin
elif (inF == 'era5_r360x180.nc'):
#%% load era5_r360x180.nc
fid = nc.Dataset(inF)
lon = np.array(fid.variables["lon"])
lat = np.array(fid.variables["lat"])
T = np.array(fid.variables["t2m"])
tim = np.array(fid.variables["time"])
Td = np.array(fid.variables["d2m"])
sst = np.array(fid.variables["sst"])
sst = np.where(sst < -100, np.nan, sst)
p = np.array(fid.variables["msl"])/100 # to set hPa
lw = np.array(fid.variables["strd"])/60/60
sw = np.array(fid.variables["ssrd"])/60/60
u = np.array(fid.variables["u10"])
v = np.array(fid.variables["v10"])
lsm = np.array(fid.variables["lsm"])
fid.close()
spd = np.sqrt(np.power(u, 2)+np.power(v, 2))
del u, v, fid
lsm = np.where(lsm > 0, np.nan, 1) # reverse 0 on land 1 over ocean
hin = np.array([10, 2, 2])
latIn = np.tile(lat, (len(lon), 1)).T.reshape(len(lon)*len(lat))
date = np.copy(tim)
#%% run AirSeaFluxCode
res = np.zeros((len(tim),len(lon)*len(lat), 36))
# reshape input and run code
for x in range(len(tim)):
temp = AirSeaFluxCode(spd.reshape(len(tim), len(lon)*len(lat))[x, :],
T.reshape(len(tim), len(lon)*len(lat))[x, :],
sst.reshape(len(tim), len(lon)*len(lat))[x, :],
lat=latIn,
hum=['Td', Td.reshape(len(tim), len(lon)*len(lat))[x, :]],
P=p.reshape(len(tim), len(lon)*len(lat))[x, :],
hin=hin,
Rs=sw.reshape(len(tim), len(lon)*len(lat))[x, :],
Rl=lw.reshape(len(tim), len(lon)*len(lat))[x, :],
gust=gustIn, cskin=cskinIn, tol=tolIn, qmeth='WMO',
meth=meth, n=30, L="ecmwf")
a = temp.loc[:,"tau":"Rnl"]
a = a.to_numpy()
del temp
res[x, :, :] = a
del a
if (outF[-3:] == '.nc'):
if (inF == 'era5_r360x180.nc'):
#%% save NetCDF4
fid = nc.Dataset(outF,'w', format='NETCDF4')
fid.createDimension('lon', len(lon))
fid.createDimension('lat', len(lat))
fid.createDimension('time', None)
longitude = fid.createVariable('lon', 'f4', 'lon')
latitude = fid.createVariable('lat', 'f4', 'lat')
Date = fid.createVariable('Date', 'i4', 'time')
tau = fid.createVariable('tau', 'f4', ('time','lat','lon'))
sensible = fid.createVariable('shf', 'f4', ('time','lat','lon'))
latent = fid.createVariable('lhf', 'f4', ('time','lat','lon'))
monob = fid.createVariable('MO', 'f4', ('time','lat','lon'))
cd = fid.createVariable('cd', 'f4', ('time','lat','lon'))
cdn = fid.createVariable('cdn', 'f4', ('time','lat','lon'))
ct = fid.createVariable('ct', 'f4', ('time','lat','lon'))
ctn = fid.createVariable('ctn', 'f4', ('time','lat','lon'))
cq = fid.createVariable('cq', 'f4', ('time','lat','lon'))
cqn = fid.createVariable('cqn', 'f4', ('time','lat','lon'))
tsrv = fid.createVariable('tsrv', 'f4', ('time','lat','lon'))
tsr = fid.createVariable('tsr', 'f4', ('time','lat','lon'))
qsr = fid.createVariable('qsr', 'f4', ('time','lat','lon'))
usr = fid.createVariable('usr', 'f4', ('time','lat','lon'))
psim = fid.createVariable('psim', 'f4', ('time','lat','lon'))
psit = fid.createVariable('psit', 'f4', ('time','lat','lon'))
psiq = fid.createVariable('psiq', 'f4', ('time','lat','lon'))
u10n = fid.createVariable('u10n', 'f4', ('time','lat','lon'))
t10n = fid.createVariable('t10n', 'f4', ('time','lat','lon'))
tv10n = fid.createVariable('tv10n', 'f4', ('time','lat','lon'))
q10n = fid.createVariable('q10n', 'f4', ('time','lat','lon'))
zo = fid.createVariable('zo', 'f4', ('time','lat','lon'))
zot = fid.createVariable('zot', 'f4', ('time','lat','lon'))
zoq = fid.createVariable('zoq', 'f4', ('time','lat','lon'))
urefs = fid.createVariable('uref', 'f4', ('time','lat','lon'))
trefs = fid.createVariable('tref', 'f4', ('time','lat','lon'))
qrefs = fid.createVariable('qref', 'f4', ('time','lat','lon'))
itera = fid.createVariable('iter', 'i4', ('time','lat','lon'))
dter = fid.createVariable('dter', 'f4', ('time','lat','lon'))
dqer = fid.createVariable('dqer', 'f4', ('time','lat','lon'))
dtwl = fid.createVariable('dtwl', 'f4', ('time','lat','lon'))
qair = fid.createVariable('qair', 'f4', ('time','lat','lon'))
qsea = fid.createVariable('qsea', 'f4', ('time','lat','lon'))
Rl = fid.createVariable('Rl', 'f4', ('time','lat','lon'))
Rs = fid.createVariable('Rs', 'f4', ('time','lat','lon'))
Rnl = fid.createVariable('Rnl', 'f4', ('time','lat','lon'))
longitude[:] = lon
latitude[:] = lat
Date[:] = tim
tau[:] = res[:, :, 0].reshape((len(tim), len(lat), len(lon)))*lsm
sensible[:] = res[:, :, 1].reshape((len(tim), len(lat), len(lon)))*lsm
latent[:] = res[:, :, 2].reshape((len(tim), len(lat), len(lon)))*lsm
monob[:] = res[:, :, 3].reshape((len(tim), len(lat), len(lon)))*lsm
cd[:] = res[:, :, 4].reshape((len(tim), len(lat), len(lon)))*lsm
cdn[:] = res[:, :, 5].reshape((len(tim), len(lat), len(lon)))*lsm
ct[:] = res[:, :, 6].reshape((len(tim), len(lat), len(lon)))*lsm
ctn[:] = res[:, :, 7].reshape((len(tim), len(lat), len(lon)))*lsm
cq[:] = res[:, :, 8].reshape((len(tim), len(lat), len(lon)))*lsm
cqn[:] = res[:, :, 9].reshape((len(tim), len(lat), len(lon)))*lsm
tsrv[:] = res[:, :, 10].reshape((len(tim), len(lat), len(lon)))*lsm
tsr[:] = res[:, :, 11].reshape((len(tim), len(lat), len(lon)))*lsm
qsr[:] = res[:, :, 12].reshape((len(tim), len(lat), len(lon)))*lsm
usr[:] = res[:, :, 13].reshape((len(tim), len(lat), len(lon)))*lsm
psim[:] = res[:, :, 14].reshape((len(tim), len(lat), len(lon)))*lsm
psit[:] = res[:, :, 15].reshape((len(tim), len(lat), len(lon)))*lsm
psiq[:] = res[:, :, 16].reshape((len(tim), len(lat), len(lon)))*lsm
u10n[:] = res[:, :, 17].reshape((len(tim), len(lat), len(lon)))*lsm
t10n[:] = res[:, :, 18].reshape((len(tim), len(lat), len(lon)))*lsm
tv10n[:] = res[:, :, 19].reshape((len(tim), len(lat), len(lon)))*lsm
q10n[:] = res[:, :, 20].reshape((len(tim), len(lat), len(lon)))*lsm
zo[:] = res[:, :, 21].reshape((len(tim), len(lat), len(lon)))*lsm
zot[:] = res[:, :, 22].reshape((len(tim), len(lat), len(lon)))*lsm
zoq[:] = res[:, :, 23].reshape((len(tim), len(lat), len(lon)))*lsm
urefs[:] = res[:, :, 24].reshape((len(tim), len(lat), len(lon)))*lsm
trefs[:] = res[:, :, 25].reshape((len(tim), len(lat), len(lon)))*lsm
qrefs[:] = res[:, :, 26].reshape((len(tim), len(lat), len(lon)))*lsm
itera[:] = res[:, :, 27].reshape((len(tim), len(lat), len(lon)))*lsm
dter[:] = res[:, :, 28].reshape((len(tim), len(lat), len(lon)))*lsm
dqer[:] = res[:, :, 29].reshape((len(tim), len(lat), len(lon)))*lsm
dtwl[:] = res[:, :, 30].reshape((len(tim), len(lat), len(lon)))*lsm
qair[:] = res[:, :, 31].reshape((len(tim), len(lat), len(lon)))*lsm
qsea[:] = res[:, :, 32].reshape((len(tim), len(lat), len(lon)))*lsm
Rl = res[:, :, 33].reshape((len(tim), len(lat), len(lon)))
Rs = res[:, :, 34].reshape((len(tim), len(lat), len(lon)))
Rnl = res[:, :, 35].reshape((len(tim), len(lat), len(lon)))
longitude.long_name = 'Longitude'
longitude.units = 'degrees East'
latitude.long_name = 'Latitude'
latitude.units = 'degrees North'
Date.long_name = "calendar date"
Date.units = "YYYYMMDD UTC"
tau.long_name = 'Wind stress'
tau.units = 'N/m^2'
sensible.long_name = 'Sensible heat fluxe'
sensible.units = 'W/m^2'
latent.long_name = 'Latent heat flux'
latent.units = 'W/m^2'
monob.long_name = 'Monin-Obukhov length'
monob.units = 'm'
cd.long_name = 'Drag coefficient'
cd.units = ''
cdn.long_name = 'Neutral Drag coefficient'
cdn.units = ''
ct.long_name = 'Heat exchange coefficient'
ct.units = ''
ctn.long_name = 'Neutral Heat exchange coefficient'
ctn.units = ''
cq.long_name = 'Moisture exchange coefficient'
cq.units = ''
cqn.long_name = 'Neutral Moisture exchange coefficient'
cqn.units = ''
tsrv.long_name = 'star virtual temperature'
tsrv.units = 'degrees Celsius'
tsr.long_name = 'star temperature'
tsr.units = 'degrees Celsius'
qsr.long_name = 'star specific humidity'
qsr.units = 'gr/kgr'
usr.long_name = 'friction velocity'
usr.units = 'm/s'
psim.long_name = 'Momentum stability function'
psit.long_name = 'Heat stability function'
u10n.long_name = '10m neutral wind speed'
u10n.units = 'm/s'
t10n.long_name = '10m neutral temperature'
t10n.units = 'degrees Celsius'
tv10n.long_name = '10m neutral virtual temperature'
tv10n.units = 'degrees Celsius'
q10n.long_name = '10m neutral specific humidity'
q10n.units = 'gr/kgr'
zo.long_name = 'momentum roughness length'
zo.units = 'm'
zot.long_name = 'temperature roughness length'
zot.units = 'm'
zoq.long_name = 'moisture roughness length'
zoq.units = 'm'
urefs.long_name = 'wind speed at ref height'
urefs.units = 'm/s'
trefs.long_name = 'temperature at ref height'
trefs.units = 'degrees Celsius'
qrefs.long_name = 'specific humidity at ref height'
qrefs.units = 'gr/kgr'
qair.long_name = 'specific humidity of air'
qair.units = 'gr/kgr'
qsea.long_name = 'specific humidity over water'
qsea.units = 'gr/kgr'
itera.long_name = 'number of iterations'
fid.close()
#%% delete variables
del longitude, latitude, Date, tau, sensible, latent, monob, cd, cdn
del ct, ctn, cq, cqn, tsrv, tsr, qsr, usr, psim, psit, psiq, u10n, t10n
del tv10n, q10n, zo, zot, zoq, urefs, trefs, qrefs, itera, dter, dqer
del qair, qsea, Rl, Rs, Rnl, dtwl
del tim, T, Td, p, lw, sw, lsm, spd, hin, latIn
else:
#%% save NetCDF4
fid = nc.Dataset(outF,'w', format='NETCDF4')
fid.createDimension('lon', len(lon))
fid.createDimension('lat', len(lat))
fid.createDimension('time', None)
longitude = fid.createVariable('lon', 'f4', 'lon')
latitude = fid.createVariable('lat', 'f4', 'lat')
Date = fid.createVariable('Date', 'i4', 'time')
tau = fid.createVariable('tau', 'f4', 'time')
sensible = fid.createVariable('shf', 'f4', 'time')
latent = fid.createVariable('lhf', 'f4', 'time')
monob = fid.createVariable('MO', 'f4', 'time')
cd = fid.createVariable('cd', 'f4', 'time')
cdn = fid.createVariable('cdn', 'f4', 'time')
ct = fid.createVariable('ct', 'f4', 'time')
ctn = fid.createVariable('ctn', 'f4', 'time')
cq = fid.createVariable('cq', 'f4', 'time')
cqn = fid.createVariable('cqn', 'f4', 'time')
tsrv = fid.createVariable('tsrv', 'f4', 'time')
tsr = fid.createVariable('tsr', 'f4', 'time')
qsr = fid.createVariable('qsr', 'f4', 'time')
usr = fid.createVariable('usr', 'f4', 'time')
psim = fid.createVariable('psim', 'f4', 'time')
psit = fid.createVariable('psit', 'f4', 'time')
psiq = fid.createVariable('psiq', 'f4', 'time')
u10n = fid.createVariable('u10n', 'f4', 'time')
t10n = fid.createVariable('t10n', 'f4', 'time')
tv10n = fid.createVariable('tv10n', 'f4', 'time')
q10n = fid.createVariable('q10n', 'f4', 'time')
zo = fid.createVariable('zo', 'f4', 'time')
zot = fid.createVariable('zot', 'f4', 'time')
zoq = fid.createVariable('zoq', 'f4', 'time')
urefs = fid.createVariable('uref', 'f4', 'time')
trefs = fid.createVariable('tref', 'f4', 'time')
qrefs = fid.createVariable('qref', 'f4', 'time')
itera = fid.createVariable('iter', 'i4', 'time')
dter = fid.createVariable('dter', 'f4', 'time')
dqer = fid.createVariable('dqer', 'f4', 'time')
dtwl = fid.createVariable('dter', 'f4', 'time')
qair = fid.createVariable('qair', 'f4', 'time')
qsea = fid.createVariable('qsea', 'f4', 'time')
Rl = fid.createVariable('Rl', 'f4', 'time')
Rs = fid.createVariable('Rs', 'f4', 'time')
Rnl = fid.createVariable('Rnl', 'f4', 'time')
longitude[:] = lon
latitude[:] = lat
Date[:] = date
tau[:] = res[0]
sensible[:] = res[1]
latent[:] = res[2]
monob[:] = res[3]
cd[:] = res[4]
cdn[:] = res[5]
ct[:] = res[6]
ctn[:] = res[7]
cq[:] = res[8]
cqn[:] = res[9]
tsrv[:] = res[10]
tsr[:] = res[11]
qsr[:] = res[12]
usr[:] = res[13]
psim[:] = res[14]
psit[:] = res[15]
psiq[:] = res[16]
u10n[:] = res[17]
t10n[:] = res[18]
tv10n[:] = res[19]
q10n[:] = res[20]
zo[:] = res[21]
zot[:] = res[22]
zoq[:] = res[23]
urefs[:] = res[24]
trefs[:] = res[25]
qrefs[:] = res[26]
itera[:] = res[27]
dter[:] = res[28]
dqer[:] = res[29]
dtwl[:] = res[30]
qair[:] = res[31]
qsea[:] = res[32]
Rl[:] = res[33]
Rs[:] = res[34]
Rnl[:] = res[35]
longitude.long_name = 'Longitude'
longitude.units = 'degrees East'
latitude.long_name = 'Latitude'
latitude.units = 'degrees North'
Date.long_name = "calendar date"
Date.units = "YYYYMMDD UTC"
tau.long_name = 'Wind stress'
tau.units = 'N/m^2'
sensible.long_name = 'Sensible heat fluxe'
sensible.units = 'W/m^2'
latent.long_name = 'Latent heat flux'
latent.units = 'W/m^2'
monob.long_name = 'Monin-Obukhov length'
monob.units = 'm'
cd.long_name = 'Drag coefficient'
cd.units = ''
cdn.long_name = 'Neutral Drag coefficient'
cdn.units = ''
ct.long_name = 'Heat exchange coefficient'
ct.units = ''
ctn.long_name = 'Neutral Heat exchange coefficient'
ctn.units = ''
cq.long_name = 'Moisture exchange coefficient'
cq.units = ''
cqn.long_name = 'Neutral Moisture exchange coefficient'
cqn.units = ''
tsrv.long_name = 'star virtual temperature'
tsrv.units = 'degrees Celsius'
tsr.long_name = 'star temperature'
tsr.units = 'degrees Celsius'
qsr.long_name = 'star specific humidity'
qsr.units = 'gr/kgr'
usr.long_name = 'friction velocity'
usr.units = 'm/s'
psim.long_name = 'Momentum stability function'
psit.long_name = 'Heat stability function'
u10n.long_name = '10m neutral wind speed'
u10n.units = 'm/s'
t10n.long_name = '10m neutral temperature'
t10n.units = 'degrees Celsius'
tv10n.long_name = '10m neutral virtual temperature'
tv10n.units = 'degrees Celsius'
q10n.long_name = '10m neutral specific humidity'
q10n.units = 'gr/kgr'
zo.long_name = 'momentum roughness length'
zo.units = 'm'
zot.long_name = 'temperature roughness length'
zot.units = 'm'
zoq.long_name = 'moisture roughness length'
zoq.units = 'm'
urefs.long_name = 'wind speed at ref height'
urefs.units = 'm/s'
trefs.long_name = 'temperature at ref height'
trefs.units = 'degrees Celsius'
qrefs.long_name = 'specific humidity at ref height'
qrefs.units = 'gr/kgr'
qair.long_name = 'specific humidity of air'
qair.units = 'gr/kgr'
qsea.long_name = 'specific humidity over water'
qsea.units = 'gr/kgr'
itera.long_name = 'number of iterations'
fid.close()
#%% delete variables
del longitude, latitude, Date, tau, sensible, latent, monob, cd, cdn
del ct, ctn, cq, cqn, tsrv, tsr, qsr, usr, psim, psit, psiq, u10n, t10n
del tv10n, q10n, zo, zot, zoq, urefs, trefs, qrefs, itera, dter, dqer
del qair, qsea, Rl, Rs, Rnl
del t, rh, date, p, sw, spd, hin
else:
#%% save as .csv
np.savetxt(outF, np.vstack((date, lon, lat, res)).T,
delimiter=',',
header="date, lon, lat, tau, shf, lhf, L, cd, cdn, ct, ctn,"
" cq, cqn, tsrv, tsr, qsr, usr, psim, psit, psiq, u10n,"
" t10n, tv10n, q10n, zo, zot, zoq, uref, tref, qref, iter,"
" dter, dqer, dtwl, qair, qsea, Rl, Rs, Rnl")
return res, lon, lat
#%% run function
start_time = time.perf_counter()
#------------------------------------------------------------------------------
inF = input("Give input file name (data_all.csv or era5_r360x180.nc): \n")
meth = input("Give prefered method: \n")
while meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
"C40", "ecmwf","Beljaars"]:
print("method unknown")
meth = input("Give prefered method: \n")
else:
meth = meth #[meth]
ext = meth+"_"
#------------------------------------------------------------------------------
gustIn = input("Give gustiness option (to use default press enter): \n")
if (gustIn == ''):
gustIn = None
ext = ext+'gust_'
else:
gustIn = np.asarray(eval(gustIn), dtype=float)
if ((np.all(gustIn) == 0)):
ext = ext+'nogust_'
else:
ext = ext+'gust_'
#------------------------------------------------------------------------------
cskinIn = input("Give cool skin option (to use default press enter): \n")
if (cskinIn == ''):
cskinIn = None
if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96" or meth == "UA" or
meth == "LY04")):
cskinIn = 0
ext = ext+'noskin_'
elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
or meth == "ecmwf" or meth == "Beljaars")):
cskinIn = 1
ext = ext+'skin_'
else:
cskinIn = int(cskinIn)
if (cskinIn == 0):
ext = ext+'noskin_'
elif (cskinIn == 1):
ext = ext+'skin_'
#------------------------------------------------------------------------------
tolIn = input("Give tolerance option (to use default press enter): \n")
if (tolIn == ''):
tolIn = ['flux', 1e-3, 0.1, 0.1]
else:
tolIn = eval(tolIn)
ext = ext+'tol'+tolIn[0]
#------------------------------------------------------------------------------
outF = input("Give path and output file name: \n")
if ((outF == '') and (inF == "data_all.csv")):
outF = "out_"+inF[:-4]+"_"+ext+".csv"
elif ((outF == '') and (inF == "era5_r360x180.nc")):
outF = "out_"+inF[:-3]+"_"+ext+".nc"
elif ((outF[-4:] == '.csv') and (inF == 'era5_r360x180.nc')):
outF = outF[:-4]+".nc"
elif ((outF[-3:] != '.nc') and (outF[-4:] != '.csv')):
if (inF == "data_all.csv"):
outF = outF+".csv"
else:
outF = outF+".nc"
else:
outF = outF
#------------------------------------------------------------------------------
print("\n run_ASFC.py, started for method "+meth)
res, lon, lat = toy_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth)
print("run_ASFC.py took ", np.round((time.perf_counter()-start_time)/60, 2),
"minutes to run")
#%% generate flux plots
if (inF == 'era5_r360x180.nc'):
cm = plt.cm.get_cmap('RdYlBu')
ttl = ["tau (Nm$^{-2}$)", "shf (Wm$^{-2}$)", "lhf (Wm$^{-2}$)"]
for i in range(3):
plt.figure()
plt.contourf(lon, lat,
np.nanmean(res[:, :, i], axis=0).reshape(len(lat),
len(lon)),
100, cmap=cm)
plt.colorbar()
plt.tight_layout()
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title(meth+', '+ttl[i])
plt.savefig('./'+ttl[i][:3]+'_'+ext+'.png', dpi=300, bbox_inches='tight')
elif (inF == "data_all.csv"):
ttl = ["tau (Nm$^{-2}$)", "shf (Wm$^{-2}$)", "lhf (Wm$^{-2}$)"]
for i in range(3):
plt.figure()
plt.plot(res[i],'.c', markersize=1)
plt.title(meth)
plt.xlabel("points")
plt.ylabel(ttl[i])
plt.savefig('./'+ttl[i][:3]+'_'+ext+'.png', dpi=300, bbox_inches='tight')
#%% generate txt file with statistic
if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
or meth == "YT96" or meth == "UA" or
meth == "LY04")):
cskinIn = 0
elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
or meth == "ecmwf" or meth == "Beljaars")):
cskinIn = 1
if (np.all(gustIn == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
gustIn = [1, 1.2, 600]
elif (np.all(gustIn == None) and (meth == "UA" or meth == "ecmwf")):
gustIn = [1, 1, 1000]
elif np.all(gustIn == None):
gustIn = [1, 1.2, 800]
elif ((np.size(gustIn) < 3) and (gustIn == 0)):
gust = [0, 0, 0]
if (tolIn == None):
tolIn = ['flux', 0.01, 1, 1]
print("Input summary", file=open('./stats.txt', 'a'))
print('input file name: {}, \n method: {}, \n gustiness: {}, \n cskin: {},'
' \n tolerance: {}'.format(inF, meth, gustIn, cskinIn, tolIn),
file=open('./stats.txt', 'a'))
ttl = np.asarray(["tau ", "shf ", "lhf ", "L ", "cd ", "cdn ",
"ct ", "ctn ", "cq ", "cqn ", "tsrv ", "tsr ",
"qsr ", "usr ", "psim ", "psit ", "psiq ", "u10n ",
"t10n ", "tv10n", "q10n ", "zo ", "zot ", "zoq ",
"urefs", "trefs", "qrefs", "itera", "dter ", "dqer ",
"dtwl ", "qair ", "qsea ", "Rl ", "Rs ", "Rnl "])
header = ["var", "mean", "median", "min", "max", "5%", "95%"]
n = np.shape(res)
stats = np.copy(ttl)
if (inF == 'era5_r360x180.nc'):
stats = np.c_[stats, np.nanmean(res.reshape(n[0]*n[1], n[2]), axis=0)]
stats = np.c_[stats, np.nanmedian(res.reshape(n[0]*n[1], n[2]), axis=0)]
stats = np.c_[stats, np.nanmin(res.reshape(n[0]*n[1], n[2]), axis=0)]
stats = np.c_[stats, np.nanmax(res.reshape(n[0]*n[1], n[2]), axis=0)]
stats = np.c_[stats, np.nanpercentile(res.reshape(n[0]*n[1], n[2]), 5,
axis=0)]
stats = np.c_[stats, np.nanpercentile(res.reshape(n[0]*n[1], n[2]), 95,
axis=0)]
print(tabulate(stats, headers=header, tablefmt="github", numalign="left",
floatfmt=("s", "2.2e", "2.2e", "2.2e", "2.2e", "2.2e",
"2.2e")), file=open('./stats.txt', 'a'))
print('-'*79+'\n', file=open('./stats.txt', 'a'))
elif (inF == "data_all.csv"):
stats = np.c_[stats, np.nanmean(res, axis=1)]
stats = np.c_[stats, np.nanmedian(res, axis=1)]
stats = np.c_[stats, np.nanmin(res, axis=1)]
stats = np.c_[stats, np.nanmax(res, axis=1)]
stats = np.c_[stats, np.nanpercentile(res, 5, axis=1)]
stats = np.c_[stats, np.nanpercentile(res, 95, axis=1)]
print(tabulate(stats, headers=header, tablefmt="github", numalign="left",
floatfmt=("s", "2.2e", "2.2e", "2.2e", "2.2e", "2.2e",
"2.2e")), file=open('./stats.txt', 'a'))
print('-'*79+'\n', file=open('./stats.txt', 'a'))
print('input file name: {}, \n method: {}, \n gustiness: {}, \n cskin: {},'
' \n tolerance: {}, \n output is written in: {}'.format(inF, meth,
gustIn, cskinIn,
tolIn, outF),
file=open('./readme.txt', 'w'))