Source code for supy.util._ohm

# -*- coding: utf-8 -*-
Created on Mon Dec 17 11:13:44 2018

George Meachim,
Ting Sun

import numpy as np
import pandas as pd

from ._plot import plot_comp, plot_day_clm
from .._env import logger_supy

# Linear fitting of QS, QN, deltaQN/dt (entire timeseries)
[docs]def derive_ohm_coef(ser_QS, ser_QN): """ A function to linearly fit two independant variables to a dependent one. Input params: QS_Ser: QN_Ser: The first independent variable (Net all wave radiation). Pandas Series. Returns: a1, a2 coefficients and a3 (intercept) Parameters ---------- ser_QS : pd.Series The dependent variable QS (Surface heat storage). ser_QN : pd.Series The first independent variable (Net all wave radiation). Returns ------- Tuple a1, a2 coefficients and a3 (intercept) """ from sklearn.linear_model import LinearRegression # derive dt in hours dt_hr = ser_QN.index.freq / pd.Timedelta("1H") # Calculate difference between neighbouring QN values ser_delta_QN_dt = ser_QN.diff() / dt_hr # Drop NaNs and infinite values ser_QS = ser_QS.replace([np.inf, -np.inf], np.nan).dropna(how="all") ser_QN = ser_QN.loc[ser_QS.index] ser_delta_QN_dt = ser_delta_QN_dt.loc[ser_QS.index] # Create DataFrame with regression quantities and rename cols frames = [ser_QS, ser_QN, ser_delta_QN_dt] regression_df = pd.concat(frames, axis=1) regression_df.columns = ["QS", "QN", "delta_QN_dt"] # Reindex after dropping NaNs regression_df.reset_index(drop=True, inplace=True) regression_df.fillna(regression_df.mean(), inplace=True) feature_cols = ["QN", "delta_QN_dt"] X = regression_df[feature_cols].replace([np.inf, -np.inf], np.nan).dropna(how="all") y = regression_df.QS lm = LinearRegression(), y) a1 = lm.coef_[0] a2 = lm.coef_[1] a3 = lm.intercept_ return a1, a2, a3
def replace_ohm_coeffs(df_state, coefs, land_cover_type): """ This function takes as input parameters the model initial state DataFrame, the new ohm coefficients as calculated by performing linear regression on AMF Obs and the land cover type for which they were calculated. Parameters ---------- df_state : pandas.DataFrame Model state dataframe used in supy coefs : tuple new a1, a2, a3 coefficients to replace those in `df_state`; note: 1. the format should be (a1, a3, a3); 2. any of a1/a2/a3 can be either one numeric value or a list-like structure of four values for SW/SD/WW/WD conditions indicated by `ohm_threshsw` and `ohm_threshwd`. land_cover_type : str one of seven SUEWS land cover types, can be any of 1. {"Paved","Bldgs","EveTr","DecTr","Grass","BSoil","Water"} (string case does NOT matter); or 2. an integer between 0 to 6. Returns ------- pandas.DataFrame a DataFrame with updated OHM coefficients. """ land_cover_type_dict = { "paved": 0, "bldgs": 1, "evetr": 2, "dectr": 3, "grass": 4, "bsoil": 5, "water": 6, } try: lc_index = ( land_cover_type_dict.get(land_cover_type.lower()) if isinstance(land_cover_type, str) else land_cover_type ) except: list_lc = list(land_cover_type_dict.keys()) logger_supy.error( f"land_cover_type must be one of {list_lc}, instead of {land_cover_type}" ) else: # Instantiate 4x3 matrix of zeros to put old coeffs coef_matrix = np.zeros((4, 3)) coef_matrix[:, 0] = coefs[0] coef_matrix[:, 1] = coefs[1] coef_matrix[:, 2] = coefs[2] # Copy ohm_coef part of df_state_init df_ohm = df_state.loc[:, "ohm_coef"].copy() # Reshape values into matrix form values_ohm = df_ohm.values.reshape((-1, 8, 4, 3)) # Get ohm values corresponding to user specified land cover and assign to matrix values_ohm[:, lc_index] = coef_matrix # Place new ohm values into df_ohm df_ohm.loc[:, :] = values_ohm.reshape(df_ohm.shape) # Make copy of df_state_init df_state_init_copy = df_state.copy() # Replace ohm_coef part of df_state_init with new values df_state_init_copy.loc[:, "ohm_coef"] = df_ohm.values return df_state_init_copy
[docs]def sim_ohm(ser_qn: pd.Series, a1: float, a2: float, a3: float) -> pd.Series: """Calculate QS using OHM (Objective Hysteresis Model). Parameters ---------- ser_qn : pd.Series net all-wave radiation. a1 : float a1 of OHM coefficients. a2 : float a2 of OHM coefficients. a3 : float a3 of OHM coefficients. Returns ------- pd.Series heat storage flux calculated by OHM. """ # derive delta t in hour try: dt_hr = ser_qn.index.freq / pd.Timedelta("1h") except AttributeError as ex: raise RuntimeError("frequency info is missing from input `ser_qn`") from ex # Calculate rate of change of Net All-wave radiation ser_qn_dt = ser_qn.diff() / dt_hr # Derive QS from OBS quantities ser_qs = a1 * ser_qn + a2 * ser_qn_dt + a3 return ser_qs
def compare_heat_storage(ser_qn_obs, ser_qs_obs, a1, a2, a3): """This function compares the storage heat flux calculated with AMF QN and linear regression coefficients with that output by SUEWS. Input params: QN_Ser_Obs: A series of OBS net all-wave radiation. QS_Ser_SUEWS: A series of SUEWS storage heat flux values. a1, a2, a3: Linear regression coefficients from ohm_linregress_clean. Returns: MPL plot of diurnal comparison, MPL plot with 1:1 line and fitted line. """ # calculate qs using OHM ser_qs_sim = sim_ohm(ser_qn_obs, a1, a2, a3).rename("Sim") # re-organise obs and sim info one dataframe ser_qs_obs = ser_qs_obs.rename("Obs") df_qs = pd.concat([ser_qs_sim, ser_qs_obs.rename("Obs")], axis=1) # Plotting plot1 = plot_day_clm(df_qs) plot2 = plot_comp(df_qs) return plot1, plot2