Details of the pycoxmunk algorithm
Introduction
pycoxmunk uses a modified version of the mathematical equations described by Cox and Munk in their two 1954 papers describing the roughness and statistics of the sea surface derived from photographs of sunglint. A more detailed description of the method as-applied to satellite imagery is given in Sayer et al, 2010
The algorithm described here computes both the per-wavelength sea surface reflectance, \(\rho\), and the four bidirectional reflectance terms: - \(\rho_{0v}\): Solar beam to satellite view reflectances - \(\rho_{0d}\): Solar beam to diffuse reflectances - \(\rho_{dv}\): Diffuse to satellite view reflectances - \(\rho_{dd}\): Diffuse to diffuse reflectances
To estimate the reflectances, pycoxmunk requires knowledge of the sun-satellite viewing geometry through the following angles:
\(\theta_v\): The satellite zenith angle, defined as 0° when the satellite is directly overhead and 90° when the satellite is at the horizon.
\(\phi_v\): The satellite azimuth angle, defined as 0° when the satellite is directly North of the pixel, and 90° when the satellite is directly East, 180° when the satellite is directly South and 270° when the satellite is directly West.
\(\theta_s\): The solar zenith angle, defined as 0° when the sun is directly overhead and 90° when the sun is at the horizon.
\(\phi_s\): The satellite azimuth angle, defined as 0° when the sun is directly North of the pixel, and 90° when the sun is directly East, 180° when the sun is directly South and 270° when the sun is directly West.
In addition, satellite sensors typically produce images for a series of channels, where each channel has a defined range of wavelengths over which it is sensitive. This is typically called the spectral response function (SRF) of a sensor.
At present, pycoxmunk does not take into account this SRF. Instead, the central wavelength for each channel is used. The spectral properties of water at various wavelengths are defined within pycoxmunk in the CM_Constants.py file, and for a given central wavelength the two closest wavelengths in the spectral library are interpolated to provide the wavelength-specific water properties for each channel. If a channel is selected that is outside the range of wavelengths supported by the library (~450nm - 3700nm) then the closest library wavelength will be used instead.
Finally, in order to compute the surface properties a source of near-surface wind information is needed. This could be from observations, numerical weather prediction models or a reanalysis such as the ECMWF/Copernicus ERA5 products. The 10 meter wind speeds in the ERA5 dataset have been found to produce best results.
Calculation of water properties
Before performing the reflectance calculations we first compute the properties of the water being observed. Currently, this is done using fixed values from the literature (described in Sayer, 2010) but in future we plan to implement support for ESA’s Ocean Color CCI data, which will enable more a dynamic computation of water properties.
Chlorophyll-A absorption and backscattering are computed via:
where \(coef_n\) are predefined coefficients extracted from the spectral properties library by wavelength and \(chl_{conc}\) is the Chlorophyll concentration, which is currently a fixed value of 0.18.
The total absorption and backscatter are then defined by:
Where \(base_{abs}\) and \(base_{bsc}\) are extracted from the spectral properties library by wavelength.
The white cap fraction is defined by:
Calculation of water reflectance
Now, the water body reflectance is calculated via:
The underlight contribution to total water reflectance is given by:
where \(t_u\) is the upward transmission, defined as 0.52, \(t_d\) is computed from the scene geometry (see Sayer, 2010 for more details) and \(r_u\) is the upward reflectance, defined as 0.48.
The sunglint reflectance is calculated by:
Where:
\(P{slo}\) is the wave slope distribution
\(r_{sf}\) is the Fresnel reflection coefficient and
\(\beta\) is the facet tilt defined by:
with \(\Theta\) defined via:
Finally, the total reflectance is calculated using:
Calculation of bidirectional reflectance
In addition to calculating the reflectance along the sun-surface-satellite path, pycoxmunk can also calculate the bidirectional reflectance terms. These comprise of the:
\(\rho_{0v}\): Solar beam to satellite view reflectances
\(\rho_{0d}\): Solar beam to diffuse reflectances
\(\rho_{dv}\): Diffuse to satellite view reflectances
\(\rho_{dd}\): Diffuse to diffuse reflectances
Calculation of these terms is enabled optionally by the user at runtime and is disabled by default. The BRDF terms require significant additional computation and hence extend both the processing time and memory requirements for pycoxmunk.
Internally, these terms are calculated by simulating the surface reflectance across a range of angles and summing the results to gain the direct and diffuse terms listed above. The surface reflectance is simulated using the terms described in the preceeding sections.
The angles used in the calculation are computed using the
Gauss-Legendre n-point quadrature formula and by default pycoxmunk uses
four terms for both \(\theta\) and \(\phi\). This can be altered within the pycoxmunk code by editing
the n_quad_theta and n_quad_phi terms in CM_Constants.py. In the future these will be user-configurable options.
Note, however, that increasing either of these values will significantly increase runtime and memory use.
The \(\rho_{0v}\) term is simply the standard Cox-Munk sea surface reflectance, so no additional calculations are performed for this term.
The \(\rho_{0d}\) and \(\rho_{dv}\) terms use the actual solar or viewing geometry (respectively) and then simulate diffuse radiation by substituting the viewing or solar geometry with the Gauss-Legendre terms described above. The following python code offers a simplified workflow for computing these values:
def get_rho_od_dv():
qx_qw_sincos = np.cos(gauss_leg_theta_abscissas) * np.sin(gauss_leg_theta_abscissas) * gauss_leg_theta_weights
# Loop over zeniths
for i in range(0, n_quad_theta):
cur_gl_zen = gauss_leg_theta_abscissas[i]
tmp_0d = 0
tmp_dv = 0
# Loop over azimuths
for j in range(0, n_quad_phi):
cur_gl_azi = gauss_leg_phi_abscissas[j]
cur_gl_phi_w = gauss_leg_phi_weights[j]
# Compute the reflectances
tmp_0d = tmp_0d + calc_cox_munk_refl(sol_zen, sol_az, cur_gl_zen, cur_gl_azi) * cur_gl_phi_w
tmp_dv = tmp_dv + calc_cox_munk_refl(cur_gl_zen, cur_gl_azi, sat_zen, sat_azi) * cur_gl_phi_w
rho_0d = rho_0d + tmp_0d * qx_qw_sincos[i]
rho_dv = rho_dv + tmp_dv * qx_qw_sincos[i]
rho_0d = rho_0d / np.pi
rho_dv = rho_dv / np.pi
Finally, the \(\rho_{dd}\) term is calculated in a similar way, but substitutes both the solar and viewing zenith angles as opposed to just one or the other:
def get_rho_dd():
qx_qw_sincos = np.cos(gauss_leg_theta_abscissas) * np.sin(gauss_leg_theta_abscissas) * gauss_leg_theta_weights
# Initialise direct_diffuse reflectance
rho_dd = 0
# Loop over solar zeniths
for i in range(0, n_quad_theta):
outer_tmp_dd = 0
cur_zen1 = gauss_leg_theta_abscissas[i]
# Loop over satellite zeniths
for j in range(0, n_quad_theta):
inner_tmp_dd = 0
cur_zen2 = gauss_leg_theta_abscissas[j]
# Loop over relative azimuths
for k in range(0, n_quad_phi):
cur_relazi = gauss_leg_phi_abscissas[k]
cur_phi_w = gauss_leg_phi_weights[k]
# Compute the reflectances
inner_tmp_dd = inner_tmp_dd + calc_cox_munk_refl(cur_zen1, .0, cur_zen2, cur_relazi) * cur_phi_w
outer_tmp_dd = outer_tmp_dd + inner_tmp_dd * qx_qw_sincos[j]
rho_dd = rho_dd + outer_tmp_dd * qx_qw_sincos[i]
rho_dd = rho_dd / np.pi