Using the Baum Ice Crystal Database with HR

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import sasktran as sk

SASKTRAN supports calculations using the Baum ice crystal database. The core of this functionality is contained within the BaumIceCrystal optical property. In principle including ice crystals is as simple as creating a climatology to go along with the optical property and adding it to the engine, however complications involving the radiative transfer of ice clouds usually means that the HR engine must be configured in a specific way to obtain accurate calculations.

The SpeciesBaumIceCloud species combines the BaumIceCrystal optical property with a climatology that creates a Gaussian shaped cloud with a specified size and optical thickness.

baum_species = sk.SpeciesBaumIceCloud(particlesize_microns=20,
OSError                                   Traceback (most recent call last)
Cell In[2], line 1
----> 1 baum_species = sk.SpeciesBaumIceCloud(particlesize_microns=20, 
      2                                       cloud_top_m=16500,
      3                                       cloud_width_fwhm_m=500,
      4                                       vertical_optical_depth=0.01,
      5                                       vertical_optical_depth_wavel_nm=750)

File /usr/share/miniconda/envs/doc_env/lib/python3.9/site-packages/sasktran/, in SpeciesBaumIceCloud.__init__(self, particlesize_microns, cloud_top_m, cloud_width_fwhm_m, vertical_optical_depth, vertical_optical_depth_wavel_nm, altitude_resolution_m, num_sigma, name, use_delta_eddington)
    251 def __init__(self, particlesize_microns: float, cloud_top_m: float,
    252              cloud_width_fwhm_m: float, vertical_optical_depth: float,
    253              vertical_optical_depth_wavel_nm: float, altitude_resolution_m: float=10, num_sigma: int=5,
    254              name='icecloud', use_delta_eddington=True):
--> 255     optical_property = sk.BaumIceCrystal(particlesize_microns, use_delta_eddington=use_delta_eddington)
    257     cloud_sigma = cloud_width_fwhm_m / (2 * np.sqrt(2 * np.log(2)))
    259     # Cloud top is defined as middle + 2 sigma

File /usr/share/miniconda/envs/doc_env/lib/python3.9/site-packages/sasktran/, in BaumIceCrystal.__init__(self, effective_size_microns, use_delta_eddington)
    857 def __init__(self, effective_size_microns, use_delta_eddington=True):
    858     super().__init__('BAUM_ICECRYSTALS')
--> 859     self._validate_registry()
    861     self._effective_size_microns = effective_size_microns
    862     self._use_delta_eddington = use_delta_eddington

File /usr/share/miniconda/envs/doc_env/lib/python3.9/site-packages/sasktran/, in BaumIceCrystal._validate_registry()
    884 if not os.path.isfile(os.path.join(baum_folder, '')):
    885     user_config = config.user_config_file_location()
--> 886     raise IOError(
    887         'Baum Folder is invalid, please add baum_directory: xxxx to the file {}'.format(user_config))

OSError: Baum Folder is invalid, please add baum_directory: xxxx to the file /home/runner/.config/sasktran/config.yml

Which creates a cloud species that consists of 20 micron sized particles, a cloud top altitude of 16.5 km, a full width at half maximum in altitude of 500 m, and has a vertical optical depth of 0.01 at 750 nm. Note that limb path lengths can be ~100 times that of vertical path lengths, so even a vertically thin cloud may appear thick in the limb geometry. Shown below is the calculated Gaussian number density.

altitudes = np.arange(10000, 20000, 10)

cloud_dens = baum_species.climatology.get_parameter('icecloud', latitude=0, longitude=0, mjd=54372, altitudes=altitudes)

plt.plot(cloud_dens, altitudes)
plt.ylabel('Altitude [m]')
plt.xlabel('Number Density [/cm3]')
NameError                                 Traceback (most recent call last)
Cell In[3], line 3
      1 altitudes = np.arange(10000, 20000, 10)
----> 3 cloud_dens = baum_species.climatology.get_parameter('icecloud', latitude=0, longitude=0, mjd=54372, altitudes=altitudes)
      5 plt.plot(cloud_dens, altitudes)
      6 plt.ylabel('Altitude [m]')

NameError: name 'baum_species' is not defined

We can now set up the HR engine and add the species to it as we would any other species.

tanalts_km = np.arange(10, 50, 1)

# First recreate our geometry and atmosphere classes
geometry = sk.VerticalImage()
geometry.from_sza_saa(sza=60, saa=60, lat=0, lon=0, tanalts_km=tanalts_km, mjd=54372, locallook=0,
                      satalt_km=600, refalt_km=20)

atmosphere = sk.Atmosphere()

atmosphere['ozone'] = sk.Species(sk.O3OSIRISRes(), sk.Labow())
atmosphere['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())
atmosphere['cloud'] = baum_species

# And now make the engine
engine = sk.EngineHR(geometry=geometry, atmosphere=atmosphere)

# Choose some wavelengths to do the calculation at
engine.wavelengths = [340, 800]
NameError                                 Traceback (most recent call last)
Cell In[4], line 12
     10 atmosphere['ozone'] = sk.Species(sk.O3OSIRISRes(), sk.Labow())
     11 atmosphere['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())
---> 12 atmosphere['cloud'] = baum_species
     14 # And now make the engine
     15 engine = sk.EngineHR(geometry=geometry, atmosphere=atmosphere)

NameError: name 'baum_species' is not defined

But before we run the calculation we should modify some settings in the engine to account for the cloud accurately.

Setting the model parameters

The default vertical resolution of SASKTRAN-HR is ~1000 m, however clouds typically have fine structures and the cloud included above is calculated at a 10 m resolution by default. We could set the grid spacing of HR to be 10 m, however this would result in a very slow calculation. It also is not necessary because the cloud is localized over a few km, we do not need 10 m resolution throughout the entire atmosphere. There are also several settings related to the adaptive integration of rays that should be used when including clouds. For convenience the configure_for_cloud method will set most of these things.

engine.configure_for_cloud(cloud_altitudes=baum_species.altitude_grid, grid_spacing_m=1000)

rad = engine.calculate_radiance()
plt.plot(rad.T, tanalts_km)
plt.ylabel('Altitude [km]')
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 engine.configure_for_cloud(cloud_altitudes=baum_species.altitude_grid, grid_spacing_m=1000)
      3 rad = engine.calculate_radiance()
      4 plt.plot(rad.T, tanalts_km)

NameError: name 'engine' is not defined