# 6S atmospheric correction results wrong

I'm using 6s to correct sentinel data, but there is a little deviation that tells me that I'm doing something wrong.

This is the result for example of one scenario:

``````*******************************************************************************
*                        atmospheric correction result                        *
*                        -----------------------------                        *
*       input apparent reflectance            :    0.100                      *
*       measured radiance [w/m2/sr/mic]       :   37.699                      *
*       atmospherically corrected reflectance                                 *
*       Lambertian case :      0.02879                                        *
*       BRDF       case :      0.02879                                        *
*       coefficients xa xb xc                 :  0.00346  0.10170  0.14026    *
*******************************************************************************
``````

The problem is that if I calculate the atmospherically corrected reflectance (acr) using the formulas provided with a measured radiance of 37.699, I get a value of 0.02862 which does not match with the presented result = 0.02879. I know it's just a 0,6% error, but it's math... there should be no error.

I guess I'm doing something wrong and need guidance in terms of how to use these values to correct sentinel data, as it seems I am missing a step maybe, which i don't find in 6s' documentation. I need to make the calculation pixel by pixel, so running the algorithm n times (one for each pixel) is not an option.

Short answer: It's a numerical precision issue: the `xa`, `xb` and `xc` coefficients are not displayed in the 6S output file with enough precision to produce the correct answer when using the formula. It is possible to modify 6S to print the coefficients with a higher precision, which solves the problem.

Long answer: I was intrigued by this question, as I've used 6S for a long time (and written a Python interface to it) and never noticed this before...strangely, I'd never thought to check!

I thought originally that it might be an issue with your 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 *
*******************************************************************************
``````

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 `xa` coefficient 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               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               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   *
*******************************************************************************
``````

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 you 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 difference_between_formula_and_model(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()
• another thing to consider is the data precision type. Skimming thru `main.f` all the types are `real` / single-precision / `float32`. Python natively uses double precision, so you will always get a better answer with higher precision. May 15, 2018 at 3:13