<jupyter_start><jupyter_text>dset = pd.DataFrame()train_featuresM = train_featuresMdset['attr'] = train_featuresM.columnsdset['importance'] = rfM.feature_importances_dset = dset.sort_values(by='importance', ascending=False)plt.figure(figsize=(20, 14))plt.barh(y=dset['attr'], width=dset['importance'], color='1976D2')plt.title('Random Forest Predictor Importances for Moisture Scale', fontsize=22, fontweight='bold', pad=20)plt.rc('xtick', labelsize=20) plt.rc('ytick', labelsize=20) plt.xlabel('Importance', fontsize=20, labelpad=10)plt.ylabel('Predictor', fontsize=20, labelpad=10)plt.savefig('RF-predictor-importance-MoistureScale-Idaho-30minflux.png')plt.show()dset = pd.DataFrame()train_featuresT = train_featuresTdset['attr'] = train_featuresT.columnsdset['importance'] = rfT.feature_importances_dset = dset.sort_values(by='importance', ascending=False)plt.figure(figsize=(20, 14))plt.barh(y=dset['attr'], width=dset['importance'], color='1976D2')plt.title('Random Forest Predictor Importances for Temperature Scale', fontsize=20, fontweight='bold', pad=20)plt.rc('xtick', labelsize=20) plt.rc('ytick', labelsize=20) plt.xlabel('Importance', fontsize=20, labelpad=10)plt.ylabel('Predictor', fontsize=20, labelpad=10)plt.savefig('RF-predictor-importance-TemperatureScale-Cabauw-30minflux.png')plt.show()dset = pd.DataFrame()train_featuresF = train_featuresFdset['attr'] = train_featuresF.columnsdset['importance'] = rfF.feature_importances_dset = dset.sort_values(by='importance', ascending=False)plt.figure(figsize=(20, 14))plt.barh(y=dset['attr'], width=dset['importance'], color='1976D2')plt.title('Random Forest Predictor Importances for Friction Velocity', fontsize=20, fontweight='bold', pad=20)plt.rc('xtick', labelsize=20) plt.rc('ytick', labelsize=20) plt.xlabel('Importance', fontsize=20, labelpad=10)plt.savefig('RF-predictor-importance-FrictionVelocity-Cabauw-30minflux.png')plt.ylabel('Predictor', fontsize=20, labelpad=10)plt.show()<jupyter_code>## Update this to train ANN rather than RF for different regimes
import tensorflow as tf
from tensorflow import keras
def build_model(train_data):
    model = keras.Sequential([
        keras.layers.Dense(64, activation=tf.nn.relu,
                          input_shape=(train_data.shape[1],)),
        keras.layers.Dense(64, activation=tf.nn.relu,
                          input_shape=(train_data.shape[1],)),
        keras.layers.Dense(64, activation=tf.nn.relu,
                          input_shape=(train_data.shape[1],)),
        keras.layers.Dense(1)
    ])
    optimizer = tf.keras.optimizers.Adam(lr=0.001)
    batch_size=32
    model.compile(loss='mae',
                 optimizer=optimizer,
                 metrics=['mae'])
    return model

# Train ANN on negative BRN predictions
EPOCHS = 25

######### Configure and train the Random Forest ##############
# Import the Random Forest Model
# Instantiate model 
ANN_M = build_model(np.asarray(train_featuresM))
ANN_T = build_model(np.asarray(train_featuresT))
ANN_F = build_model(np.asarray(train_featuresF))

# Train the model on training data
historyM = ANN_M.fit(np.asarray(train_featuresM), np.asarray(train_labelsM), epochs=EPOCHS, validation_split=0.2,verbose=0)
historyF = ANN_F.fit(np.asarray(train_featuresF), np.asarray(train_labelsF), epochs=EPOCHS, validation_split=0.2,verbose=0)
historyT = ANN_T.fit(np.asarray(train_featuresT), np.asarray(train_labelsT), epochs=EPOCHS, validation_split=0.2,verbose=0)

# Use the forest's predict method on the test data
ANNpredictionsM = ANN_M.predict(np.asarray(test_featuresM))
ANNpredictionsT = ANN_T.predict(np.asarray(test_featuresT))
ANNpredictionsF = ANN_F.predict(np.asarray(test_featuresF))

ANNpredictionsM = ANNpredictionsM.reshape(len(ANNpredictionsM,))
ANNpredictionsT = ANNpredictionsT.reshape(len(ANNpredictionsT,))
ANNpredictionsF = ANNpredictionsF.reshape(len(ANNpredictionsF,))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsM = abs(ANNpredictionsM - test_labelsM)
ANNerrorsT = abs(ANNpredictionsT - test_labelsT)
ANNerrorsF = abs(ANNpredictionsF - test_labelsF)
# Print out the mean absolute error (mae)
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsM), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsT), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsF), 4))
data = pd.read_csv("/Volumes/d1/mccandle/mlsurfacelayer/cabauw_derived_data_20191021_30min_raw.csv", index_col="Time", parse_dates=["Time"])
data = data.dropna()
split_date = pd.Timestamp("2011-01-01")
train_data = data.loc[:split_date]
test_data = data.loc[split_date:]


train_data["moisture availability_3 cm_"] = moisture_availability(train_data["soil water content_3 cm_m3 m-3"], field_capacity=0.47)
train_data['0m Temp K'] = train_data['temperature scale_surface_K']

test_data["moisture availability_3 cm_"] = moisture_availability(test_data["soil water content_3 cm_m3 m-3"], field_capacity=0.47)
test_data['0m Temp K'] = test_data['temperature scale_surface_K']


def skin_temperature(upward_longwave_irradiance_W_m2, emissivity=1.0):
    """
    Calculates the radiative skin temperature from the upward longwave irradiance.
    Args:
        upward_longwave_irradiance_W_m2: The upward longwave irradiance from the surface
        emissivity: How much energy is emitted in comparison to a blackbody. Ranges from 0 to 1.
    Returns:
        The skin temperature in Kelvin.
    """
    sigma = 5.673067e-8
    return (upward_longwave_irradiance_W_m2 / (emissivity * sigma)) ** 0.25
data.columns
input_cols = ['global horizontal irradiance_0 m_W m-2', 'zenith_0 m_degrees',
       'azimuth_0 m_degrees',
       'pressure_2 m_hPa', 'potential temperature_10 m_K',
       'potential temperature_2 m_K', 'mixing ratio_10 m_g kg-1',
       'relative humidity_10 m_%',
       'air density_10 m_kg m-3', 'wind speed_10 m_m s-1',
       'wind speed_20 m_m s-1',
       'wind speed_40 m_m s-1',
       'mixing ratio_2 m_g kg-1',
       'relative humidity_2 m_%', 'bulk richardson_10-2 m_']
output_col_1 = 'friction velocity_surface_m s-1'
# Test adding in Bulk Richardson Number as a Predictor
datasetC = pd.DataFrame()
# Create datasets for training and testing.  
# Predictors to use:
datasetC['10m Potential Temp K'] = train_data['potential temperature_10 m_K']
datasetC['LowLevel Wind Speed m/s'] = train_data['wind speed_10 m_m s-1']
#datasetC['LowLevel Wind Dir deg'] = train_data['wind direction_10 m_degrees']
datasetC['MidLevel Wind Speed m/s'] = train_data['wind speed_20 m_m s-1']
#datasetC['MidLevel Wind Dir deg'] = train_data['wind direction_20 m_degrees']
datasetC['TopLevel Wind Speed m/s'] = train_data['wind speed_40 m_m s-1']
#datasetC['TopLevel Wind Dir deg'] = train_data['wind direction_40 m_degrees']
datasetC['2m RH %'] = train_data['relative humidity_2 m_%']
datasetC['0m GHI w/m^2'] = train_data['global horizontal irradiance_0 m_W m-2']
datasetC['2m Pressure hPa'] = train_data['pressure_2 m_hPa']
datasetC['0m Solar Zenith Angle deg'] = train_data['zenith_0 m_degrees']
datasetC['TopLevel Soil Water Content m3/m3'] = train_data['soil water content_3 cm_m3 m-3']
datasetC['TopLevel Soil Temp K'] = train_data['soil temperature_4 cm_K']
datasetC['10m Bulk Richardson'] = train_data['bulk richardson_10-2 m_']
# These are potentially needed for other computations, but not as predictors
datasetC['10m Temp K'] = train_data['temperature_10 m_K']
datasetC['2m Temp K'] = train_data['temperature_2 m_K']
datasetC['10m U-Wind m/s'] = train_data['u wind_10 m_m s-1']
datasetC['10m V-Wind m/s'] = train_data['v wind_10 m_m s-1']
#datasetC['potential temperature skin change_10 m_K m-1'] = train_data['potential temperature skin change_10 m_K m-1']
#datasetC['potential temperature skin change_20 m_K m-1'] = train_data['potential temperature skin change_20 m_K m-1']
datasetC['moisture availability_3 cm_'] = train_data['moisture availability_3 cm_']
datasetC['2m Temp K'] = train_data['temperature_2 m_K']
datasetC['mixing_ratio'] = train_data['mixing ratio_2 m_g kg-1']
datasetC['0m Temp K']  = train_data['skin temperature_0 m_K']
#datasetC['mixing ratio skin change_10 m_g kg-1 m-1'] = train_data['mixing ratio skin change_10 m_g kg-1 m-1']
#datasetC['mixing ratio skin change_20 m_g kg-1 m-1'] = train_data['mixing ratio skin change_20 m_g kg-1 m-1']


datasetTestC = pd.DataFrame()
# Create datasets for training and testing.  
# Predictors to use:
datasetTestC['10m Potential Temp K'] = test_data['potential temperature_10 m_K']
datasetTestC['LowLevel Wind Speed m/s'] = test_data['wind speed_10 m_m s-1']
#datasetTestC['LowLevel Wind Dir deg'] = test_data['wind direction_10 m_degrees']
datasetTestC['MidLevel Wind Speed m/s'] = test_data['wind speed_20 m_m s-1']
#datasetTestC['MidLevel Wind Dir deg'] = test_data['wind direction_20 m_degrees']
datasetTestC['TopLevel Wind Speed m/s'] = test_data['wind speed_40 m_m s-1']
#datasetTestC['TopLevel Wind Dir deg'] = test_data['wind direction_40 m_degrees']
datasetTestC['2m RH %'] = test_data['relative humidity_2 m_%']
datasetTestC['0m GHI w/m^2'] = test_data['global horizontal irradiance_0 m_W m-2']
datasetTestC['2m Pressure hPa'] = test_data['pressure_2 m_hPa']
datasetTestC['0m Solar Zenith Angle deg'] = test_data['zenith_0 m_degrees']
datasetTestC['TopLevel Soil Water Content m3/m3'] = test_data['soil water content_3 cm_m3 m-3']
datasetTestC['TopLevel Soil Temp K'] = test_data['soil temperature_4 cm_K']
datasetTestC['10m Bulk Richardson'] = test_data['bulk richardson_10-2 m_']
# These are potentially needed for other computations, but not as predictors
datasetTestC['10m Temp K'] = test_data['temperature_10 m_K']
datasetTestC['2m Temp K'] = test_data['temperature_2 m_K']
datasetTestC['10m U-Wind m/s'] = test_data['u wind_10 m_m s-1']
datasetTestC['10m V-Wind m/s'] = test_data['v wind_10 m_m s-1']
#datasetTestC['potential temperature skin change_10 m_K m-1'] = test_data['potential temperature skin change_10 m_K m-1']
#datasetTestC['potential temperature skin change_20 m_K m-1'] = test_data['potential temperature skin change_20 m_K m-1']
datasetTestC['moisture availability_3 cm_'] = test_data['moisture availability_3 cm_']
datasetTestC['2m Temp K'] = test_data['temperature_2 m_K'] 
datasetTestC['mixing_ratio'] = test_data['mixing ratio_2 m_g kg-1']
datasetTestC['0m Temp K']  = test_data['skin temperature_0 m_K']
#datasetTestC['mixing ratio skin change_10 m_g kg-1 m-1'] = test_data['mixing ratio skin change_10 m_g kg-1 m-1']
#datasetTestC['mixing ratio skin change_20 m_g kg-1 m-1'] = test_data['mixing ratio skin change_20 m_g kg-1 m-1']



featuresFC = ['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']
featuresTC = ['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']
featuresMC = ['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']


MoistureScaleTrainC = train_data['moisture scale_surface_g kg-1']
TemperatureScaleTrainC = train_data['temperature scale_surface_K']
FrictionVelocityTrainC = train_data['friction velocity_surface_m s-1']

MoistureScaleTestC = test_data['moisture scale_surface_g kg-1']
TemperatureScaleTestC = test_data['temperature scale_surface_K']
FrictionVelocityTestC = test_data['friction velocity_surface_m s-1']

# Put in nomenclature so that either yearly or random split can be used
# Put in nomenclature so that either yearly or random split can be used
train_featuresMC = datasetC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
train_labelsMC = MoistureScaleTrainC
train_featuresFC = datasetC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
train_labelsFC = FrictionVelocityTrainC
train_featuresTC = datasetC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
train_labelsTC = TemperatureScaleTrainC

test_featuresMC = datasetTestC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
test_labelsMC = MoistureScaleTestC
test_featuresFC = datasetTestC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
test_labelsFC = FrictionVelocityTestC
test_featuresTC = datasetTestC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
test_labelsTC = TemperatureScaleTestC

######### Configure and train the Random Forest ##############
# Import the Random Forest Model
# Instantiate model 
rfMC = RandomForestRegressor(n_estimators= 200, random_state=2, max_leaf_nodes = 200, n_jobs=4,criterion='mse',min_samples_split=50,min_samples_leaf=50)
rfTC = RandomForestRegressor(n_estimators= 200, random_state=2, max_leaf_nodes = 200, n_jobs=4,criterion='mse',min_samples_split=50,min_samples_leaf=50)
rfFC = RandomForestRegressor(n_estimators= 200, random_state=2, max_leaf_nodes = 200, n_jobs=4,criterion='mse',min_samples_split=50,min_samples_leaf=50)
# Train the model on training data
rfMC.fit(train_featuresMC, train_labelsMC);
rfFC.fit(train_featuresFC, train_labelsFC);
rfTC.fit(train_featuresTC, train_labelsTC);

# Use the forest's predict method on the test data
predictionsMC = rfMC.predict(test_featuresMC)
predictionsTC = rfTC.predict(test_featuresTC)
predictionsFC = rfFC.predict(test_featuresFC)

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsMC = abs(predictionsMC - test_labelsMC)
errorsTC = abs(predictionsTC - test_labelsTC)
errorsFC = abs(predictionsFC - test_labelsFC)
# Print out the mean absolute error (mae)
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsMC), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsTC), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsFC), 4))
print('Moisture Scale Mean: ',np.mean(MoistureScaleTestC))
print('Moisture Scale Standard Deviation: ',np.std(MoistureScaleTestC))
print('Temperature Scale Mean: ',np.mean(TemperatureScaleTestC))
print('Temperature Scale Standard Deviation: ',np.std(TemperatureScaleTestC))
print('Friction Velocity Mean: ',np.mean(FrictionVelocityTestC))
print('Friction Velocity Standard Deviation: ',np.std(FrictionVelocityTestC))

# Feature selection
importancesMC = list(rfMC.feature_importances_)
feature_listMC = featuresMC
feature_importancesMC = [(featuresMC,round(importancesMC,2)) for featuresMC, importancesMC in zip(feature_listMC,importancesMC)]
feature_importancesMC = sorted(feature_importancesMC,key = lambda x: x[1],reverse=True)
print("Features importances for predicting moisture scale")
[print('Variable:{:20} Importance: {}'.format(*pair)) for pair in feature_importancesMC];

# Feature selection
importancesTC = list(rfTC.feature_importances_)
feature_listTC = featuresTC
feature_importancesTC = [(featuresTC,round(importancesTC,2)) for featuresTC, importancesTC in zip(feature_listTC,importancesTC)]
feature_importancesTC = sorted(feature_importancesTC,key = lambda x: x[1],reverse=True)
print("Features importances for predicting temperature scale")
[print('Variable:{:20} Importance: {}'.format(*pair)) for pair in feature_importancesTC];

# Feature selection
importancesFC = list(rfFC.feature_importances_)
feature_listFC = featuresFC
feature_importancesFC = [(featuresFC,round(importancesFC,2)) for featuresFC, importancesFC in zip(feature_listFC,importancesFC)]
feature_importancesFC = sorted(feature_importancesFC,key = lambda x: x[1],reverse=True)
print("Features importances for predicting friction velocity")
[print('Variable:{:20} Importance: {}'.format(*pair)) for pair in feature_importancesFC];
dset = pd.DataFrame()
tf_MC = train_featuresMC
dset['attr'] = tf_MC.columns
dset['importance'] = rfMC.feature_importances_

dset = dset.sort_values(by='importance', ascending=False)


plt.figure(figsize=(20, 14))
plt.barh(y=dset['attr'], width=dset['importance'], color='#1976D2')
plt.title('Random Forest Predictor Importances for Moisture Scale', fontsize=24, fontweight='bold', pad=20)
plt.rc('xtick', labelsize=24) 
plt.rc('ytick', labelsize=24) 
plt.xlabel('Importance', fontsize=24, labelpad=10)
#plt.ylabel('Predictor', fontsize=20, labelpad=10)
plt.savefig('RF-predictor-importance-MoistureScale-Cabauw-30minflux.png')
plt.show()

dset = pd.DataFrame()
tf_TC = train_featuresTC
dset['attr'] = tf_TC.columns
dset['importance'] = rfTC.feature_importances_

dset = dset.sort_values(by='importance', ascending=False)


plt.figure(figsize=(20, 14))
plt.barh(y=dset['attr'], width=dset['importance'], color='#1976D2')
plt.title('Random Forest Predictor Importances for Temperature Scale', fontsize=24, fontweight='bold', pad=20)
plt.rc('xtick', labelsize=24) 
plt.rc('ytick', labelsize=24) 
plt.xlabel('Importance', fontsize=24, labelpad=10)
#plt.ylabel('Predictor', fontsize=20, labelpad=10)
plt.savefig('RF-predictor-importance-TemperatureScale-Cabauw-30minflux.png')
plt.show()

dset = pd.DataFrame()
tf_FC = train_featuresFC
dset['attr'] = tf_FC.columns
dset['importance'] = rfFC.feature_importances_

dset = dset.sort_values(by='importance', ascending=False)


plt.figure(figsize=(20, 14))
plt.barh(y=dset['attr'], width=dset['importance'], color='#1976D2')
plt.title('Random Forest Predictor Importances for Friction Velocity', fontsize=24, fontweight='bold', pad=20)
plt.rc('xtick', labelsize=24) 
plt.rc('ytick', labelsize=24) 
plt.xlabel('Importance', fontsize=24, labelpad=10)
#plt.ylabel('Predictor', fontsize=20, labelpad=10)

plt.savefig('RF-predictor-importance-FrictionVelocity-Cabauw-30minflux.png')
plt.show()
## Update this to train ANN rather than RF for different regimes
import tensorflow as tf
from tensorflow import keras
def build_model(train_data):
    model = keras.Sequential([
        keras.layers.Dense(64, activation=tf.nn.relu,
                          input_shape=(train_data.shape[1],)),
        keras.layers.Dense(64, activation=tf.nn.relu,
                          input_shape=(train_data.shape[1],)),
        keras.layers.Dense(64, activation=tf.nn.relu,
                          input_shape=(train_data.shape[1],)),
        keras.layers.Dense(1)
    ])
    optimizer = tf.keras.optimizers.Adam(lr=0.001)
    batch_size=32
    model.compile(loss='mae',
                 optimizer=optimizer,
                 metrics=['mae'])
    return model

# Train ANN on negative BRN predictions
EPOCHS = 25

######### Configure and train the Random Forest ##############
# Import the Random Forest Model
# Instantiate model 
ANN_MC = build_model(np.asarray(train_featuresMC))
ANN_TC = build_model(np.asarray(train_featuresTC))
ANN_FC = build_model(np.asarray(train_featuresFC))

# Train the model on training data
historyMC = ANN_MC.fit(np.asarray(train_featuresMC), np.asarray(train_labelsMC), epochs=EPOCHS, validation_split=0.2,verbose=0)
historyFC = ANN_FC.fit(np.asarray(train_featuresFC), np.asarray(train_labelsFC), epochs=EPOCHS, validation_split=0.2,verbose=0)
historyTC = ANN_TC.fit(np.asarray(train_featuresTC), np.asarray(train_labelsTC), epochs=EPOCHS, validation_split=0.2,verbose=0)

# Use the forest's predict method on the test data
ANNpredictionsMC = ANN_MC.predict(np.asarray(test_featuresMC))
ANNpredictionsTC = ANN_TC.predict(np.asarray(test_featuresTC))
ANNpredictionsFC = ANN_FC.predict(np.asarray(test_featuresFC))

ANNpredictionsMC = ANNpredictionsMC.reshape(len(ANNpredictionsMC,))
ANNpredictionsTC = ANNpredictionsTC.reshape(len(ANNpredictionsTC,))
ANNpredictionsFC = ANNpredictionsFC.reshape(len(ANNpredictionsFC,))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsMC = abs(ANNpredictionsMC - test_labelsMC)
ANNerrorsTC = abs(ANNpredictionsTC - test_labelsTC)
ANNerrorsFC = abs(ANNpredictionsFC - test_labelsFC)
# Print out the mean absolute error (mae)
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsMC), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsTC), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsFC), 4))
from math import log, atan, sqrt
import numpy as np
from numba import jit
@jit(nopython=True)
#def mo_similarity(u10, v10, tsk, t2, qsfc, q2, psfc, z0=0.1, zt0=2, z10=10.0, z2=10.0):
#@jit(nopython=True)
def mo_similarity(u10, v10, tsk, t2, qsfc, q2, psfc, mavail=1, z0=0.01, zt0=0.001, z10=10.0, z2=2.0):
    """
    Calculate flux information based on Monin-Obukhov similarity theory.

    Args:
        u10: 10 m level zonal wind speed [m/s]
        v10: 10 m level meridional wind speed [m/s]
        tsk: surface skin temperature [K]
        t2: 2 m level temperature [K]
        qsfc: ground mixing ratio [kg /kg]
        q2: 2 m level mixing ratio [kg / kg]
        psfc: surface pressure [hPa]
        z0: momentum roughness length
        zt0: heat flux roughness length
        z10: Height of "10 m" layer
        z2: Height of "2 m" layer

    Returns:
        ustar: friction velocity
        tstar: temperature scale
        wthv0: -ustar * tstar
        zeta10: z/L for a z of 10 m. Used for momentum flux
        phim10: momentum universal function at 10 m
        zeta2: z/L for a z of 2 m.
        phih2: sensible heat flux universal function at 2 m.
    """
    if z0 > 0:
        z10oz0 = z10 / z0
        z2oz0 = z2 / z0
        z2ozt0 = z2 / zt0
    else:
        raise ValueError("z0 must be greater than 0")
    # else:
    #    sys.exit("Surface roughnes, z0, must be greter than 0.!")
    #
    # Gravitational acceleration
    g = 9.81
    #
    # Gas constant over spcific heat capacity at constant pressure
    r = 287.058
    cp = 1005.
    rocp = r / cp
    #
    # Reference pressure and temperature
    p0 = 1000.
    t0 = 300.
    #
    # Set M-O parameters based on Dyer 1974 paper
    karman = 0.4
    beta = 5.0
    gamma = 16.0
    #
    # Air density
    rho = psfc / (r * t2)
    #
    # Potential temperature speed at level 1
    th2 = t2 * (p0 / psfc) ** rocp
    th0 = tsk * (p0 / psfc) ** rocp
    #
    # Virtual potential temperature
    thv2 = th2 * (1. + 0.61 * q2)
    thv0 = th0 * (1. + 0.61 * qsfc)
    #
    # Small number
    epsilon = 1.e-6
    #
    # Initial values of drag coefficients - neutrally stratified case
    cd = karman ** 2 / ((log(z10oz0)) ** 2)
    ch = karman ** 2 / ((log(z2ozt0)) ** 2)
    cq = karman ** 2 / ((log(z2oz0)) ** 2)
    #
    # Initial values of surface friction velocity, temperature scale, and
    # heat flux
    wind_speed = sqrt(u10 * u10 + v10 * v10)
    if wind_speed < 0.1:
        wind_speed = 0.1
        u10 = 0.1
    tauxz = cd * wind_speed * u10
    tauyz = cd * wind_speed * v10
    ustar = (tauxz ** 2 + tauyz ** 2) ** 0.25

    tstar = -ch / ustar * wind_speed * (thv0 - thv2)
    qstar = cq / ustar * mavail * wind_speed * (qsfc - q2)
    wthv0 = -ustar * tstar
    #
    # Set stopping criterion
    diff = 1.
    #
    # Set stability functions
    psim10 = 0.
    psim2 = 0.
    psih2 = 0.
    psiq2 = 0.
    phim10 = 0.
    phih2 = 0
    #
    zeta10 = 0.
    zeta2 = 0.
    count = 0
    while diff > epsilon and count < 100:
        #
        # Surface friction velocity and temperature scale
        tauxz = cd * wind_speed * u10
        tauyz = cd * wind_speed * v10
        ustar = (tauxz ** 2 + tauyz ** 2) ** 0.25
        if ustar < 0.01:
            ustar = 0.01
        wspd2 = ustar / karman * (log(z2oz0) - psim2)
        tstar = -ch / ustar * wspd2 * (thv0 - thv2)
        qstar = cq / ustar * mavail * wspd2 * (qsfc - q2)
        wthv0 = -ustar * tstar
        #
        # Compute drag coefficients
        cdold = cd
        chold = ch
        cqold = cq
        #
        # Neutrally stratified case
        if wthv0 == 0:
            zeta10 = 0.
            zeta2 = 0.
            psim10 = 0.
            psim2 = 0.
            psih2 = 0.
            psiq2 = 0.
            phim10 = 1.
            phih2 = 1.
            cd = karman ** 2 / ((log(z10oz0)) ** 2)
            ch = karman ** 2 / ((log(z2ozt0)) ** 2)
            cq = karman ** 2 / ((log(z2oz0)) ** 2)
        elif abs(wthv0) > 0:
            #
            # Obukhov length scale
            olength = -ustar ** 3 / (karman * g / t0 * wthv0)
            if abs(olength) < 10 and olength > 0:
                olength = z10
            elif abs(olength) < 10 and olength < 0:
                olength = -z10
            #
            # Free convection
            # if (olength == 0.):
            #    sys.exit("Free convection!")
            #
            # Monin-Obukhov stability parameter
            zeta10 = z10 / olength
            zeta2 = z2 / olength
            #
            # Convective case
            if (zeta2 < -epsilon) & (zeta10 >= -2.):
                xi10 = 1. / ((1. - gamma * zeta10) ** 0.25)
                xi2 = 1. / ((1. - gamma * zeta2) ** 0.25)
                psim10 = log(0.5 * (1.0 + xi10 ** 2) * (0.5 * (1.0 + xi10)) ** 2) \
                    - 2. * atan(xi10) + 0.5 * np.pi
                psim2 = log(0.5 * (1.0 + xi2 ** 2) * (0.5 * (1.0 + xi2)) ** 2) \
                    - 2. * atan(xi2) + 0.5 * np.pi
                psih2 = 2.0 * log(0.5 * (1.0 + xi2 ** 2))
                psiq2 = 2.0 * log(0.5 * (1.0 + xi2 ** 2))
                phim10 = 1. / ((1. - gamma * zeta10) ** 0.25)
                phih2 = 1. / ((1. - gamma * zeta2) ** 0.25)
            #
            # Stably stratified case
            elif (zeta2 > epsilon) & (zeta2 <= 1.):
                psim10 = - beta * zeta10
                psim2 = - beta * zeta2
                psih2 = - beta * zeta2
                psiq2 = - beta * zeta2
                phim10 = (1. + beta * zeta10)
                phih2 = (1. + beta * zeta2)
            #
            # Neutrally stratified case
            elif (zeta2 <= epsilon) & (zeta2 >= -epsilon):
                psim10 = 0.
                psim2 = 0.
                psih2 = 0.
                psiq2 = 0.
                phim10 = 1.
                phih2 = 1.
            #
            cd = karman ** 2 / ((log(z10oz0) - psim10) ** 2)
            ch = karman ** 2 / ((log(z2ozt0) - psim2) * (log(z2ozt0) - psih2))
            cq = karman ** 2 / ((log(z2oz0) - psim2) * (log(z2oz0) - psiq2))
        #
        diff = abs(cd - cdold) + abs(ch - chold) + abs(cq - cqold)
        count += 1
    #
    return ustar, tstar, qstar, wthv0, zeta10, phim10, zeta2, phih2
# Test the MO Similarity on the test dataset for the RF Predictions #
#--- Inputs
# Wind velocity at 10 m
u10 = datasetTest['10m U-Wind m/s']
v10 = datasetTest['10m V-Wind m/s']

# Skin temperature and temperature at 2 m
tsk = datasetTest['0m Temp K']
t2  = datasetTest['10m Temp K'] 
#
# Ground mixing ratio and mixing ratio at 2 m
#qsfc, vap_pressure, sat_vap_pressure = mixing_ratio(datasetTest['0m Temp K']-273.15, 
#                    datasetTest['2m RH %'], datasetTest['2m Pressure hPa'])
q2, vap_pressure, sat_vap_pressure    = mixing_ratio(datasetTest['10m Temp K']-273.15, 
                    datasetTest['2m RH %'], datasetTest['2m Pressure hPa'])

def saturation_vapor_pressure(temperature_K):
    """
    Calculate saturation vapor pressure from Clausius-Clapeyron equation in A First Course in
    Atmospheric Thermodynamics by <NAME>
    Args:
        temperature_K: Temperature in Kelvin
    Returns:
        saturation vapor pressure in hPa/q2//q2
    """
    Rv = 461.0
    T0 = 273.0
    es0 = 6.11
    L = 2.5e6
    return es0 * np.exp(L / Rv * (1 / T0 - 1 / temperature_K))

qsfc = saturation_mixing_ratio(datasetTest['0m Temp K'], datasetTest['2m Pressure hPa'])


qsfc = qsfc / 1000
q2 = q2 /1000 #* datasetTest['2m RH %'] / 100
#
# Surface pressure 
psfc = datasetTest['2m Pressure hPa']
#
# Momentum roughness length, heat flux roughness and reference level
z0=0.01 # might want to go back to 0.1 
zt0=0.001 
z10=10.0 
z2=2.0

# Outside the range of M-O Similarity Theory
phim=[]
zolm=[]
phih=[]
zolh=[]
t_2_out=[]
t_sk_out=[]
ustar_out=[]
tstar_out=[]
pred_ustar_out = []
pred_tstar_out = []
actual_ustar_out = []
actual_tstar_out = []
qstar_out = []

# Within the range of M-O Similarity Theory
phimO=[]
zolmO=[]
phihO=[]
zolhO=[]
t_2_in=[]
t_sk_in=[]
ustar_in=[]
tstar_in=[]
pred_ustar_in = []
pred_tstar_in = []
pred_qstar_in = []
actual_ustar_in = []
actual_tstar_in = []
actual_qstar_in = []
qstar_in = []

#ustar,tstar,wthv0,zeta10,phim10,zeta2,phih2 = mo_similarity(u10,v10,tsk,t2,qsfc,q2,psfc,z0,zt0)
#
for i, val in enumerate(u10):
    ustar,tstar,qstar,wthv0,zeta10,phim10,zeta2,phih2 = mo_similarity(u10.iloc[i],v10.iloc[i],tsk.iloc[i],t2.iloc[i],qsfc.iloc[i],q2.iloc[i],psfc.iloc[i],mavail=1, z0=0.01, zt0=0.001, z10=10.0, z2=2.0)
    if ((zeta10 < -2.) | (zeta10 > 1.)):
        phimO.append(phim10)
        zolmO.append(zeta10)
        phihO.append(phih2)
        zolhO.append(zeta2)
        t_2_out.append(t2.iloc[i])
        t_sk_out.append(tsk.iloc[i])
        ustar_out.append(ustar)
        tstar_out.append(tstar)
        qstar_out.append(qstar)
        pred_ustar_out.append(predictionsF[i])
        pred_tstar_out.append(predictionsT[i])
        pred_qstar_out.append(predictionsM[i])
        actual_ustar_out.append(test_labelsF.iloc[i])
        actual_tstar_out.append(test_labelsT.iloc[i])
        actual_qstar_out.append(test_labelsM.iloc[i])
        
    else:
        phim.append(phim10)
        zolm.append(zeta10)
        phih.append(phih2)
        zolh.append(zeta2)
        t_2_in.append(t2.iloc[i])
        t_sk_in.append(tsk.iloc[i])
        ustar_in.append(ustar)
        tstar_in.append(tstar)
        qstar_in.append(qstar)
        pred_ustar_in.append(predictionsF[i])
        pred_tstar_in.append(predictionsT[i])
        pred_qstar_in.append(predictionsM[i])
        actual_ustar_in.append(test_labelsF.iloc[i])
        actual_tstar_in.append(test_labelsT.iloc[i])
        actual_qstar_in.append(test_labelsM.iloc[i])
#
p=plt.plot(zolm,phim,'ro')
plt.title('Plot of Phi M vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi M', size = 16)
plt.show()
q=plt.plot(zolh,phih,'ro')
plt.title('Plot of Phi H vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi H', size = 16)
plt.show()
p=plt.plot(zolmO,phimO,'ro')
plt.title('Plot of Phi M vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi M', size = 16)
plt.show()
q=plt.plot(zolhO,phihO,'ro')
plt.title('Plot of Phi H vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi H', size = 16)
plt.show()

datasetTest['ustar_in'] = ustar_in
datasetTest['tstar_in'] = tstar_in
datasetTest['qstar_in'] = qstar_in
plt.hist(qsfc)
plt.hist(q2)
plt.hist(sat_vap_pressure)
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsT
y = test_labelsT
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
#plt.ylim(-4.5, 0.9)
#plt.xlim(-1.2,1)
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
#plt.text(-1,-3,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF Predicted vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_Observed_tscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsT
y = test_labelsT
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
#plt.ylim(-4.5, 0.9)
#plt.xlim(-1.2,1)
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
#plt.text(-1,-3,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN Predicted vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_Observed_tscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
y = np.array(actual_tstar_in)
x = np.array(tstar_in) 
#x = x[y>-5] # Filter out a bad data point
#y = x[x>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
plt.ylim(-4.5, 0.9)
plt.xlim(-1.2,1)
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
#plt.text(-1,-3,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('M-O Similarity vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('M-O Similarity', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('MO_Observed_tscale_idaho_30min.png')
plt.show()

print("MAE for MO Similarity on T-Scale...")
errorsMO_TI_C = abs(x - y)
print(np.mean(errorsMO_TI_C))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsM
y = test_labelsM
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF Predicted vs Observed Moisture Scale for Idaho ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Observed_mscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsM
y = test_labelsM
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN Predicted vs Observed Moisture Scale for Idaho ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Observed_mscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]

x = np.array(qstar_in)
y = np.array(actual_qstar_in) 

y = y[x>-5] # Filter out a bad data point
x = x[x>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('M-O Similarity vs Observed Moisture Scale for Idaho', size = 24)
plt.xlabel('M-O Similarity', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('MO_Observed_mscale_idaho_30min.png')
plt.show()

print("MAE for MO Similarity on Moisture Scale...")
errorsMO_FVI_C = abs(x - y)
print(np.mean(errorsMO_FVI_C))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsF
y = test_labelsF
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF Predicted vs Observed Friction Velocity for Idaho', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Observed_friction-velocity_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsF
y = test_labelsF
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN Predicted vs Observed Friction Velocity for Idaho', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Observed_friction-velocity_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = np.array(actual_ustar_in)
y = np.array(ustar_in) 
y = y[x>-5] # Filter out a bad data point
x = x[x>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('M-O Similarity vs Observed Friction Velocity for Idaho', size = 24)
plt.xlabel('M-O Similarity', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('MO_Observed_friction-velocity_idaho_30min.png')
plt.show()

print("MAE for MO Similarity on Friction Velocity...")
errorsMO_FVI_C = abs(x - y)
print(np.mean(errorsMO_FVI_C))
# Test the MO Similarity on the test dataset for the RF Predictions #
#--- Inputs
# Wind velocity at 10 m


u10 = datasetTestC['10m U-Wind m/s']
v10 = datasetTestC['10m V-Wind m/s']

# Skin temperature and temperature at 2 m
tsk = datasetTestC['0m Temp K']
t2  = datasetTestC['10m Temp K'] 
#
# Ground mixing ratio and mixing ratio at 2 m
#qsfc, vap_pressure, sat_vap_pressure = mixing_ratio(datasetTestC['0m Temp K']-273.15, 
#                    datasetTestC['2m RH %'], datasetTestC['2m Pressure hPa'])
q2, vap_pressure, sat_vap_pressure = mixing_ratio(datasetTestC['10m Temp K']-273.15, 
                    datasetTestC['2m RH %'], datasetTestC['2m Pressure hPa'])

def saturation_vapor_pressure(temperature_K):
    """
    Calculate saturation vapor pressure from Clausius-Clapeyron equation in A First Course in
    Atmospheric Thermodynamics by <NAME>
    Args:
        temperature_K: Temperature in Kelvin
    Returns:
        saturation vapor pressure in hPa/q2//q2
    """
    Rv = 461.0
    T0 = 273.0
    es0 = 6.11
    L = 2.5e6
    return es0 * np.exp(L / Rv * (1 / T0 - 1 / temperature_K))

qsfc = saturation_mixing_ratio(datasetTestC['0m Temp K'], datasetTestC['2m Pressure hPa'])

qsfc = qsfc / 1000
q2 = q2 /1000
#
# Surface pressure 
psfc = datasetTestC['2m Pressure hPa']
#
# Momentum roughness length, heat flux roughness and reference level
z0=0.01 # might want to go back to 0.1 
zt0=0.001 
z10=10.0 
z2=2.0

# Outside the range of M-O Similarity Theory
phim=[]
zolm=[]
phih=[]
zolh=[]
t_2_out=[]
t_sk_out=[]
ustar_out=[]
tstar_out=[]
pred_ustar_out = []
pred_tstar_out = []
actual_ustar_out = []
actual_tstar_out = []
qstar_out = []

# Within the range of M-O Similarity Theory
phimO=[]
zolmO=[]
phihO=[]
zolhO=[]
t_2_in=[]
t_sk_in=[]
ustar_in=[]
tstar_in=[]
pred_ustar_in = []
pred_tstar_in = []
pred_qstar_in = []
actual_ustar_in = []
actual_tstar_in = []
actual_qstar_in = []
qstar_in = []

#ustar,tstar,wthv0,zeta10,phim10,zeta2,phih2 = mo_similarity(u10,v10,tsk,t2,qsfc,q2,psfc,z0,zt0)
#
for i, val in enumerate(u10):
    ustar,tstar,qstar,wthv0,zeta10,phim10,zeta2,phih2 = mo_similarity(u10.iloc[i],v10.iloc[i],tsk.iloc[i],t2.iloc[i],qsfc.iloc[i],q2.iloc[i],psfc.iloc[i],mavail=1, z0=0.01, zt0=0.001, z10=10.0, z2=2.0)
    if ((zeta10 < -2.) | (zeta10 > 1.)):
        phimO.append(phim10)
        zolmO.append(zeta10)
        phihO.append(phih2)
        zolhO.append(zeta2)
        t_2_out.append(t2.iloc[i])
        t_sk_out.append(tsk.iloc[i])
        ustar_out.append(ustar)
        tstar_out.append(tstar)
        qstar_out.append(qstar)
        pred_ustar_out.append(predictionsFC[i])
        pred_tstar_out.append(predictionsTC[i])
        pred_qstar_out.append(predictionsMC[i])
        actual_ustar_out.append(test_labelsFC.iloc[i])
        actual_tstar_out.append(test_labelsTC.iloc[i])
        actual_qstar_out.append(test_labelsMC.iloc[i])
        
    else:
        phim.append(phim10)
        zolm.append(zeta10)
        phih.append(phih2)
        zolh.append(zeta2)
        t_2_in.append(t2.iloc[i])
        t_sk_in.append(tsk.iloc[i])
        ustar_in.append(ustar)
        tstar_in.append(tstar)
        qstar_in.append(qstar)
        pred_ustar_in.append(predictionsFC[i])
        pred_tstar_in.append(predictionsTC[i])
        pred_qstar_in.append(predictionsMC[i])
        actual_ustar_in.append(test_labelsFC.iloc[i])
        actual_tstar_in.append(test_labelsTC.iloc[i])
        actual_qstar_in.append(test_labelsMC.iloc[i])
#
p=plt.plot(zolm,phim,'ro')
plt.title('Plot of Phi M vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi M', size = 16)
plt.show()
q=plt.plot(zolh,phih,'ro')
plt.title('Plot of Phi H vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi H', size = 16)
plt.show()
p=plt.plot(zolmO,phimO,'ro')
plt.title('Plot of Phi M vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi M', size = 16)
plt.show()
q=plt.plot(zolhO,phihO,'ro')
plt.title('Plot of Phi H vs M-O Stability Parameter for Idaho', size = 16)
plt.xlabel('M-O Stability Parameter', size = 16)
plt.ylabel('Phi H', size = 16)
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsTC
y = test_labelsTC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF Predicted vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_Observed_tscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsTC
y = test_labelsTC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN Predicted vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_Observed_tscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = np.array(actual_tstar_in)
y = np.array(tstar_in) 
y = y[x>-5] # Filter out a bad data point
x = x[x>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('M-O Similarity vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('M-O Similarity', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('MO_Observed_tscale_cabauw_30min.png')
plt.show()

print("MAE for MO Similarity on T-Scale...")
errorsMO_T_C = abs(x - y)
print(np.mean(errorsMO_T_C))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsMC
y = test_labelsMC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF Predicted vs Observed Moisture Scale for Cabauw ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Observed_mscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsMC
y = test_labelsMC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN Predicted vs Observed Moisture Scale for Cabauw ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Observed_mscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]

x = np.array(qstar_in)
y = np.array(actual_qstar_in) 

y = y[x>-5] # Filter out a bad data point
x = x[x>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('M-O Similarity vs Observed Moisture Scale for Cabauw', size = 24)
plt.xlabel('M-O Similarity', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('MO_Observed_mscale_cabauw_30min.png')
plt.show()

print("MAE for MO Similarity on Moisture Scale...")
errorsMO_FVI_C = abs(x - y)
print(np.mean(errorsMO_FVI_C))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsFC
y = test_labelsFC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF Predicted vs Observed Friction Velocity for Cabauw', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Observed_friction-velocity_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsFC
y = test_labelsFC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN Predicted vs Observed Friction Velocity for Cabauw', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Observed_friction-velocity_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = np.array(actual_ustar_in)
y = np.array(ustar_in) 
y = y[x>-5] # Filter out a bad data point
x = x[x>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('M-O Similarity vs Observed Friction Velocity for Cabauw', size = 24)
plt.xlabel('M-O Similarity', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('MO_Observed_friction-velocity_cabauw_30min.png')
plt.show()

print("MAE for MO Similarity on Friction Velocity...")
errorsMO_FV_C = abs(x - y)
print(np.mean(errorsMO_FV_C))

plt.rcParams['figure.figsize'] = [10, 10]
for i in range(1,3):
    
    # Set up the plot
    ax = plt.subplot(2, 1, i)
    
    # Draw the plot
    if i == 1:
        ax.hist(datasetTest['10m U-Wind m/s'],bins=50)
        ax.set_title('10m U-Wind at Idaho', size = 10)
    elif i == 2:
        ax.hist(datasetTestC['10m U-Wind m/s'],bins=50)
        ax.set_title('10m U-Wind at Cabauw', size = 10)
    
    # Title and labels
#    ax.set_title('Histogram with Binwidth = %d' % binwidth, size = 30)
#    ax.set_xlabel('Delay (min)', size = 22)
#    ax.set_ylabel('Flights', size= 22)

plt.tight_layout()
plt.show()
#MoistureScaleTestC = test_data['moisture scale_surface_g kg-1']
#TemperatureScaleTestC = test_data['temperature scale_surface_K']
#FrictionVelocityTestC = test_data['friction velocity_surface_m s-1']

plt.rcParams['figure.figsize'] = [10, 10]
for i in range(1,3):
    
    # Set up the plot
    ax = plt.subplot(2, 1, i)
    
    # Draw the plot
    if i == 1:
        ax.hist(MoistureScaleTest,range =(-2, 2),bins=50)
        ax.set_title('Moisture Scale at Idaho', size = 16)
    elif i == 2:
        ax.hist(MoistureScaleTestC,range =(-2, 2),bins=50)
        ax.set_title('Moisture Scale at Cabauw', size = 16)
    
    # Title and labels
#    ax.set_title('Histogram with Binwidth = %d' % binwidth, size = 30)
#    ax.set_xlabel('Delay (min)', size = 22)
#    ax.set_ylabel('Flights', size= 22)

plt.tight_layout()
plt.show()
#MoistureScaleTestC = test_data['moisture scale_surface_g kg-1']
#TemperatureScaleTestC = test_data['temperature scale_surface_K']
#FrictionVelocityTestC = test_data['friction velocity_surface_m s-1']

plt.rcParams['figure.figsize'] = [10, 10]
for i in range(1,3):
    
    # Set up the plot
    ax = plt.subplot(2, 1, i)
    
    # Draw the plot
    if i == 1:
        ax.hist(TemperatureScaleTest,range =(-2.5, 1),bins=50)
        ax.set_title('Temperature Scale at Idaho', size = 16)
    elif i == 2:
        ax.hist(TemperatureScaleTestC,range =(-2.5, 1),bins=50)
        ax.set_title('Temperature Scale at Cabauw', size = 16)
    
    # Title and labels
#    ax.set_title('Histogram with Binwidth = %d' % binwidth, size = 30)
#    ax.set_xlabel('Delay (min)', size = 22)
#    ax.set_ylabel('Flights', size= 22)

plt.tight_layout()
plt.show()
#MoistureScaleTestC = test_data['moisture scale_surface_g kg-1']
#TemperatureScaleTestC = test_data['temperature scale_surface_K']
#FrictionVelocityTestC = test_data['friction velocity_surface_m s-1']

plt.rcParams['figure.figsize'] = [10, 10]
for i in range(1,3):
    
    # Set up the plot
    ax = plt.subplot(2, 1, i)
    
    # Draw the plot
    if i == 1:
        ax.hist(FrictionVelocityTest,range =(0, 1.2),bins=50)
        ax.set_title('Friction Velocity at Idaho', size = 16)
    elif i == 2:
        ax.hist(FrictionVelocityTestC,range =(0, 1.2),bins=50)
        ax.set_title('Friction Velocity at Cabauw', size = 16)
    
    # Title and labels
#    ax.set_title('Histogram with Binwidth = %d' % binwidth, size = 30)
#    ax.set_xlabel('Delay (min)', size = 22)
#    ax.set_ylabel('Flights', size= 22)

plt.tight_layout()
plt.show()
# Use the forest's predict method on the test data
predictionsMI_C = rfM.predict(test_featuresMC)
predictionsTI_C = rfT.predict(test_featuresTC)
predictionsFI_C = rfF.predict(test_featuresFC)

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsMI_C = abs(predictionsMI_C - test_labelsMC)
errorsTI_C = abs(predictionsTI_C - test_labelsTC)
errorsFI_C = abs(predictionsFI_C - test_labelsFC)
# Print out the mean absolute error (mae)
print('Errors for the random forest trained on Idaho and applied to Cabauw')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsMI_C), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsTI_C), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsFI_C), 4))

# Use the ANN's predict method on the test data
ANNpredictionsMI_C = ANN_M.predict(np.asarray(test_featuresMC))
ANNpredictionsTI_C = ANN_T.predict(np.asarray(test_featuresTC))
ANNpredictionsFI_C = ANN_F.predict(np.asarray(test_featuresFC))

ANNpredictionsMI_C = ANNpredictionsMI_C.reshape(len(ANNpredictionsMI_C,))
ANNpredictionsTI_C = ANNpredictionsTI_C.reshape(len(ANNpredictionsTI_C,))
ANNpredictionsFI_C = ANNpredictionsFI_C.reshape(len(ANNpredictionsFI_C,))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsMI_C = abs(ANNpredictionsMI_C - test_labelsMC)
ANNerrorsTI_C = abs(ANNpredictionsTI_C - test_labelsTC)
ANNerrorsFI_C = abs(ANNpredictionsFI_C - test_labelsFC)
# Print out the mean absolute error (mae)
print('Errors for the ANN trained on Idaho and applied to Cabauw')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsMI_C), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsTI_C), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsFI_C), 4))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsTI_C
y = test_labelsTC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho Trained - Cabauw RF vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_Idaho-on-Cabauw_Observed_tscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsTI_C
y = test_labelsTC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho Trained - Cabauw ANN vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_Idaho-on-Cabauw_Observed_tscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsMI_C
y = test_labelsMC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho Trained - Cabauw RF vs Observed Moisture Scale for Cabauw ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Idaho-on-Cabauw_Observed_mscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsMI_C
y = test_labelsMC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho Trained - Cabauw ANN vs Observed Moisture Scale for Cabauw ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
plt.savefig('ANN_Idaho-on-Cabauw_Observed_mscale_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsFI_C
y = test_labelsFC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho Trained - Cabauw RF vs Observed Friction Velocity for Cabauw', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Idaho-on-Cabauw_Observed_friction-velocity_cabauw_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsFI_C
y = test_labelsFC
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho Trained - Cabauw ANN vs Observed Friction Velocity for Cabauw', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Idaho-on-Cabauw_Observed_friction-velocity_cabauw_30min.png')
plt.show()
# Use the forest's predict method on the test data
predictionsMC_I = rfMC.predict(test_featuresM)
predictionsTC_I = rfTC.predict(test_featuresT)
predictionsFC_I = rfFC.predict(test_featuresF)

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsMC_I = abs(predictionsMC_I - test_labelsM)
errorsTC_I = abs(predictionsTC_I - test_labelsT)
errorsFC_I = abs(predictionsFC_I - test_labelsF)
# Print out the mean absolute error (mae)
print('Cabauw Trained and Idaho Predicted Random Forest')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsMC_I), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsTC_I), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsFC_I), 4))

# Use the ANN's predict method on the test data
ANNpredictionsMC_I = ANN_MC.predict(np.asarray(test_featuresM))
ANNpredictionsTC_I = ANN_TC.predict(np.asarray(test_featuresT))
ANNpredictionsFC_I = ANN_FC.predict(np.asarray(test_featuresF))

ANNpredictionsMC_I = ANNpredictionsMC_I.reshape(len(ANNpredictionsMC_I,))
ANNpredictionsTC_I = ANNpredictionsTC_I.reshape(len(ANNpredictionsTC_I,))
ANNpredictionsFC_I = ANNpredictionsFC_I.reshape(len(ANNpredictionsFC_I,))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsMC_I = abs(ANNpredictionsMC_I - test_labelsM)
ANNerrorsTC_I = abs(ANNpredictionsTC_I - test_labelsT)
ANNerrorsFC_I = abs(ANNpredictionsFC_I - test_labelsF)
# Print out the mean absolute error (mae)
print('Cabauw Trained and Idaho Predicted ANN')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsMC_I), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsTC_I), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsFC_I), 4))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsTC_I
y = test_labelsT
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Cabauw Trained - Idaho RF vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_Cabauw-on-Idaho_Observed_tscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsTC_I
y = test_labelsT
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Cabauw Trained - Idaho ANN vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_Cabauw-on-Idaho_Observed_tscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsMC_I
y = test_labelsM
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Cabauw Trained - Idaho RF vs Observed Moisture Scale for Idaho ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Cabauw-on-Idaho_Observed_mscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsMC_I
y = test_labelsM
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Cabauw Trained - Idaho ANN vs Observed Moisture Scale for Idaho ($q*$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Cabauw-on-Idaho_Observed_mscale_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsFC_I
y = test_labelsF
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Cabauw Trained - Idaho RF vs Observed Friction Velocity for Idaho', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('RF_Cabauw-on-Idaho_Observed_friction-velocity_idaho_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsFC_I
y = test_labelsF
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Cabauw Trained - Idaho ANN vs Observed Friction Velocity for Idaho', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(0, 1.5), np.linspace(0, 1.5), 'k--')
plt.savefig('ANN_Cabauw-on-Idaho_Observed_friction-velocity_idaho_30min.png')
plt.show()
## Add training datasets together and predicted on combined test datasets
trainingM = [train_featuresM,train_featuresMC]
trainM = pd.concat(trainingM)
trainingF = [train_featuresF,train_featuresFC]
trainF = pd.concat(trainingF)
trainingT = [train_featuresT,train_featuresTC]
trainT = pd.concat(trainingT)

testingM = [test_featuresM,test_featuresMC]
testM = pd.concat(testingM)
testingF = [test_featuresF,test_featuresFC]
testF = pd.concat(testingF)
testingT = [test_featuresT,test_featuresTC]
testT = pd.concat(testingT)

training_labelsM = [train_labelsM,train_labelsMC]
train_labelsMA = pd.concat(training_labelsM)
training_labelsF = [train_labelsF,train_labelsFC]
train_labelsFA = pd.concat(training_labelsF)
training_labelsT = [train_labelsT,train_labelsTC]
train_labelsTA = pd.concat(training_labelsT)

testing_labelsM = [test_labelsM,test_labelsMC]
test_labelsMA = pd.concat(testing_labelsM)
testing_labelsF = [test_labelsF,test_labelsFC]
test_labelsFA = pd.concat(testing_labelsF)
testing_labelsT = [test_labelsT,test_labelsTC]
test_labelsTA = pd.concat(testing_labelsT)
######### Configure and train the Random Forest ##############
# Import the Random Forest Model
# Instantiate model 
rfM_All = RandomForestRegressor(n_estimators= 200, random_state=2, max_leaf_nodes = 200, n_jobs=4,criterion='mse',min_samples_split=50,min_samples_leaf=50)
rfT_All = RandomForestRegressor(n_estimators= 200, random_state=2, max_leaf_nodes = 200, n_jobs=4,criterion='mse',min_samples_split=50,min_samples_leaf=50)
rfF_All = RandomForestRegressor(n_estimators= 200, random_state=2, max_leaf_nodes = 200, n_jobs=4,criterion='mse',min_samples_split=50,min_samples_leaf=50)
# Train the model on training data
rfM_All.fit(trainM, train_labelsMA);
rfF_All.fit(trainF, train_labelsFA);
rfT_All.fit(trainT, train_labelsTA);

# Use the forest's predict method on the test data
predictionsM_All = rfM_All.predict(testM)
predictionsT_All = rfT_All.predict(testT)
predictionsF_All = rfF_All.predict(testF)

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsM_All = abs(predictionsM_All - test_labelsMA)
errorsT_All = abs(predictionsT_All - test_labelsTA)
errorsF_All = abs(predictionsF_All - test_labelsFA)
# Print out the mean absolute error (mae)
print('For the RF trained on all...')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsM_All), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsT_All), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsF_All), 4))

######### Configure and train the Random Forest ##############
# Import the Random Forest Model
# Instantiate model 
ANN_MC_All = build_model(np.asarray(trainM))
ANN_TC_All = build_model(np.asarray(trainT))
ANN_FC_All = build_model(np.asarray(trainF))

# Train the model on training data
historyMC_All = ANN_MC_All.fit(np.asarray(trainM), np.asarray(train_labelsMA), epochs=EPOCHS, validation_split=0.2,verbose=0)
historyFC_All = ANN_FC_All.fit(np.asarray(trainF), np.asarray(train_labelsFA), epochs=EPOCHS, validation_split=0.2,verbose=0)
historyTC_All = ANN_TC_All.fit(np.asarray(trainT), np.asarray(train_labelsTA), epochs=EPOCHS, validation_split=0.2,verbose=0)

# Use the forest's predict method on the test data
ANNpredictionsMC_All = ANN_MC_All.predict(np.asarray(testM))
ANNpredictionsTC_All = ANN_TC_All.predict(np.asarray(testT))
ANNpredictionsFC_All = ANN_FC_All.predict(np.asarray(testF))

ANNpredictionsMC_All = ANNpredictionsMC_All.reshape(len(ANNpredictionsMC_All,))
ANNpredictionsTC_All = ANNpredictionsTC_All.reshape(len(ANNpredictionsTC_All,))
ANNpredictionsFC_All = ANNpredictionsFC_All.reshape(len(ANNpredictionsFC_All,))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsMC_All = abs(ANNpredictionsMC_All - test_labelsMA)
ANNerrorsTC_All = abs(ANNpredictionsTC_All - test_labelsTA)
ANNerrorsFC_All = abs(ANNpredictionsFC_All - test_labelsFA)
# Print out the mean absolute error (mae)
print('For the ANN trained on all...')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsMC_All), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsTC_All), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsFC_All), 4))
# Use the RF-All to predict the Idaho and Cabauw Test Datasets
predictionsM_Idaho_All = rfM_All.predict(test_featuresM)
predictionsT_Idaho_All = rfT_All.predict(test_featuresT)
predictionsF_Idaho_All = rfF_All.predict(test_featuresF)

predictionsM_Cabauw_All = rfM_All.predict(test_featuresMC)
predictionsT_Cabauw_All = rfT_All.predict(test_featuresTC)
predictionsF_Cabauw_All = rfF_All.predict(test_featuresFC)

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsM_All_C = abs(predictionsM_Cabauw_All - test_labelsMC)
errorsT_All_C = abs(predictionsT_Cabauw_All - test_labelsTC)
errorsF_All_C = abs(predictionsF_Cabauw_All - test_labelsFC)
# Print out the mean absolute error (mae)
print('RF trained on all and applied to Cabauw')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsM_All_C), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsT_All_C), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsF_All_C), 4))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsM_All_I = abs(predictionsM_Idaho_All - test_labelsM)
errorsT_All_I = abs(predictionsT_Idaho_All - test_labelsT)
errorsF_All_I = abs(predictionsF_Idaho_All - test_labelsF)
# Print out the mean absolute error (mae)
print('RF trained on all and applied to Idaho')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsM_All_I), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsT_All_I), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsF_All_I), 4))


######## Move to using the ANN.... ############
# Use the ANN-All to predict the Idaho and Cabauw Test Datasets
ANNpredictionsM_Idaho_All = ANN_MC_All.predict(np.asarray(test_featuresM))
ANNpredictionsT_Idaho_All = ANN_TC_All.predict(np.asarray(test_featuresT))
ANNpredictionsF_Idaho_All = ANN_FC_All.predict(np.asarray(test_featuresF))

ANNpredictionsM_Idaho_All = ANNpredictionsM_Idaho_All.reshape(len(ANNpredictionsM_Idaho_All,))
ANNpredictionsT_Idaho_All = ANNpredictionsT_Idaho_All.reshape(len(ANNpredictionsT_Idaho_All,))
ANNpredictionsF_Idaho_All = ANNpredictionsF_Idaho_All.reshape(len(ANNpredictionsF_Idaho_All,))

ANNpredictionsM_Cabauw_All = ANN_MC_All.predict(np.asarray(test_featuresMC))
ANNpredictionsT_Cabauw_All = ANN_TC_All.predict(np.asarray(test_featuresTC))
ANNpredictionsF_Cabauw_All = ANN_FC_All.predict(np.asarray(test_featuresFC))

ANNpredictionsM_Cabauw_All = ANNpredictionsM_Cabauw_All.reshape(len(ANNpredictionsM_Cabauw_All,))
ANNpredictionsT_Cabauw_All = ANNpredictionsT_Cabauw_All.reshape(len(ANNpredictionsT_Cabauw_All,))
ANNpredictionsF_Cabauw_All = ANNpredictionsF_Cabauw_All.reshape(len(ANNpredictionsF_Cabauw_All,))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsM_All_C = abs(ANNpredictionsM_Cabauw_All - test_labelsMC)
ANNerrorsT_All_C = abs(ANNpredictionsT_Cabauw_All - test_labelsTC)
ANNerrorsF_All_C = abs(ANNpredictionsF_Cabauw_All - test_labelsFC)
# Print out the mean absolute error (mae)
print('ANN trained on all and applied to Cabauw')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsM_All_C), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsT_All_C), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsF_All_C), 4))
     
     ######### Compute Errors on Test Data #############
# Calculate the absolute errors
ANNerrorsM_All_I = abs(ANNpredictionsM_Idaho_All - test_labelsM)
ANNerrorsT_All_I = abs(ANNpredictionsT_Idaho_All - test_labelsT)
ANNerrorsF_All_I = abs(ANNpredictionsF_Idaho_All - test_labelsF)
# Print out the mean absolute error (mae)
print('ANN trained on all and applied to Idaho')
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(ANNerrorsM_All_I), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(ANNerrorsT_All_I), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(ANNerrorsF_All_I), 4))

######### Compute Errors on Test Data #############
# Calculate the absolute errors
errorsM_All_I = abs(predictionsM_Idaho_All - test_labelsM)
errorsT_All_I = abs(predictionsT_Idaho_All - test_labelsT)
errorsF_All_I = abs(predictionsF_Idaho_All - test_labelsF)
# Print out the mean absolute error (mae)
print('Mean Absolute Error for predicting Moisture Scale:', round(np.mean(errorsM_All_I), 4))
print('Mean Absolute Error for predicting Temperature Scale:', round(np.mean(errorsT_All_I), 4))
print('Mean Absolute Error for predicting Friction Velocity:', round(np.mean(errorsF_All_I), 4))
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsT_Idaho_All
y = test_labelsT
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained RF vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_all-on-idaho_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsT_Idaho_All
y = test_labelsT
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained ANN vs Observed Temperature Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_all-on-idaho_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsM_Idaho_All
y = test_labelsM
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained RF vs Observed Moisture Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_all-on-idaho_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsM_Idaho_All
y = test_labelsM
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained ANN vs Observed Moisture Scale for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_all-on-idaho_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsF_Idaho_All
y = test_labelsF
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained RF vs Observed Friction Velocity for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_all-on-idaho_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsF_Idaho_All
y = test_labelsF
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained ANN vs Observed Friction Velocity for Idaho ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_all-on-idaho_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsT_Cabauw_All
y = test_labelsTC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained RF vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_all-on-cabauw_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsT_Cabauw_All
y = test_labelsTC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained ANN Predictions vs Observed Temperature Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_all-on-cabauw_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsM_Cabauw_All
y = test_labelsMC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained RF vs Observed Moisture Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_all-on-cabauw_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsM_Cabauw_All
y = test_labelsMC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained ANN vs Observed Moisture Scale for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_all-on-cabauw_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predictionsF_Cabauw_All
y = test_labelsFC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained RF vs Observed Friction Velocity for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_all-on-cabauw_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredictionsF_Cabauw_All
y = test_labelsFC
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('Idaho + Cabuaw Trained ANN vs Observed Friction Velocity for Cabauw ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_all-on-cabauw_Observed_friction-velocity_30min.png')
plt.show()
# Split results by Bulk Richardson Number
TestCF = datasetTestC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
TestCF['FrictionVelocityTestC'] = FrictionVelocityTestC
TestCM = datasetTestC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
TestCM['MoistureScaleTestC'] = MoistureScaleTestC
TestCT = datasetTestC[['10m Potential Temp K','LowLevel Wind Speed m/s','MidLevel Wind Speed m/s',
             '2m RH %','0m GHI w/m^2','0m Solar Zenith Angle deg','moisture availability_3 cm_',
             '0m Temp K','10m Bulk Richardson']]
TestCT['TemperatureScaleTestC'] = TemperatureScaleTestC

test_featuresMC_Neg = TestCM.loc[TestCM['10m Bulk Richardson']<0]
test_labelsMC_Neg = test_featuresMC_Neg['MoistureScaleTestC']
test_featuresMC_Neg = test_featuresMC_Neg.drop(columns=['MoistureScaleTestC'])
test_featuresMC_Pos = TestCM.loc[TestCM['10m Bulk Richardson']>0]
test_labelsMC_Pos = test_featuresMC_Pos['MoistureScaleTestC']
test_featuresMC_Pos = test_featuresMC_Pos.drop(columns=['MoistureScaleTestC'])
predMC_Neg = rfMC.predict(test_featuresMC_Neg) 
predMC_Pos = rfMC.predict(test_featuresMC_Pos)
ANNpredMC_Neg = ANN_MC.predict(np.asarray(test_featuresMC_Neg)) 
ANNpredMC_Pos = ANN_MC.predict(np.asarray(test_featuresMC_Pos))

ANNpredMC_Neg = ANNpredMC_Neg.reshape(len(ANNpredMC_Neg,))
ANNpredMC_Pos = ANNpredMC_Pos.reshape(len(ANNpredMC_Pos,))


test_featuresTC_Neg = TestCT.loc[TestCT['10m Bulk Richardson']<0]
test_labelsTC_Neg = test_featuresTC_Neg['TemperatureScaleTestC']
test_featuresTC_Neg = test_featuresTC_Neg.drop(columns=['TemperatureScaleTestC'])
test_featuresTC_Pos = TestCT.loc[TestCT['10m Bulk Richardson']>0]
test_labelsTC_Pos = test_featuresTC_Pos['TemperatureScaleTestC']
test_featuresTC_Pos = test_featuresTC_Pos.drop(columns=['TemperatureScaleTestC'])
predTC_Neg = rfTC.predict(test_featuresTC_Neg)  
predTC_Pos = rfTC.predict(test_featuresTC_Pos) 
ANNpredTC_Neg = ANN_TC.predict(np.asarray(test_featuresTC_Neg))  
ANNpredTC_Pos = ANN_TC.predict(np.asarray(test_featuresTC_Pos))  

ANNpredTC_Neg = ANNpredTC_Neg.reshape(len(ANNpredTC_Neg,))
ANNpredTC_Pos = ANNpredTC_Pos.reshape(len(ANNpredTC_Pos,))

test_featuresFC_Neg = TestCF.loc[TestCF['10m Bulk Richardson']<0]
test_labelsFC_Neg = test_featuresFC_Neg['FrictionVelocityTestC']
test_featuresFC_Neg = test_featuresFC_Neg.drop(columns=['FrictionVelocityTestC'])
test_featuresFC_Pos = TestCF.loc[TestCF['10m Bulk Richardson']>0]
test_labelsFC_Pos = test_featuresFC_Pos['FrictionVelocityTestC']
test_featuresFC_Pos = test_featuresFC_Pos.drop(columns=['FrictionVelocityTestC'])
predFC_Neg = rfFC.predict(test_featuresFC_Neg.iloc[:,0:9])  
predFC_Pos = rfFC.predict(test_featuresFC_Pos.iloc[:,0:9])  
ANNpredFC_Neg = ANN_FC.predict(np.asarray(test_featuresFC_Neg))  
ANNpredFC_Pos = ANN_FC.predict(np.asarray(test_featuresFC_Pos))  

ANNpredFC_Neg = ANNpredFC_Neg.reshape(len(ANNpredFC_Neg,))
ANNpredFC_Pos = ANNpredFC_Pos.reshape(len(ANNpredFC_Pos,))

## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predFC_Pos
y = test_labelsFC_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF vs Observed Friction Velocity for Cabauw +BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_cabauw_pos_BRN_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredFC_Pos
y = test_labelsFC_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN vs Observed Friction Velocity for Cabauw +BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_cabauw_pos_BRN_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predFC_Neg
y = test_labelsFC_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF vs Observed Friction Velocity for Cabauw -BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_cabauw_neg_BRN_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredFC_Neg
y = test_labelsFC_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN vs Observed Friction Velocity for Cabauw -BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_cabauw_neg_BRN_Observed_friction-velocity_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predMC_Pos
y = test_labelsMC_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF vs Observed Moisture Scale for Cabauw +BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_cabauw_pos_BRN_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredMC_Pos
y = test_labelsMC_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN vs Observed Moisture Scale for Cabauw +BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_cabauw_pos_BRN_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predMC_Neg
y = test_labelsMC_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF vs Observed Moisture Scale for Cabauw -BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_cabauw_neg_BRN_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredMC_Neg
y = test_labelsMC_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN vs Observed Moisture Scale for Cabauw -BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_cabauw_neg_BRN_Observed_mscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predTC_Pos
y = test_labelsTC_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF vs Observed Temperature Scale for Cabauw +BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_cabauw_pos_BRN_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredTC_Pos
y = test_labelsTC_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN vs Observed Temperature Scale for Cabauw +BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_cabauw_pos_BRN_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predTC_Neg
y = test_labelsTC_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('RF vs Observed Temperature Scale for Cabauw -BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('RF_cabauw_neg_BRN_Observed_tscale_30min.png')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = ANNpredTC_Neg
y = test_labelsTC_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
print("The R^2 value is...",r2)
plt.axhline(0,color='gray',linestyle = '--') # x = 0
plt.axvline(0,color='gray',linestyle = '--') # y = 0
plt.title('ANN vs Observed Temperature Scale for Cabauw -BRN ($\\theta *$)', size = 24)
plt.xlabel('Predicted', size = 24)
plt.ylabel('Observed', size = 24)
#plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.savefig('ANN_cabauw_neg_BRN_Observed_tscale_30min.png')
plt.show()
# Split results by Bulk Richardson Number
TestF = test_featuresF
TestF['FrictionVelocityTest'] = FrictionVelocityTest
TestM = test_featuresM
TestM['MoistureScaleTest'] = MoistureScaleTest
TestT = test_featuresT
TestT['TemperatureScaleTest'] = TemperatureScaleTest

test_featuresM_Neg = TestM.loc[TestM['10m Bulk Richardson']<0]

test_labelsM_Neg = test_featuresM_Neg['MoistureScaleTest']
test_featuresM_Neg = test_featuresM_Neg.drop(columns=['MoistureScaleTest'])
test_featuresM_Pos = TestM.loc[TestM['10m Bulk Richardson']>0]

test_labelsM_Pos = test_featuresM_Pos['MoistureScaleTest']
test_featuresM_Pos = test_featuresM_Pos.drop(columns=['MoistureScaleTest'])
predM_Neg = rfM.predict(test_featuresM_Neg) 
predM_Pos = rfM.predict(test_featuresM_Pos)

test_featuresT_Neg = TestT.loc[TestT['10m Bulk Richardson']<0]

test_labelsT_Neg = test_featuresT_Neg['TemperatureScaleTest']
test_featuresT_Neg = test_featuresT_Neg.drop(columns=['TemperatureScaleTest'])
test_featuresT_Pos = TestT.loc[TestT['10m Bulk Richardson']>0]

test_labelsT_Pos = test_featuresT_Pos['TemperatureScaleTest']
test_featuresT_Pos = test_featuresT_Pos.drop(columns=['TemperatureScaleTest'])
predT_Neg = rfT.predict(test_featuresT_Neg)  
predT_Pos = rfT.predict(test_featuresT_Pos)  

test_featuresF_Neg = TestF.loc[TestF['10m Bulk Richardson']<0]

test_labelsF_Neg = test_featuresF_Neg['FrictionVelocityTest']
test_featuresF_Neg = test_featuresF_Neg.drop(columns=['FrictionVelocityTest'])
test_featuresF_Pos = TestF.loc[TestF['10m Bulk Richardson']>0]

test_labelsF_Pos = test_featuresF_Pos['FrictionVelocityTest']
test_featuresF_Pos = test_featuresF_Pos.drop(columns=['FrictionVelocityTest'])
predF_Neg = rfF.predict(test_featuresF_Neg)  
predF_Pos = rfF.predict(test_featuresF_Pos)  
# Split results by Bulk Richardson Number
TestF = test_featuresF
TestF['FrictionVelocityTest'] = FrictionVelocityTest
TestM = test_featuresM
TestM['MoistureScaleTest'] = MoistureScaleTest
TestT = test_featuresT
TestT['TemperatureScaleTest'] = TemperatureScaleTest

test_featuresM_Neg = TestM.loc[TestM['10m Bulk Richardson']<0]

#df = df[df.hc == 2]
tstar_in2 = pd.DataFrame(tstar_in)

tstar_in_Neg = tstar_in2[TestM['10m Bulk Richardson']<0]
test_labelsM_Neg = test_featuresM_Neg['MoistureScaleTest']
test_featuresM_Neg = test_featuresM_Neg.drop(columns=['MoistureScaleTest'])
test_featuresM_Pos = TestM.loc[TestM['10m Bulk Richardson']>0]
tstar_in_Pos = tstar_in.loc[TestM['10m Bulk Richardson']>0]
test_labelsM_Pos = test_featuresM_Pos['MoistureScaleTest']
test_featuresM_Pos = test_featuresM_Pos.drop(columns=['MoistureScaleTest'])
predM_Neg = rfM.predict(test_featuresM_Neg) 
predM_Pos = rfM.predict(test_featuresM_Pos)

test_featuresT_Neg = TestT.loc[TestT['10m Bulk Richardson']<0]

test_labelsT_Neg = test_featuresT_Neg['TemperatureScaleTest']
test_featuresT_Neg = test_featuresT_Neg.drop(columns=['TemperatureScaleTest'])
test_featuresT_Pos = TestT.loc[TestT['10m Bulk Richardson']>0]

test_labelsT_Pos = test_featuresT_Pos['TemperatureScaleTest']
test_featuresT_Pos = test_featuresT_Pos.drop(columns=['TemperatureScaleTest'])
predT_Neg = rfT.predict(test_featuresT_Neg)  
predT_Pos = rfT.predict(test_featuresT_Pos)  

test_featuresF_Neg = TestF.loc[TestF['10m Bulk Richardson']<0]

test_labelsF_Neg = test_featuresF_Neg['FrictionVelocityTest']
test_featuresF_Neg = test_featuresF_Neg.drop(columns=['FrictionVelocityTest'])
test_featuresF_Pos = TestF.loc[TestF['10m Bulk Richardson']>0]

test_labelsF_Pos = test_featuresF_Pos['FrictionVelocityTest']
test_featuresF_Pos = test_featuresF_Pos.drop(columns=['FrictionVelocityTest'])
predF_Neg = rfF.predict(test_featuresF_Neg)  
predF_Pos = rfF.predict(test_featuresF_Pos)  
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predF_Pos
y = test_labelsF_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
plt.text(0.5,1.0,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R value is...",r_value)
plt.title('RF Predictions vs Observed Friction Velocity for Idaho +BRN ($\\theta *$)', size = 20)
plt.xlabel('Predicted', size = 20)
plt.ylabel('Observed', size = 20)
plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predF_Neg
y = test_labelsF_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
plt.text(0.5,1.2,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R value is...",r_value)
plt.title('RF Predictions vs Observed Friction Velocity for Idaho -BRN ($\\theta *$)', size = 20)
plt.xlabel('Predicted', size = 20)
plt.ylabel('Observed', size = 20)
plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predM_Pos
y = test_labelsM_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
plt.text(0,-0.15,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R value is...",r_value)
plt.title('RF Predictions vs Observed Moisture Scale for Idaho +BRN ($\\theta *$)', size = 20)
plt.xlabel('Predicted', size = 20)
plt.ylabel('Observed', size = 20)
plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predM_Neg
y = test_labelsM_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
plt.text(0.5,0,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R value is...",r_value)
plt.title('RF Predictions vs Observed Moisture Scale for Idaho -BRN ($\\theta *$)', size = 20)
plt.xlabel('Predicted', size = 20)
plt.ylabel('Observed', size = 20)
plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predT_Pos
y = test_labelsT_Pos
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
plt.text(-0.1,0.8,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R value is...",r_value)
plt.title('RF Predictions vs Observed Temperature Scale for Idaho +BRN ($\\theta *$)', size = 20)
plt.xlabel('Predicted', size = 20)
plt.ylabel('Observed', size = 20)
plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.show()
## Plots below after flux and met data are joined together...
from matplotlib.colors import LogNorm
from scipy import stats
plt.rcParams['figure.figsize'] = [10, 10]
x = predT_Neg
y = test_labelsT_Neg
x = x[y>-5] # Filter out a bad data point
y = y[y>-5] # Filter out a bad data point
plt.hist2d(x, y, bins=100, cmin=1, norm=LogNorm())
slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
line = slope*x+intercept
r2 = np.corrcoef(x, y)[0, 1] ** 2
plt.text(-0.3,-2,"$R^2$: {0:0.2f}".format(r2),fontsize=16)
print("The R value is...",r_value)
plt.title('RF Predictions vs Observed Temperature Scale for Idaho -BRN ($\\theta *$)', size = 20)
plt.xlabel('Predicted', size = 20)
plt.ylabel('Observed', size = 20)
plt.plot(np.linspace(2, -15), np.linspace(2, -15), 'k--')
plt.show()
















# Predict whether the temperature scale is postive or negative, and then train separate Regime-Dependent-RF's (RD-RFs)
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
t_scale_class = RandomForestClassifier(n_estimators = 10, n_jobs=4, max_features="sqrt",
                                       min_samples_split=2,min_samples_leaf=1, class_weight="balanced")
t_scale_class.fit(train_featuresT.iloc[:,0:19], np.where(train_labelsT >= 0, 1, -1))
# Compute the confusion matrix
from sklearn.metrics import mean_absolute_error, r2_score, confusion_matrix
t_scale_sign_preds = t_scale_class.predict(test_featuresT.iloc[:,0:19])
test_sign_labels = np.where(test_labelsT >= 0, 1, -1)
print("Confusion matrix comparing test data labels and predictions")
confusion_matrix(test_sign_labels, t_scale_sign_preds)
# Train random forest on negative temperature scales
t_scale_neg_reg = RandomForestRegressor(n_estimators= 50, random_state=2, n_jobs=4,criterion='mse',
                                        min_samples_split=4,min_samples_leaf=2)
neg_t_scale = train_labelsT < 0
train_data = train_featuresT.iloc[:,0:19]
input_cols = featuresF
t_scale_neg_reg.fit(train_data.loc[neg_t_scale, input_cols], 
                    np.log(-train_labelsT[neg_t_scale]))
# Train random forest on positive temperature scales
t_scale_pos_reg = RandomForestRegressor(n_estimators= 50, random_state=2, n_jobs=4,criterion='mse',
                                        min_samples_split=4,min_samples_leaf=2)
pos_t_scale = train_labelsT > 0
t_scale_pos_reg.fit(train_data.loc[pos_t_scale, input_cols], 
                    np.log(train_labelsT[pos_t_scale]))
# Make predictions on test data
test_data = test_featuresF.iloc[:,0:19]
t_scale_pos_preds = t_scale_pos_reg.predict(test_data[input_cols])
t_scale_neg_preds = t_scale_neg_reg.predict(test_data[input_cols])
# Combine both positive and negative predictions
t_scale_combined_preds = np.where(t_scale_sign_preds > 0, np.exp(t_scale_pos_preds), -np.exp(t_scale_neg_preds))
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.hist2d(t_scale_combined_preds, test_labelsT, cmin=1, bins=np.linspace(-2,1,31), norm=LogNorm())
plt.xlim(-2, 1)
plt.ylim(-2, 1)
plt.plot(np.linspace(-2, 1), np.linspace(-2, 1), 'k--')
plt.plot(np.zeros(50), np.linspace(-2, 1), 'k--')
plt.plot(np.linspace(-2, 1), np.zeros(50), 'k--')
plt.colorbar()
plt.title("Multi-Stage Random Forest", fontsize=14)
plt.ylabel("Observed Temperature Scale", fontsize=14)
mae_rf_multi = mean_absolute_error(test_labelsT, t_scale_combined_preds)
r2_rf_multi = r2_score(test_labelsT, t_scale_combined_preds)
plt.xlabel("MAE: {0:0.4f} R$^2$: {1:0.3f}".format(mae_rf_multi, r2_rf_multi), fontsize=14)
plt.subplot(1, 2, 2)
rf_t_scale_preds = predictionsT
plt.hist2d(rf_t_scale_preds, test_labelsT, cmin=1, bins=np.linspace(-2,1,31), norm=LogNorm())
plt.xlim(-2, 1)
plt.ylim(-2, 1)
plt.plot(np.linspace(-2, 1), np.linspace(-2, 1), 'k--')
plt.plot(np.zeros(50), np.linspace(-2, 1), 'k--')
plt.plot(np.linspace(-2, 1), np.zeros(50), 'k--')

plt.title("Single Random Forest", fontsize=14)
plt.ylabel("Observed Temperature Scale", fontsize=14)
plt.colorbar()

mae_rf_single = mean_absolute_error(test_labelsT, rf_t_scale_preds)
r2_rf_single = r2_score(test_labelsT, rf_t_scale_preds)
plt.xlabel("MAE: {0:0.4f} R$^2$: {1:0.3f}".format(mae_rf_single, r2_rf_single), fontsize=14)
# Predict whether the temperature scale is postive or negative, and then train separate Regime-Dependent-RF's (RD-RFs)
m_scale_class = RandomForestClassifier(n_estimators = 50, n_jobs=4, max_features="sqrt",
                                       min_samples_split=2,min_samples_leaf=1, class_weight="balanced")
m_scale_class.fit(train_featuresM.iloc[:,0:19], np.where(train_labelsM >= 0, 1, -1))

# Compute the confusion matrix
m_scale_sign_preds = m_scale_class.predict(test_featuresM.iloc[:,0:19])
test_sign_labels = np.where(test_labelsM >= 0, 1, -1)
print("Confusion matrix comparing test data labels and predictions")
print(confusion_matrix(test_sign_labels, m_scale_sign_preds))

# Train random forest on negative temperature scales
m_scale_neg_reg = RandomForestRegressor(n_estimators= 50, random_state=2, n_jobs=4,criterion='mse',
                                        min_samples_split=4,min_samples_leaf=2)
neg_m_scale = train_labelsM < 0
train_data = train_featuresM.iloc[:,0:19]
input_cols = featuresF
m_scale_neg_reg.fit(train_data.loc[neg_m_scale, input_cols], 
                    np.log(-train_labelsM[neg_m_scale]))

# Train random forest on positive temperature scales
m_scale_pos_reg = RandomForestRegressor(n_estimators= 50, random_state=2, n_jobs=4,criterion='mse',
                                        min_samples_split=4,min_samples_leaf=2)
pos_m_scale = train_labelsM > 0
m_scale_pos_reg.fit(train_data.loc[pos_m_scale, input_cols], 
                    np.log(train_labelsM[pos_m_scale]))

# Make predictions on test data
test_data = test_featuresM.iloc[:,0:19]
m_scale_pos_preds = m_scale_pos_reg.predict(test_data[input_cols])
m_scale_neg_preds = m_scale_neg_reg.predict(test_data[input_cols])
# Combine both positive and negative predictions
m_scale_combined_preds = np.where(m_scale_sign_preds > 0, np.exp(m_scale_pos_preds), -np.exp(m_scale_neg_preds))

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.hist2d(m_scale_combined_preds, test_labelsM, cmin=1, bins=np.linspace(-2,1,31), norm=LogNorm())
plt.xlim(-2, 1)
plt.ylim(-2, 1)
plt.plot(np.linspace(-2, 1), np.linspace(-2, 1), 'k--')
plt.plot(np.zeros(50), np.linspace(-2, 1), 'k--')
plt.plot(np.linspace(-2, 1), np.zeros(50), 'k--')
plt.colorbar()
plt.title("Multi-Stage Random Forest", fontsize=14)
plt.ylabel("Observed Temperature Scale", fontsize=14)
mae_rf_multi = mean_absolute_error(test_labelsM, m_scale_combined_preds)
r2_rf_multi = r2_score(test_labelsM, m_scale_combined_preds)
plt.xlabel("MAE: {0:0.4f} R$^2$: {1:0.3f}".format(mae_rf_multi, r2_rf_multi), fontsize=14)
plt.subplot(1, 2, 2)
rf_m_scale_preds = predictionsM
plt.hist2d(rf_m_scale_preds, test_labelsM, cmin=1, bins=np.linspace(-2,1,31), norm=LogNorm())
plt.xlim(-2, 1)
plt.ylim(-2, 1)
plt.plot(np.linspace(-2, 1), np.linspace(-2, 1), 'k--')
plt.plot(np.zeros(50), np.linspace(-2, 1), 'k--')
plt.plot(np.linspace(-2, 1), np.zeros(50), 'k--')

plt.title("Single Random Forest", fontsize=14)
plt.ylabel("Observed Temperature Scale", fontsize=14)
plt.colorbar()

mae_rf_single = mean_absolute_error(test_labelsM, rf_m_scale_preds)
r2_rf_single = r2_score(test_labelsM, rf_m_scale_preds)
plt.xlabel("MAE: {0:0.4f} R$^2$: {1:0.3f}".format(mae_rf_single, r2_rf_single), fontsize=14)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15166795965865013, 'fraction_numerical': 0.03491921610955577, 'mean_word_length': 4.847872183046501, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 122, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '539611', 'license': 'MIT', 'n_tokens_mistral': 48783, 'n_tokens_neox': 44770, 'n_words': 10361, 'path': 'notebooks/ANN_RF_Final-30min_averages_WRF_predictors.ipynb', 'repo_name': 'NCAR/mlsurfacelayer', 'size': '118604'}