#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
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 or csv
and statistics in "outS".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 toy_ASFC(inF, outF, outS, sst_fl, gustIn, cskinIn, tolIn, meth, qmIn, LIn,
stdIn):
"""
Example routine of how to run AirSeaFluxCode with the test data given
and save output either as .csv or NetCDF
Parameters
----------
inF : str
input filename either data_all.csv or era5_r360x180.nc
outF : str
output filename
outS : str
output statistics 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, 1e-05, 1e-3, 0.1, 0.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
"""
out_var = ("tau", "sensible", "latent", "u10n", "t10n", "q10n")
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"])
spd = spd+np.random.normal(0, stdIn[0], spd.shape)
t = np.asarray(inDt["Air temperature"])
t = t+np.random.normal(0, stdIn[1], t.shape)
sst = np.asarray(inDt["SST"])
sst = sst+np.random.normal(0, stdIn[2], sst.shape)
rh = np.asarray(inDt["RH"])
rh = rh+np.random.normal(0, stdIn[3], rh.shape)
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
res = AirSeaFluxCode(spd, t, sst, sst_fl, lat=lat, hum=['rh', rh], P=p,
hin=hin, Rs=sw, tol=tolIn, gust=gustIn,
cskin=cskinIn, meth=meth, qmeth=qmIn, L=LIn,
maxiter=10, out_var = out_var)
flg = res["flag"]
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"])
tim = np.array(fid.variables["time"])
lsm = np.array(fid.variables["lsm"])
icon = np.array(fid.variables["siconc"])
lsm = np.where(lsm > 0, np.nan, 1) # reverse 0 on land 1 over ocean
icon = np.where(icon == 0, 1, np.nan) # keep only ice-free regions
msk = lsm*icon
T = np.array(fid.variables["t2m"])*msk
T = T+np.random.normal(0, stdIn[1], T.shape)
Td = np.array(fid.variables["d2m"])*msk
Td = Td+np.random.normal(0, stdIn[3], T.shape)
sst = np.array(fid.variables["sst"])
sst = np.where(sst < -100, np.nan, sst)*msk
sst = sst+np.random.normal(0, stdIn[2], sst.shape)
p = np.array(fid.variables["msl"])*msk/100 # to set hPa
lw = np.array(fid.variables["strd"])*msk/60/60
sw = np.array(fid.variables["ssrd"])*msk/60/60
u = np.array(fid.variables["u10"])
v = np.array(fid.variables["v10"])
fid.close()
spd = np.sqrt(np.power(u, 2)+np.power(v, 2))*msk
spd = spd+np.random.normal(0, stdIn[0], spd.shape)
del u, v, fid
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), 6))
flg = np.empty((len(tim),len(lon)*len(lat)), dtype="object")
# reshape input and run code
for x in range(len(tim)):
a = 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, :],
sst_fl, 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,
meth=meth, qmeth=qmIn, maxiter=30, L=LIn,
out_var = out_var)
temp = a.iloc[:, :-1]
temp = temp.to_numpy()
flg[x, :] = a["flag"]
res[x, :, :] = temp
del a, temp
n = np.shape(res)
res = np.asarray([res[:, :, i]*msk.reshape(n[0], n[1])
for i in range(6)])
res = np.moveaxis(res, 0, -1)
flg = np.where(np.isnan(msk.reshape(len(tim), len(lon)*len(lat))),
'm', flg)
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'))
u10n = fid.createVariable('u10n', 'f4', ('time','lat','lon'))
t10n = fid.createVariable('t10n', 'f4', ('time','lat','lon'))
q10n = fid.createVariable('q10n', 'f4', ('time','lat','lon'))
flag = fid.createVariable('flag', 'U1', ('time','lat','lon'))
longitude[:] = lon
latitude[:] = lat
Date[:] = tim
tau[:] = res[:, :, 0].reshape((len(tim), len(lat), len(lon)))*msk
sensible[:] = res[:, :, 1].reshape((len(tim), len(lat), len(lon)))*msk
latent[:] = res[:, :, 2].reshape((len(tim), len(lat), len(lon)))*msk
u10n[:] = res[:, :, 3].reshape((len(tim), len(lat), len(lon)))*msk
t10n[:] = res[:, :, 4].reshape((len(tim), len(lat), len(lon)))*msk
q10n[:] = res[:, :, 5].reshape((len(tim), len(lat), len(lon)))*msk
flag[:] = flg.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 = "gregorian"
Date.units = "hours since 1900-01-01 00:00:00.0"
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'
u10n.long_name = '10m neutral wind speed'
u10n.units = 'm/s'
t10n.long_name = '10m neutral temperature'
t10n.units = 'degrees Celsius'
q10n.long_name = '10m neutral specific humidity'
q10n.units = 'kgr/kgr'
flag.long_name = ('flag "n" normal, "u": u10n < 0, "q": q10n < 0,'
'"l": zol<0.01, "m": missing, "i": points that'
'have not converged')
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, icon, msk
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')
u10n = fid.createVariable('u10n', 'f4', 'time')
t10n = fid.createVariable('t10n', 'f4', 'time')
q10n = fid.createVariable('q10n', 'f4', 'time')
flag = fid.createVariable('flag', 'U1', 'time')
longitude[:] = lon
latitude[:] = lat
Date[:] = date
tau[:] = res["tau"]
sensible[:] = res["shf"]
latent[:] = res["lhf"]
u10n[:] = res["u10n"]
t10n[:] = res["t10n"]
q10n[:] = res["q10n"]
flag[:] = res["flag"]
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'
u10n.long_name = '10m neutral wind speed'
u10n.units = 'm/s'
t10n.long_name = '10m neutral temperature'
t10n.units = 'degrees Celsius'
q10n.long_name = '10m neutral specific humidity'
q10n.units = 'kgr/kgr'
flag.long_name = ('flag "n" normal, "u": u10n < 0, "q": q10n < 0,'
'"l": zol<0.01, "m": missing, "i": points that'
'have not converged')
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, ug, rh, Rib
# del t, date, p, sw, spd, hin, sst
else:
#%% save as .csv
res.insert(loc=0, column='date', value=date)
res.insert(loc=1, column='lon', value=lon)
res.insert(loc=2, column='lat', value=lat)
res.to_csv(outF, float_format='%.3f')
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", "NCAR", "C30", "C35",
"ecmwf","Beljaars"]:
print("method unknown")
meth = input("Give prefered method: \n")
else:
meth = meth #[meth]
ext = meth+"_"
#------------------------------------------------------------------------------
sst_fl = input("Give SST flag (bulk/sin): \n")
#------------------------------------------------------------------------------
qmIn = input("Give prefered method for specific humidity: \n")
if (qmIn == ''):
qmIn = 'Buck2' # default
while qmIn not in ["HylandWexler", "Hardy", "Preining", "Wexler",
"GoffGratch", "WMO", "MagnusTetens", "Buck", "Buck2",
"WMO2018", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
print("method unknown")
meth = input("Give prefered method: \n")
else:
qmIn = qmIn
#------------------------------------------------------------------------------
gustIn = input("Give gustiness option (to switch it off enter 0;"
" to set your own input use the form [1, B, zi, ugmin]"
" i.e. [1, 1, 800, 0.01] or "
"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 == "NCAR")):
cskinIn = 0
ext = ext+'noskin_'
elif ((cskinIn == None) and (meth == "C30" or meth == "C35"
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 = ['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1]
else:
tolIn = eval(tolIn)
ext = ext+'tol'+tolIn[0]
#------------------------------------------------------------------------------
LIn = input("Give prefered method for L (tsrv or Rb): \n")
if (LIn == ''):
LIn = 'tsrv' # default
elif LIn not in ["tsrv", "Rb"]:
LIn = 'tsrv' # default
else:
LIn = LIn
#------------------------------------------------------------------------------
stdIn = input("Give noise std for spd, T, SST and Td/rh \n (e.g. [0.01, 0, 0, 0]"
" adds noise only to spd): \n")
if (stdIn == ''):
stdIn = [0, 0, 0, 0] # no noise added
else:
stdIn = eval(stdIn)
#------------------------------------------------------------------------------
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
#------------------------------------------------------------------------------
outS = input("Give path and statistics file name: \n")
if ((outS == '') and (inF == "data_all.csv")):
outS = "RV_"+ext+"_stats.txt"
elif ((outS == '') and (inF == "era5_r360x180.nc")):
outS = "era5_"+ext+"_stats.txt"
elif (outS[-4:] != '.txt'):
outF = outS+".txt"
#------------------------------------------------------------------------------
print("\n run_ASFC.py, started for method "+meth)
res, lon, lat = toy_ASFC(inF, outF, outS, sst_fl, gustIn, cskinIn, tolIn,
meth, qmIn, LIn, stdIn)
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[ttl[i][:3]],'.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 statistics
if (cskinIn == None) and (meth in ["S80", "S88", "LP82", "YT96", "UA", "NCAR"]):
cskinIn = 0
elif (cskinIn == None) and (meth in ["C30", "C35", "ecmwf", "Beljaars"]):
cskinIn = 1
if np.all(gustIn == None) and (meth in ["C30", "C35"]):
gustIn = [1, 1.2, 600, 0.2]
elif np.all(gustIn == None) and (meth in ["UA", "ecmwf"]):
gustIn = [1, 1, 1000, 0.01]
elif np.all(gustIn == None):
gustIn = [1, 1.2, 800, 0.01]
elif (np.size(gustIn) < 4) and (gustIn == 0):
gust = [0, 0, 0, 0]
if tolIn == None:
tolIn = ['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1]
print("Input summary", file=open('./'+outS, 'a'))
print('input file name: {}, \n method: {}, \n gustiness: {}, \n cskin: {},'
' \n tolerance: {}, \n qmethod: {}, \n L: {}'.format(inF, meth, gustIn,
cskinIn, tolIn,
qmIn, LIn),
file=open('./'+outS, 'a'))
ttl = np.asarray(["tau ", "shf ", "lhf ",
"u10n ", "t10n ", "q10n "])
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('./'+outS, 'a'))
print('-'*79+'\n', file=open('./'+outS, 'a'))
elif (inF == "data_all.csv"):
a = res.loc[:, "tau":"q10n"].to_numpy(dtype="float64").T
stats = np.c_[stats, np.nanmean(a, axis=1)]
stats = np.c_[stats, np.nanmedian(a, axis=1)]
stats = np.c_[stats, np.nanmin(a, axis=1)]
stats = np.c_[stats, np.nanmax(a, axis=1)]
stats = np.c_[stats, np.nanpercentile(a, 5, axis=1)]
stats = np.c_[stats, np.nanpercentile(a, 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('./'+outS, 'a'))
print('-'*79+'\n', file=open('./'+outS, 'a'))
del 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'))