Merge pull request #87 from jmcook1186/add-easy-albedo-func

Add get_albedo() function

app/streamlit/app.py

@@ -95,7 +95,6 @@ def run_snicar(
         layer = 1
     else:
         raise ValueError("invalid layer type")
-
 
     # first build classes from config file and validate their contents
     (
@@ -121,7 +120,7 @@ def run_snicar(
     illumination.solzen = solar_zenith_angle
     illumination.calculate_irradiance()
     # validate inputs to ensure no invalid combinations have been chosen
-    status = validate_inputs(ice, rt_config, model_config, illumination, impurities)
+    status = validate_inputs(ice, illumination, impurities)
 
     # now get the optical properties of the ice column
     ssa_snw, g_snw, mac_snw = get_layer_OPs(ice, model_config)

src/biosnicar/adding_doubling_solver.py

@@ -33,7 +33,6 @@
 from scipy.signal import savgol_filter
 
 from biosnicar.classes import Outputs
-from biosnicar.setup_snicar import build_impurities_array
 
 
 def adding_doubling_solver(tau, ssa, g, L_snw, ice, illumination, model_config):
@@ -108,11 +107,9 @@ def adding_doubling_solver(tau, ssa, g, L_snw, ice, illumination, model_config):
     # Delta-Eddington computation for that layer is done.
 
     for lyr in np.arange(0, ice.nbr_lyr, 1):  # loop through layers
-
         # condition: if current layer is above fresnel layer or the
         # top layer is a Fresnel layer
         if lyr < lyrfrsnl:
-
             mu0n = mu0
 
         else:
@@ -147,7 +144,6 @@ def adding_doubling_solver(tau, ssa, g, L_snw, ice, illumination, model_config):
         )
 
         if lyr == lyrfrsnl:
-
             (
                 rdif_a,
                 rdif_b,
@@ -301,7 +297,8 @@ def calc_reflectivity_transmittivity(
 
     # coefficient for delta eddington solution for all layers
     # Eq. 50: Briegleb and Light 2007
-    ts = (1 - (wtot * ftot)) * tautot  # layer delta-scaled extinction optical depth
+    # layer delta-scaled extinction optical depth
+    ts = (1 - (wtot * ftot)) * tautot
     ws = ((1 - ftot) * wtot) / (
         1 - (wtot * ftot)
     )  # layer delta-scaled single scattering albedo
@@ -324,7 +321,8 @@ def calc_reflectivity_transmittivity(
     rdif_a[:, lyr] = (
         (ue**2 - 1) * (1 / extins - extins) / ne
     )  # R BAR = layer reflectivity to DIFFUSE radiation
-    tdif_a[:, lyr] = 4 * ue / ne  # T BAR layer transmissivity to DIFFUSE radiation
+    # T BAR layer transmissivity to DIFFUSE radiation
+    tdif_a[:, lyr] = 4 * ue / ne
 
     # evaluate rdir, tdir for direct beam
     trnlay[:, lyr] = np.maximum(
@@ -784,7 +782,6 @@ def calc_correction_fresnel_layer(
     critical_angle = np.arcsin(ref_indx)
 
     for wl in np.arange(0, model_config.nbr_wvl, 1):
-
         if np.arccos(illumination.mu_not) < critical_angle[wl]:
             # in this case, no total internal reflection
 
@@ -825,7 +822,8 @@ def calc_correction_fresnel_layer(
         # Eq. 25  Brigleb and light 2007
         # diffuse reflection of flux arriving from above
 
-        Rf_dif_a = ice.fl_r_dif_a[wl]  # reflection from diffuse unpolarized radiation
+        # reflection from diffuse unpolarized radiation
+        Rf_dif_a = ice.fl_r_dif_a[wl]
         Tf_dif_a = 1 - Rf_dif_a  # transmission from diffuse unpolarized radiation
 
         # diffuse reflection of flux arriving from below
@@ -837,7 +835,8 @@ def calc_correction_fresnel_layer(
         # the fresnel (refractive) layer, always taken to be above
         # the present layer lyr (i.e. be the top interface):
 
-        rintfc = 1 / (1 - Rf_dif_b * rdif_a[wl, lyr])  # denom interface scattering
+        # denom interface scattering
+        rintfc = 1 / (1 - Rf_dif_b * rdif_a[wl, lyr])
 
         # layer transmissivity to DIRECT radiation
         # Eq. B7  Briegleb & Light 2007
@@ -914,7 +913,6 @@ def calc_reflection_below(
     for lyr in np.arange(
         ice.nbr_lyr - 1, -1, -1
     ):  # starts at the bottom and works its way up to the top layer
-
         # Eq. B5  Briegleb and Light 2007
         # interface scattering
         refkp1 = 1 / (1 - rdif_b[:, lyr] * rupdif[:, lyr + 1])
@@ -957,7 +955,6 @@ def trans_refl_at_interfaces(
     dfdir,
     dfdif,
 ):
-
     """Calculates transmission and reflection at layer interfaces.
 
     Args:
@@ -987,7 +984,6 @@ def trans_refl_at_interfaces(
     puny = 1e-10  # not sure how should we define this
 
     for lyr in np.arange(0, ice.nbr_lyr + 1, 1):
-
         # Eq. 52  Briegleb and Light 2007
         # interface scattering
         refk = 1 / (1 - rdndif[:, lyr] * rupdif[:, lyr])
@@ -1025,7 +1021,6 @@ def trans_refl_at_interfaces(
         )
 
         if np.max(dfdir[:, lyr]) < puny:
-
             dfdir[:, lyr] = np.zeros(
                 (model_config.nbr_wvl,), dtype=int
             )  # echmod necessary?
@@ -1034,10 +1029,9 @@ def trans_refl_at_interfaces(
         dfdif[:, lyr] = trndif[:, lyr] * (1 - rupdif[:, lyr]) * refk
 
         if np.max(dfdif[:, lyr]) < puny:
-
             dfdif[:, lyr] = np.zeros(
                 (model_config.nbr_wvl,), dtype=int
-            )  #!echmod necessary?
+            )  # !echmod necessary?
 
     return fdirup, fdifup, fdirdn, fdifdn
 
@@ -1081,7 +1075,6 @@ def calculate_fluxes(
     """
 
     for n in np.arange(0, ice.nbr_lyr + 1, 1):
-
         F_up[:, n] = (
             fdirup[:, n] * (illumination.Fs * illumination.mu_not * np.pi)
             + fdifup[:, n] * illumination.Fd
@@ -1115,7 +1108,6 @@ def calculate_fluxes(
     F_btm_net = -F_net[:, ice.nbr_lyr]
 
     for i in np.arange(0, ice.nbr_lyr, 1):  # [0,1,2,3,4]
-
         F_abs_vis[i] = sum(F_abs[0 : model_config.vis_max_idx, i])
         F_abs_nir[i] = sum(
             F_abs[model_config.vis_max_idx : model_config.nir_max_idx, i]
@@ -1152,7 +1144,6 @@ def conservation_of_energy_check(illumination, F_abs, F_btm_net, F_top_pls):
     energy_conservation_error = sum(abs(energy_sum))
 
     if energy_conservation_error > 1e-10:
-
         raise ValueError(f"energy conservation error: {energy_conservation_error}")
     else:
         pass
@@ -1177,7 +1168,8 @@ def get_outputs(illumination, albedo, model_config, L_snw, F_abs, F_btm_net):
 
     # Radiative heating rate:
     F_abs_slr = np.sum(F_abs, axis=0)
-    heat_rt = F_abs_slr / (L_snw * 2117)  # [K/s] 2117 = specific heat ice (J kg-1 K-1)
+    # [K/s] 2117 = specific heat ice (J kg-1 K-1)
+    heat_rt = F_abs_slr / (L_snw * 2117)
     outputs.heat_rt = heat_rt * 3600  # [K/hr]
 
     # Spectral albedo

src/biosnicar/classes.py

@@ -21,7 +21,6 @@
 """
 import math
 import os
-
 import numpy as np
 import xarray as xr
 import yaml
@@ -47,7 +46,6 @@ class Impurity:
     """
 
     def __init__(self, file, coated, cfactor, unit, name, conc):
-
         self.name = name
         self.cfactor = cfactor
         self.unit = unit
@@ -94,7 +92,6 @@ class Ice:
     """
 
     def __init__(self, input_file):
-
         # use config to calculate refractive indices
         with open(input_file, "r") as ymlfile:
             inputs = yaml.load(ymlfile, Loader=yaml.FullLoader)
@@ -178,7 +175,6 @@ class Illumination:
     """
 
     def __init__(self, input_file):
-
         with open(input_file, "r") as ymlfile:
             inputs = yaml.load(ymlfile, Loader=yaml.FullLoader)
 
@@ -259,7 +255,6 @@ class RTConfig:
     """
 
     def __init__(self, input_file):
-
         with open(input_file, "r") as ymlfile:
             inputs = yaml.load(ymlfile, Loader=yaml.FullLoader)
 
@@ -289,7 +284,6 @@ class ModelConfig:
     """
 
     def __init__(self, input_file):
-
         with open(input_file, "r") as ymlfile:
             inputs = yaml.load(ymlfile, Loader=yaml.FullLoader)
 
@@ -366,7 +360,6 @@ class PlotConfig:
     """
 
     def __init__(self, input_file):
-
         with open(input_file, "r") as ymlfile:
             inputs = yaml.load(ymlfile, Loader=yaml.FullLoader)
 
@@ -457,7 +450,6 @@ class BioOpticalConfig:
     """
 
     def __init__(self, input_file):
-
         with open(input_file, "r") as ymlfile:
             inputs = yaml.load(ymlfile, Loader=yaml.FullLoader)
 
@@ -501,7 +493,6 @@ def __init__(self, input_file):
         self.validate_biooptical_inputs()
 
     def validate_biooptical_inputs(self):
-
         if self.Mie:
             assert self.GO == False
         if self.GO:

src/biosnicar/column_OPs.py

@@ -515,7 +515,7 @@ def mix_in_impurities(ssa_snw, g_snw, mac_snw, ice, impurities, model_config):
         if impurity.unit == 1:
 
             mss_aer[0 : ice.nbr_lyr, i] = (
-                np.array(impurity.conc) / 917 * 10 ** 6
+                np.array(impurity.conc) / 917 * 10**6
             ) * impurity.cfactor
 
         else:
@@ -575,5 +575,5 @@ def mix_in_impurities(ssa_snw, g_snw, mac_snw, ice, impurities, model_config):
     return tau, ssa, g, L_snw
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     pass

src/biosnicar/display.py

@@ -5,7 +5,6 @@
 
 
 def setup_axes(plot_config):
-
     rc = {
         "figure.figsize": (8, 6),
         "axes.facecolor": str(plot_config.facecolor),
@@ -25,7 +24,6 @@ def setup_axes(plot_config):
 
 
 def plot_albedo(plot_config, model_config, albedo):
-
     rc = setup_axes(plot_config)
     plt.style.use("seaborn")
     sns.set_style("white")
@@ -49,7 +47,6 @@ def plot_albedo(plot_config, model_config, albedo):
 
 
 def display_out_data(outputs):
-
     print("\n** OUTPUT DATA **")
     print("Broadband albedo: ", np.round(outputs.BBA, 4))
 
@@ -66,7 +63,6 @@ def display_out_data(outputs):
 
 
 def calculate_band_ratios(albedo):
-
     I2DBA = albedo[51] / albedo[46]
     I3DBA = (albedo[46] - albedo[50]) / albedo[55]
     NDCI = ((albedo[50] - albedo[48]) - (albedo[55] - albedo[48])) * (

src/biosnicar/geometric_optics_ice.py

@@ -41,10 +41,8 @@
 
 NOTE: The extinction coefficient in the current implementation is 2 for all size
 parameters as assumed in the conventional geometric optics approximation.
-
-
-
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
@@ -57,7 +55,6 @@
 
 
 def preprocess_RI(ri_source, path_to_ri):
-
     """Preprocessing of wavelength and RI data.
 
     Preprocessing function that ensures the wavelengths and real/imaginary
@@ -153,7 +150,6 @@ def calc_optical_params(
     delta = 0.3
 
     for i in np.arange(0, len(wavelengths), 1):
-
         mr = reals[i]
         mi = imags[i]
         wl = wavelengths[i]

src/biosnicar/get_albedo.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from pathlib import Path
+from validate_inputs import validate_inputs
+from adding_doubling_solver import adding_doubling_solver
+from column_OPs import get_layer_OPs, mix_in_impurities
+from display import display_out_data, plot_albedo
+from setup_snicar import setup_snicar
+from toon_rt_solver import toon_solver
+
+
+def get_albedo(solver, plot, validate):
+    (
+        ice,
+        illumination,
+        rt_config,
+        model_config,
+        plot_config,
+        impurities,
+    ) = setup_snicar("default")
+
+    if validate:
+        validate_inputs(ice, illumination, impurities)
+
+    # now get the optical properties of the ice column
+    ssa_snw, g_snw, mac_snw = get_layer_OPs(ice, model_config)
+    tau, ssa, g, L_snw = mix_in_impurities(
+        ssa_snw, g_snw, mac_snw, ice, impurities, model_config
+    )
+    # now run one or both of the radiative transfer solvers
+    if solver == "toon":
+        print("\nRunning biosnicar with the Toon solver\n")
+        outputs = toon_solver(
+            tau, ssa, g, L_snw, ice, illumination, model_config, rt_config
+        )
+    elif solver == "adding-doubling":
+        print("\nRunning biosnicar with the adding-doubling solver\n")
+        outputs = adding_doubling_solver(
+            tau, ssa, g, L_snw, ice, illumination, model_config
+        )
+    else:
+        return "solver not recognized, please choose toon or adding-doubling"
+
+    if plot:
+        plot_albedo(plot_config, model_config, outputs.albedo)
+    display_out_data(outputs)
+    return outputs.albedo

src/biosnicar/inputs.yaml

@@ -67,7 +67,7 @@ IMPURITIES:
     CFACTOR: 1
     COATED: False
     UNIT: 1
-    CONC: [1000000 ,0]
+    CONC: [0 ,0]
 
 PLOT:
     FIG_SIZE: (8,6)

src/biosnicar/main.py

@@ -0,0 +1,9 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from pathlib import Path
+from get_albedo import get_albedo
+
+# call easy albedo func
+albedo = get_albedo("adding-doubling", plot=True, validate=True)