Bandpass averaged quantities using Python and pandas

The OceanColor website provides a description of how the bandpass averaged quantities used in ocean color data processing were computed. The following is an example of how to compute these using Python and the pandas module. With this example you should be able to replicate the Rayleigh values for VIIRS-SNPP reported at: https://oceancolor.gsfc.nasa.gov/resources/docs/rsr/rsr_tables/

Import a couple modules

import csv
import pandas as pd

Download and read in the required data sets

VIIRS-SNPP spectral response functions

rsr_headers = ["wavelength","M1","M2","M3","M4","M5","M6","M7","M8","M10","M11"]
rsr = pd.read_csv('VIIRSN_IDPSv3_RSRs.txt',names=rsr_headers,skiprows=5, delimiter=' ',skipinitialspace=True,index_col="wavelength")

Bodhaine Rayleigh optical thickness

rayleigh_headers=["wavelength","tau_r","dpol"]
rayleigh = pd.read_csv('rayleigh_bodhaine.txt',names=rayleigh_headers,skiprows=16, delimiter=' ',index_col="wavelength")

Thuillier solar spectrum

solar_headers=["wavelength","irradiance"]
solar = pd.read_csv('f0.txt',names=solar_headers,skiprows=13, delimiter=' ',index_col="wavelength")

Merge the datasets

merged = pd.merge(pd.merge(rsr,solar,right_index=True, left_index=True),rayleigh,right_index=True, left_index=True)

Compute the bandpass averaged Rayleigh optical thickness for band M1 (410nm) using the following equation:

where X is the Rayleigh optical thickness, W is the weighting function (in this case the Thuillier solar spectrum) and RSR is the spectral response function.

merged['numerator'] = merged.apply(lambda row: (row['M1']*row['tau_r']*row['irradiance']), axis=1)
merged['denominator'] = merged.apply(lambda row: (row['M1']*row['irradiance']), axis=1)

print(merged.numerator.sum()/merged.denominator.sum())
0.31756155339404074