Earthdata Forum

Dr. :
I computing rayleigh radiance using a code that is similar to "rayrad.pro" https://oceancolor.gsfc.nasa.gov (/cgi/idllibrary.cgi?dir=atmocorr). However, the result is different with seadas product "Lr". the scatter in the following attachment file shows difference at eight bands. I find，after analysis in ENVI，the difference increase from center to edge of the scan line. the figure and their description and my code are as follows. Does anyone know the reason?

# -*- coding: utf-8 -*-
import numpy as np
from netCDF4 import Dataset
from scipy import interpolate
import general
def rayleigh(rayleigh_lut_path=None, sza=None, vza=None, saa=None, vaa=None, windspeed=None, pressure=None, F0=None):
"""
Args:
rayleigh_lut_path ():
sza (): solar zenith
vza (): sensor zenith
vaa (): sensor azimuth
saa (): solar azimuth
windspeed (): wind speed
pressure (): atmospheric pressure
Returns:
Lr_i
"""
reaa=vaa-180-saa
reaa[reaa<-180]=reaa[reaa<-180]+360
reaa[reaa > 180] = reaa[reaa > 180] - 360
rayleigh_lut = general.get_filelist(rayleigh_lut_path, 'rayleigh', 'iqu.hdf')
windspeed[windspeed > 30] = 30
sza[sza > 88] = 88
mu0 = np.cos(sza / 180 * np.pi)
mu = np.cos(vza / 180 * np.pi)
airmass = 1 / mu0 + 1 / mu
Taur = np.zeros(shape=(rayleigh_lut.__len__()))
Lr_i = np.zeros(shape=(sza.shape[0], sza.shape[1], rayleigh_lut.__len__()))

for i in range(rayleigh_lut.__len__()):
rayDtset = Dataset(rayleigh_lut[i])
taur = rayDtset.variables['taur'][:]
depol = rayDtset.variables['depol'][:]
senz = rayDtset.variables['senz'][:]
solz = rayDtset.variables['solz'][:]
wind = rayDtset.variables['wind'][:]
# sigma = rayDtset.variables['sigma'][:]
i_ray = rayDtset.variables['i_ray'][:]
# q_ray = rayDtset.variables['q_ray'][:]
# u_ray = rayDtset.variables['u_ray'][:]
Taur[i] = taur

Norder0 = np.zeros_like(sza)
Norder1 = np.ones_like(vza)
Norder2 = Norder0 + 2

ray_i0 = interpolate.interpn(
(wind.reshape(-1), solz.reshape(-1), np.arange(3).reshape(-1), senz.reshape(-1)), i_ray,
np.stack([windspeed, sza, Norder0, vza], axis=2))
ray_i1 = interpolate.interpn(
(wind.reshape(-1), solz.reshape(-1), np.arange(3).reshape(-1), senz.reshape(-1)), i_ray,
np.stack([windspeed, sza, Norder1, vza], axis=2))
ray_i2 = interpolate.interpn(
(wind.reshape(-1), solz.reshape(-1), np.arange(3).reshape(-1), senz.reshape(-1)), i_ray,
np.stack([windspeed, sza, Norder2, vza], axis=2))

ray_i = ray_i0 + ray_i1 * np.cos(reaa / 180 * np.pi) + ray_i2 * np.cos(2 * reaa / 180 * np.pi)

# from Wang menghua ;is from Seadas
p0 = 1013.25 # hpa
x = (-(0.6543 - 1.608 * taur) + (0.8192 - 1.2541 * taur) * np.log(airmass)) * taur * airmass
fac = ((1.0 - np.exp(-x * pressure / p0)) / (1.0 - np.exp(-x)))
Lr_i[:, :, i] = ray_i * fac * F0[i]

return Lr_i

Hi, the Rayleigh table coefficients i_ray have the dimensions of wave mean square slope *sigma* and not *windspeed*. Before the interpolation convert windspeed to sigma using sigma = sqrt(0.00534*WS).