Atmospheric correction with RTTOV

[Andrew Prata]

I’m interested in testing whether or not RTTOV is suitable for atmospheric correction and BRDF retrieval.

I’m starting by assuming a Lambertian surface. According to Vermote et al. (1997) (i.e. the MODIS algorithm), for a Lambertian surface, we may write:

rho_TOA = rho_a + (rho_s * tau_dn * tau_up) / (1 – rho_s * S)

where,
rho_TOA is the top-of-atmosphere reflectance (I will get this from the satellite measurement).
rho_a is the intrinsic atmospheric reflectance (due to Rayleigh and aerosol scattering).
rho_s is the surface reflectance (what I want to retrieve).
tau_dn is the total (direct + diffuse) downward transmission from sun to surface (function of sun zenith angle).
tau_up is the total (direct + diffuse) upward transmission from surface to satellite (function of sat zenith angle).
S is the spherical albedo of the atmosphere.

I assume that if I set surfemisrefl(1,:,:) = 0 (i.e. rho_s = 0) then RTTOV should give me back rho_TOA = rho_a (i.e. the intrinsic atmospheric reflectance).

And according to the user guide, tausun_total_path2 is the “Transmittance for combined sun-surface-satellite path”. Assuming that the rttov_transmission structure contains total transmittances (i.e. direct + diffuse terms) then we have:
tausun_total_path2 = tau_dn * tau_up

Which leaves S to be determined. Can I output the spherical albedo from RTTOV somehow? Or would you recommend a different approach for atmospheric correction with RTTOV?

Thanks in advance.

Cheers
Andrew

Vermote, E. F., El Saleous, N., Justice, C. O., Kaufman, Y. J., Privette, J. L., Remer, L., Roger, J. C., and Tanré, D.: Atmospheric correction of visible to middle‐infrared EOS‐MODIS data over land surfaces: Background, operational algorithm and validation, J. Geophys. Res., 102, 17131–17141, https://doi.org/10.1029/97JD00201, 1997.

[James Hocking]

Hi Andrew,

If you set the surface reflectance to zero in RTTOV, then for visible/near-IR channels (with no thermal component) the only contribution to the TOA radiance is atmospheric Rayleigh scattering, as you require.

I’m not sure I understand when you talk about “direct+diffuse transmittance”. The transmittances output by RTTOV represent the level-to-space extinction by gas absorption and, where relevant, also by Rayleigh scattering exintinction (i.e., Rayleigh scattering out of the sensor line-of-sight). The tausun_total_path2 transmittances are the combined transmittances along the full sun-surface-satellite path. By contrast, the tausun_total_path1 transmittances are those along the surface-satellite path. These are transmittances for radiation passing through the layer(s) along the given paths.

Regarding the spherical albedo of the atmosphere, RTTOV cannot compute this directly. I wonder if it is possible to compute by running multiple simulations with differing solar geometry and integrating the results over the hemisphere?

I’m afraid I don’t know about any specific methods for retrieving BRDF.

Best wishes,
James

[Andrew Prata]

Hi James

Thank you for the helpful response.

I have been testing RTTOV output against Py6S (https://py6s.readthedocs.io/en/latest/introduction.html) for AVHRR channels 1 and 2 to see if I can get out the variables needed to implement the Vermote et al. approach and I’m finding some interesting similarities and differences.

By setting surfemisrefl(1,:,:) = 0 in RTTOV, I get reflectances of:
Ch1 = 0.040551
Ch2 = 0.015841

and for Py6S (which gives me atmospheric intrinsic reflectance as a variable by default), I get:
Ch1 = 0.0451549
Ch2 = 0.0176200

This seems to be an acceptable agreement between the two RTMs and in line with our discussion above.

However, I’m having trouble reconciling the differences between the transmittance output from RTTOV and the output I get from Py6S.

To clarify:
T(µs) is the total transmittance from the top of the atmosphere to the ground along the path of the incoming solar beam (I wrote this previously as tau_dn).

T(µv) is the total transmittance from the ground to the top of the atmosphere in the view direction of
the satellite (I wrote this previously as tau_up).

where, µs is the cosine of the solar zenith angle, µv is the cosine of the satellite view angle.

The total transmittance (for sun and view angles) can be broken down into direct and diffuse components:

T(µs) = e^-τ/µs + td(µs)
T(µv) = e^-τ/µv + td(µv)

where τ is the total optical thickness and td the diffuse transmittance.

I am currently assuming that, in RTTOV,

tausun_total_path1 = T(µv)
and
tausun_total_path2 = T(µs) * T(µv)

But perhaps I’m misinterpreting what you mean when you say “combined”.

More details on the above formulation can be found here: https://modis.gsfc.nasa.gov/data/atbd/atbd_mod08.pdf (see Eqs. 2 and 3, page 21)

If you could help me figure out what these variables are in RTTOV-terms that would be extremely useful.

Thanks again.
Andrew

[James Hocking]

Hi Andrew,

Your assumptions about tausun_total_path1/2 in RTTOV are correct. The only minor caveat is that we are dealing with transmittances for a finite spectral bandwidth and so in theory multiplying transmittances together is not strictly correct. Nevertheless it is a reasonable approximation in the visible. RTTOV computes the gas absorption optical depth regression for tausun_total_path2, and incorporates the Rayleigh extinction, and then uses this to obtain tausun_total_path1 by scaling the corresponding optical depths by the relevant path lengths (which is an approximation as it is equivalent to multiplying or dividing transmittances).

In order to be a fast model, RTTOV computes transmittances for one particular geometry (i.e., the sun-surface-satellite geometry for visible channels). The resulting transmittances are used for all calculations.

Rayleigh extinction is considered along the entire direct solar beam on the sun-surface-satellite path. For Rayleigh single-scattering towards the sensor, we consider the direct downward solar beam along the sun-surface path, and include the radiation scattered upwards towards the sensor (assumed to be along the surface-satellite path) and radiation scattered downwards along this same path that undergoes specular reflection at the surface and subsequently reaches the satellite. The Rayleigh extinction is of course included for this single-scattered radiation as well.

For Rayleigh multiple scattering via the DOM solver, the computed sun-surface-satellite transmittances are scaled by the relevant path lengths for each of the “ordinates” along which the radiation is resolved in the Discrete Ordinates Method. This is of course again an approximation.

I’m afraid I still don’t know what diffuse transmittance is. I’m not familiar with the idea of adding transmittances together or what that represents. Perhaps the above descriptions might help understand the differences between RTTOV and Py6S?

Best wishes,
James

[Andrew Prata]

Hi James

Apologies in advance (this is a long response).

Thank you for those clarifications – this additional information certainly helps my understanding.

To understand what the diffuse transmittance means in this context, consider a sensor at the ground viewing the sun through a clear atmosphere (analogous to the sun-surface path). Some fraction of solar radiation incident at the top of the atmosphere will be received by the ground-based sensor and can be broken down into four components that must sum to 1 (to account for all radiation incident at TOA).

The four possible fates of these photons are:
1) They are transmitted directly through the atmosphere (along the sun-surface path) and received by the ground-base sensor (i.e. undergo no scattering events). This fraction of the incident radiation is the ‘direct’ transmittance (t_dir).
2) They are scattered one or more times in the atmosphere and then are received (along the sun-surface path) by the ground-based sensor. This fraction of the incident radiation is the ‘diffuse transmittance’ (t_diff).
3) They are reflected back to space by the atmosphere and are not received at the surface (r).
4) What ever is not transmitted nor reflected must be absorbed (a).

These terms must sum to 1 and so we have:
t_dir + t_diff + r + a = 1.

This is where the idea of the ‘total’ transmittance comes from i.e. t = t_dir + t_diff
and so t + r + a = 1. Based on your description, it sounds like the diffuse and direct transmittance components are taken together in RTTOV. However, I’m still getting unexpected results using RTTOV (unexpected because I’m probably misunderstanding some of the settings I’m using in RTTOV).

For example, to simulate the TOA reflectance (along the surface-satellite path) for an homogeneous Lambertian surface with a defined surface reflectance of say 0.50, in pyrttov, I am doing:

# Create avhrr instance
avhrrRttov = pyrttov.Rttov()

# Set RTTOV options
rttov_installdir = ‘/home/565/ap3706/rttov13/’
avhrr_platform = ‘noaa_14_avhrr’
avhrrRttov.FileCoef = ‘{}/{}’.format(rttov_installdir, “rtcoef_rttov13/rttov13pred54L/rtcoef_” + avhrr_platform + “_o3co2.dat”)
avhrrRttov.Options.AddInterp = True
avhrrRttov.Options.AddSolar = True
avhrrRttov.Options.VerboseWrapper = True
avhrrRttov.Options.StoreTrans = True
avhrrRttov.Options.StoreRad = True
avhrrRttov.Options.StoreRad2 = True
avhrrRttov.Options.StoreEmisTerms = True
avhrrRttov.Options.OzoneData = True
avhrrRttov.Options.UseQ2m = False

# Load instrument
avhrrRttov.loadInst()

# Read in atmospheric profiles
P, T, Q, O3 = pd.read_csv(path_to_data)
nprofiles = 1
nlevels = len(P)
myProfiles = pyrttov.Profiles(nprofiles, nlevels)
myProfiles.GasUnits = 1
myProfiles.P = [P]
myProfiles.T = [T]
myProfiles.Q = [Q]
myProfiles.O3 = [O3]

myProfiles.DateTimes = [[1999, 6, 1, 6, 0, 0]]
myProfiles.Angles = [[7.694666971847099, 258.4655527580682, 78.22247653925383, 306.85557335774774]]
myProfiles.SurfGeom = [[-33.76589584350586, 142.7560272216797, 0.0]]
myProfiles.SurfType = [[0, 0]]
myProfiles.S2m = [[1017.3526, 288.77625, 0.0, -1.9953531, 1.4665465, 0.0]]
myProfiles.Skin = [[290.4469, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]

# Assign myProfiles to avhrr instance
avhrrRttov.Profiles = myProfiles

# Set surface parameters
nchan_avhrr = 5
surfemisrefl_avhrr = np.zeros((5, nprofiles, nchan_avhrr), dtype=np.float64)
avhrrRttov.SurfEmisRefl = surfemisrefl_avhrr
avhrrRttov.SurfEmisRefl[0,:,:] = 1.
avhrrRttov.SurfEmisRefl[1,:,:] = 0.50

# Run direct model
avhrrRttov.runDirect()

# Results
avhrrRttov.Refl
array([[1.01585354, 1.09382368, 1.30480583, 0. , 0. ]])

What is strange to me is that I’m getting reflectances greater than 1. If the surface reflectance is 0.5 (for a Lambertian reflector) then I would expected the TOA reflectance to be < 0.5.

I also tried to force the DOM solver for clear sky by adding the following options to the above code:
avhrrRttov.FileSccld = ‘{}/{}’.format(rttov_installdir, “rtcoef_rttov13/cldaer_visir/sccldcoef_” + avhrr_platform + “.dat”)
avhrrRttov.Options.DomRayleigh = True
avhrrRttov.Options.AddClouds = True

This gives me:
avhrrRttov.Refl
array([[0.78584153, 0.74329933, 0.98264457, 0. , 0. ]])

These results seem a bit more reasonable but still I would expect: TOA reflectance < surface reflectance for the solar channels.

Do you have any ideas on how to explain these results?

Also if you’d like to try and reproduce my results I can send you the atmospheric profiles I’m using (they’re derived from ERA5 data).

Again – thank you for spending time on this problem.

Cheers
Andrew

[James Hocking]

Hi Andrew,

Thanks for the explanation, that is useful.

One thing to point out immediately: the surface “reflectance” you are specifying as input to RTTOV is a BRDF (bi-directional reflectance function) value for the given sun-surface-satellite geometry while the output reflectance from RTTOV is a BRF (bi-directional reflectance factor) consistent with what is measured by satellites.

For a Lambertian surface with reflectance (albedo) 0.5, the BRDF should be set to 0.5/pi (because the incoming radiation is reflected equally in all directions with this reflectance).

In the absence of any atmosphere at all, the output BRF is equal to pi*BRDF.

And, just for clarity, physical BRDFs and BRFs may be greater than 1.

To run the DOM solver you would additionally need to set this as the solver used for VIS channels:
avhrrRttov.Options.VisScattModel = 1

Let me know whether the above makes sense of your results.

Best wishes,
James

[Andrew Prata]

Hi James

Thank you for pointing out the BRDF issue. Using 0.5/pi I get much better results. I’m nearly there; however, there is still something that seems inconsistent.

As a self-consistency check, I should be able to get back the surface reflectance that I put into RTTOV if I assume the surface is homogenous Lambertian.

Starting with:
rho_TOA = rho_a + (rho_s * T_tot) / (1 – rho_s * S),
where T_tot = T(µs) * T(µv).

I can solve for the surface reflectance:
rho_s = C / (1. + C * S)
where, C = (rho_TOA – rho_a) / T_tot

Now, I can get the values of rho_a and S from RTTOV as follows:

For rho_a, the intrinsic atmospheric reflectance, I set avhrrRttov.SurfEmisRefl[1,:,:] = 0., re-run RTTOV and then extract the result of:
rho_a = avhrrRttov.Refl[0][0]
rho_a = 0.04463625 # For AVHRR Channel 1

For S, I use the parameterisation given in the 6SV User Guide:

from scipy.special import expn

def compute_spherical_albedo(tau_R):
S = (1./(4 + 3 * tau_R)) * (3 * tau_R – 4 * expn(3, tau_R) + 6 * expn(4, tau_R))
return S

where tau_R is the Rayleigh optical depth. To get tau_R from RTTOV, I set avhrrRttov.Options.DoOpdepCalc = False, re-run RTTOV, and extract the result of:
tau_R = -1. * np.log(avhrrRttov.TauSunTotalPath1) * np.cos(zenangle * np.pi/180.)
tau_R = 0.0529807 # For AVHRR Channel 1

Then the spherical albedo from RTTOV is:
S = compute_spherical_albedo(tau_R)
S = 0.047595370 # For AVHRR Channel 1

Finally, I re-run RTTOV with the following settings:
avhrrRttov.Options.DoOpdepCalc = True
avhrrRttov.SurfEmisRefl[1,:,:] = 0.5/np.pi

Now I have:
rho_TOA = avhrrRttov.Refl[0][0]
T_tot = avhrrRttov.TauSunTotalPath2[0][0]

So I can evaluate:
C = (rho_TOA – rho_a) / T_tot
which gives me back rho_s
rho_s = C / (1. + C * S)

In theory, this should be rho_s = 0.5.
However, when I do the above, I get rho_s = 0.58654.

This is close but is different enough to make me wonder if I’m missing something. I think it must be to do with the transmittance, T_tot, as defined above, because the other values (i.e. tau_R, S and rho_a) are comparable to the values I get by other means (e.g. from Py6S output).

Any further ideas on my above settings/assumptions?

Thanks again.
Andrew

[Author]

Posts