3. Roughness parameters¶
[1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path
from platypus.core import Problem
from platypus.types import Real,random
from platypus.algorithms import NSGAIII
import warnings
warnings.filterwarnings('ignore')
This is a custom package specifically designed for this analysis. It contains various functions for reading and computing and plotting.
[2]:
from z0_util import cal_vap_sat, cal_dens_dry, cal_dens_vap, cal_cpa, cal_dens_air, cal_Lob
3.1. Function to calculate Neutral condition¶
[13]:
def cal_neutral(df_val,z_meas,h_sfc):
# calculate Obukhov length
ser_Lob = df_val.apply(
lambda ser: cal_Lob(ser.H, ser.USTAR, ser.TA, ser.RH, ser.PA * 10), axis=1)
# zero-plane displacement: estimated using rule f thumb `d=0.7*h_sfc`
z_d = 0.7 * h_sfc
if z_d >= z_meas:
print(
'vegetation height is greater than measuring height. Please fix this before continuing'
)
# calculate stability scale
ser_zL = (z_meas - z_d) / ser_Lob
# determine periods under quasi-neutral conditions
limit_neutral = 0.01
ind_neutral = ser_zL.between(-limit_neutral, limit_neutral)
ind_neutral=ind_neutral[ind_neutral]
df_sel = df_val.loc[ind_neutral.index, ['WS', 'USTAR']].dropna()
ser_ustar = df_sel.USTAR
ser_ws = df_sel.WS
return ser_ws,ser_ustar
3.2. Function to calculate z0 and d using MO optimization¶
[14]:
def optimize_MO(df_val,z_meas,h_sfc):
ser_ws,ser_ustar=cal_neutral(df_val,z_meas,h_sfc)
def func_uz(params):
z0=params[0]
d=params[1]
z = z_meas
k = 0.4
uz = (ser_ustar / k) * np.log((z - d) / z0)
o1=abs(1-np.std(uz)/np.std(ser_ws))
o2=np.mean(abs(uz-ser_ws))/(np.mean(ser_ws))
return [o1,o2],[uz.min(),d-z0]
problem = Problem(2,2,2)
problem.types[0] = Real(0, 10)
problem.types[1] = Real(0, h_sfc)
problem.constraints[0] = ">=0"
problem.constraints[1] = ">=0"
problem.function = func_uz
random.seed(12345)
algorithm=NSGAIII(problem, divisions_outer=50)
algorithm.run(30000)
z0s=[]
ds=[]
os1=[]
os2=[]
for s in algorithm.result:
z0s.append(s.variables[0])
ds.append(s.variables[1])
os1.append(s.objectives[0])
os2.append(s.objectives[1])
idx=os2.index(min(os2, key=lambda x:abs(x-np.mean(os2))))
z0=z0s[idx]
d=ds[idx]
return z0,d,ser_ws,ser_ustar
3.3. Loading data, cleaning and getting ready for optimization¶
[15]:
name_of_site='US-MMS'
years=[2010,2012,2016]
df_attr=pd.read_csv('all_attrs.csv')
z_meas=df_attr[df_attr.site==name_of_site].meas.values[0]
h_sfc=df_zmeas=df_attr[df_attr.site==name_of_site].height.values[0]
folder='data/data_csv_zip_clean_roughness/'
site_file = folder+'/'+name_of_site + '_clean.csv.gz'
df_data = pd.read_csv(site_file, index_col='time', parse_dates=['time'])
# Rain
bb=pd.DataFrame(~np.isin(df_data.index.date,df_data[df_data.P!=0].index.date))
bb.index=df_data.index
df_data=df_data[bb.values]
df_data=df_data[(df_data['WS']!=0)]
df_years=df_data.loc[f'{years[0]}']
for i in years[1:]:
df_years=df_years.append(df_data.loc[f'{i}'])
df_val = df_years.loc[:, ['H', 'USTAR', 'TA', 'RH', 'PA', 'WS']].dropna()
df_val.head()
[15]:
H | USTAR | TA | RH | PA | WS | |
---|---|---|---|---|---|---|
time | ||||||
2010-01-01 00:00:00 | 21.873 | 0.739 | -8.08 | 70.503 | 98.836 | 3.695 |
2010-01-01 01:00:00 | 41.819 | 0.855 | -9.17 | 72.757 | 98.880 | 3.928 |
2010-01-01 02:00:00 | -6.078 | 0.699 | -9.63 | 72.611 | 98.910 | 3.088 |
2010-01-01 03:00:00 | -16.788 | 0.581 | -10.03 | 73.868 | 98.970 | 3.623 |
2010-01-01 04:00:00 | 5.006 | 0.562 | -10.36 | 74.242 | 99.030 | 3.474 |
3.4. Calculating z0 and d¶
[16]:
z0,d,ser_ws,ser_ustar=optimize_MO(df_val,z_meas,h_sfc)
3.5. Calculating model wind speed using logarithmic law¶
[17]:
def uz(z0,d,ser_ustar,z_meas):
z = z_meas
k = 0.4
uz = (ser_ustar / k) * np.log((z - d) / z0)
return uz
ws_model=uz(z0,d,ser_ustar,z_meas)
[18]:
plt.scatter(ser_ws,ws_model,color='k')
plt.xlabel('WS-obs')
plt.ylabel('WS-model')
plt.title(f'z0={np.round(z0,2)} & d={np.round(d,2)}')
plt.plot([0,10],[0,10],color='r',linestyle='--')
plt.ylim([0,10])
plt.xlim([0,10])
[18]:
(0.0, 10.0)