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    *
*       y=xa*(measured radiance)-xb;  acr=y/(1.+xc*y)                         *
*******************************************************************************

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.

up vote 8 down vote accepted
+50

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 *
* 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 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           ' 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 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 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.

  • Wow! thanks @robintw for your time to answer this so thoroughly, i will be using Py6S (thanks also for that too). I'll try to recompile 6S following your steps to solve this! – luckyshonway May 11 at 21:08
  • 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. – Mike T May 15 at 3:13

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.