Robin's Blog

Mismatch between 6S atmospheric correction results & those from coefficients

A while back a friend on Twitter pointed me towards a question on the GIS StackExchange site about the 6S model, asking if "that was the thing you wrote". I didn’t write the 6S model (Eric Vermote and colleagues did that), but I did write a fairly well-used Python interface to the 6S model, so I know a fair amount about it.

The question was about atmospherically correcting radiance values using 6S. When you configure the atmospheric correction mode in 6S you give it a radiance value measured at the sensor and it outputs an atmospherically-corrected radiance value. Simple. However, it also outputs three coefficients: x_a, x_b and x_c which can be used to atmospherically-correct other at-sensor radiance values. These coefficients are used in the following formulae, given in the 6S output:

y=xa*(measured radiance)-xb
acr=y/(1.+xc*y) 

where acr is the atmospherically-corrected radiance.

The person asking the question had found that when he used the formula to correct the same radiance that he had corrected using 6S itself, he got a different answer. In his case, the result from 6S itself was 0.02862, but when he ran his at-sensor radiance through the formula he got a different answer: 0.02879, a difference of 0.6%.

I was intrigued by this question, as I’ve used 6S for a long time and never noticed this before…strangely, I’d never thought to check! The rest of this post is basically a copy of my answer on the StackExchange site, but with a few bits of extra explanation.

I thought originally that it might be an issue with the parameterisation of 6S – but I tried a few different parameterisations myself and came up with the same issue – I was getting a slightly different atmospherically-corrected reflectance when putting the coefficients through the formula, compared to the reflectance that was output by the 6S model directly.

The 6S manual is very detailed, but somehow never seems to answer the questions that I have – for example, it doesn’t explain anywhere how the three coefficients are calculated. It does, however, have an example output file which includes the atmospheric correction results (see the final page of Part 1 of the manual). This includes the following outputs:

*******************************************************************************
* atmospheric correction result *
* ----------------------------- *
* input apparent reflectance : 0.100 *
* measured radiance [w/m2/sr/mic] : 38.529 *
* atmospherically corrected reflectance *
* Lambertian case : 0.22180 *
* BRDF case : 0.22180 *
* coefficients xa xb xc : 0.00685 0.03885 0.06835 *
* y=xa*(measured radiance)-xb; acr=y/(1.+xc*y) *
*******************************************************************************

If you work through the calculation using the formula given you find that the result of the calculation doesn’t match the 6S output. Let me say that again: in the example provided by the 6S authors, the model output and formula don’t match! I couldn’t quite believe this…

So, I wondered if the formula was some sort of simple curve fitting to a few outputs from 6S, and would therefore be expected to have a small error compared to the actual model outputs. As mentioned earlier, the manual explains a lot of things in a huge amount of detail, but is completely silent on the calculation of these coefficients. Luckily the 6S source code is available to download. Less conveniently, the source code is in written in Fortran 77!

I am by no means an expert in Fortran 77 (in fact, I’ve never written any Fortran code in real-life), but I’ve had a dig in to the code to try and find out how the coefficients are calculated.

If you want to follow along, the code to calculate the coefficients starts at line 3382 of main.f. The actual coefficients are set in lines 3393-3397:

 xa=pi*sb/xmus/seb/tgasm/sutott/sdtott
 xap=1./tgasm/sutott/sdtott
 xb=ainr(1,1)/sutott/sdtott/tgasm
 xb=ainr(1,1)/sutott/sdtott/tgasm
 xc=sast

(strangely xb is set twice, to the same value, and another coefficient xap is set, which never seems to be used – I have no idea why!).

It’s fairly obvious from this code that there is no complicated curve fitting algorithm used – the coefficients are simply algebraic manipulations of other variables used in the model. For example, xc is set to the value of the variable sast, which, through a bit of detective work, turns out to be the total spherical albedo (see line 3354). You can check this in the 6S output: the value of xc is always the same as the total spherical albedo which is shown a few lines further up in the output file. Similarly xb is calculated based on various variables including tgasm, which is the total global gas transmittance and sdtott, which is the total downward scattering transmittance, and so on. (These variables are difficult to decode, because Fortran 77 has a limit of six characters for variable names, so they aren’t very descriptive!).

I was stumped at this point, until I thought about numerical precision. I realised that the xacoefficient has a number of zeros after the decimal point, and wondered if there might not be enough significant figures to produce an accurate output when using the formula. It turned out this was the case, but I’ll go through how I altered the 6S code to test this.

Line 3439 of main.f is responsible for writing the coefficients to the file. It consists of:

write(iwr, 944)xa,xb,xc

This tells Fortran to write the output to the file/output stream iwr using the format code specified at line 944, and write the three variables xa, xb and xc. Looking at line 944 (that is, the line given a Fortran line number of 944, which is actually line 3772 in the file…just to keep you on your toes!) we see:

  944 format(1h*,6x,40h coefficients xa xb xc                 :, 
     s           1x, 3(f8.5,1x),t79,1h*,/,1h*,6x,
     s           ' y=xa*(measured radiance)-xb;  acr=y/(1.+xc*y)',
     s               t79,1h*,/,79(1h*))

This rather complicated line explains how to format the output. The key bit is 3(f8.5,1x) which tells Fortran to write a floating point number (f) with a maximum width of 8 characters, and 5 decimal places (8.5) followed by a space (1x), and to repeat that three times (the 3(...)). We can alter this to print out more decimal places – for example, I changed it to 3(f10.8,1x), which gives us 8 decimal places. If we do this, then we find that the output runs into the *‘s that are at the end of each line, so we need to alter a bit of the rest of the line to reduce the number of spaces after the text coefficients xa xb xc. The final, working line looks like this:

  944 format(1h*,6x,35h coefficients xa xb xc            :, 
 s           1x, 3(f10.8,1x),t79,1h*,/,1h*,6x,
 s           ' y=xa*(measured radiance)-xb;  acr=y/(1.+xc*y)',
 s               t79,1h*,/,79(1h*))

If you alter this line in main.f and recompile 6S, you will see that your output looks like this:

*******************************************************************************
*                        atmospheric correction result                        *
*                        -----------------------------                        *
*       input apparent reflectance            :    0.485                      *
*       measured radiance [w/m2/sr/mic]       :  240.000                      *
*       atmospherically corrected reflectance                                 *
*       Lambertian case :      0.45439                                        *
*       BRDF       case :      0.45439                                        *
*       coefficients xa xb xc            : 0.00297362 0.20291930 0.24282509   *
*       y=xa*(measured radiance)-xb;  acr=y/(1.+xc*y)                         *
*******************************************************************************

If you then apply the formula you will find that the output of the formula, and the output of the model match – at least, to the number of decimal places of the model output.

In my tests of this, I got the following for the original 6S code:

• Model: 0.4543900000
• Formula: 0.4537049078
• Perc Diff: 0.1507718536%

(the percentage difference I was getting was smaller than the questioner found – but that will just depend on the parameterisation used)

and this for my altered 6S code:

• Model: 0.4543900000
• Formula: 0.4543942552
• Perc Diff: -0.0009364659%

A lot better!

For reference, to investigate this I used Py6S, the Python interface to the 6S model that I wrote. I used the following functions to automatically calculate the results using the formula from a Py6S SixS object, and to calculate the percentage difference automatically:

def calc_acr(radiance, xa, xb, xc):
    y = xa * radiance - xb
    acr = y/(1.0 + xc * y)

    return acr

def calc_acr_from_obj(radiance, s):
    return calc_acr(radiance, s.outputs.coef_xa, s.outputs.coef_xb, s.outputs.coef_xc)

def difference_between_formula_and_model(s):
    formula = calc_acr_from_obj(s.outputs.measured_radiance, s)
    model = s.outputs.atmos_corrected_reflectance_lambertian

    diff = model - formula

    perc_diff = (diff / model) * 100

    print("Model: %.10f" % model)
    print("Formula: %.10f" % formula)
    print("Perc Diff: %.10f%%" % perc_diff)

and my example errors above came from running Py6S using the following parameterisation:

s = SixS()
s.altitudes.set_sensor_satellite_level()
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromRadiance(240)
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)
s.run()

Just as a slight addendum, if you’re atmospherically-correcting Sentinel-2 data with 6S then you might want to consider using ARCSI – an atmospheric correction tool that uses Py6S internally, but does a lot of the hard work for you. The best way to learn ARCSI is with their tutorial document.

If you found this post useful, please consider buying me a coffee.
This post originally appeared on Robin's Blog.

Tagged with:

Categorised as: Academic, My Software, Py6S, Remote Sensing