Bug in Rayleigh transmittance if sensor is set on custom altitude?

[Juha S]

First of all, thank you for publishing openly the Py6S package. I really appreciate such a free, handy tool!

 

What I am trying to do, is to calculate atmospheric transmittance between a ground target and a remote sensing drone (at e.g. 100m altitude). While trying to calculate that I found some really strange behavior (a bug?) in Rayleigh upward transmittance result. This seems to occur when sensor altitudes is set with set_sensor_custom_altitude() instead of set_sensor_sea_level() or set_sensor_satellite_level().

 

To reproduce the error, I ran the 6S simulation using default Py6S parameters on sea level, on satellite level, and on 4 custom altitudes in between. This is how the result: “transmittance_rayleigh_scattering.upward” comes out:

 

 

 

As expected the “Sea level” (with set_sensor_sea_level()) has transmittance of 1 while the “Satellite” (with set_sensor_satellite_level()) shows about what you should expect. However, the custom altitude simulations (with set_sensor_custom_altitude()) are not in between these two! The "99km" transmittance is basically the same as the satellite, but the rest of them start to develop to wrong direction! The error is clear when noticing that the ”0.1km” curve is nowhere close to the “Sea level” curve.

 

I believe I am setting the altitude correctly and this is a bug in the Py6S/6S as this strange behavior happens only with “transmittance_rayleigh_scattering.upward”. When I run similar plots on transmittance_aerosol_scattering.upward and transmittance_global_gas.upward, they act as expected:

 

 

Here is my minimum python code reproducing the problem and plotting the Rayleigh figure above:

 

#%% Simulate atmosphere under different drone altitudes
import numpy as np
import Py6S

Wavelengths_um = np.arange(350, 1001, 50)/1000

# Altitude set to SEA LEVEL
s = Py6S.SixS()
s.altitudes = Py6S.Altitudes()
s.altitudes.set_sensor_sea_level()
s.altitudes.set_target_sea_level()
_, results_SEA = Py6S.SixSHelpers.Wavelengths.run_wavelengths(s, Wavelengths_um)

# Altitude set to 0.1km
s = Py6S.SixS()
s.altitudes = Py6S.Altitudes()
s.altitudes.set_sensor_custom_altitude(0.100)
s.altitudes.set_target_sea_level()
_, results_100m = Py6S.SixSHelpers.Wavelengths.run_wavelengths(s, Wavelengths_um)

# Altitude set to 1km
s = Py6S.SixS()
s.altitudes = Py6S.Altitudes()
s.altitudes.set_sensor_custom_altitude(1.0)
s.altitudes.set_target_sea_level()
_, results_1km = Py6S.SixSHelpers.Wavelengths.run_wavelengths(s, Wavelengths_um)
 
# Altitude set to 10km
s = Py6S.SixS()
s.altitudes = Py6S.Altitudes()
s.altitudes.set_sensor_custom_altitude(10.0)
s.altitudes.set_target_sea_level()
_, results_10km = Py6S.SixSHelpers.Wavelengths.run_wavelengths(s, Wavelengths_um)
 
# Altitude set to 100km
s = Py6S.SixS()
s.altitudes = Py6S.Altitudes()
s.altitudes.set_sensor_custom_altitude(99.0)
s.altitudes.set_target_sea_level()
_, results_99km = Py6S.SixSHelpers.Wavelengths.run_wavelengths(s, Wavelengths_um)
 
# Altitude set to SATELLITE
s = Py6S.SixS()
s.altitudes = Py6S.Altitudes()
s.altitudes.set_sensor_satellite_level()
s.altitudes.set_target_sea_level()
_, results_SAT = Py6S.SixSHelpers.Wavelengths.run_wavelengths(s, Wavelengths_um)
 
#%% Plot transmittance spectra
def transmittance_spectrum(results):
    return [result.transmittance_rayleigh_scattering.upward for result in results]
 
import matplotlib.pyplot as plt
 
plt.figure()
 
plt.plot(Wavelengths_um, transmittance_spectrum(results_SEA), 'ko:', label='Sea level')
plt.plot(Wavelengths_um, transmittance_spectrum(results_100m), 'kx-', label='0.1km')
plt.plot(Wavelengths_um, transmittance_spectrum(results_1km), 'bx-', label='1km')
plt.plot(Wavelengths_um, transmittance_spectrum(results_10km), 'gx-', label='10km')
plt.plot(Wavelengths_um, transmittance_spectrum(results_99km), 'rx-', label='99km')
plt.plot(Wavelengths_um, transmittance_spectrum(results_SAT), 'ro:', label='Satellite')
 
plt.ylim([0.45, 1.05])
plt.grid()
plt.legend()
plt.xlabel('Wavelength (um)')
plt.ylabel('Rayleigh transmittance')

 

My Py6S version is "1.7.2_py_1" that was installed this week from conda-forge.

 

Am I doing something wrong or have I discovered a bug in simulation/in output of results?

 

Best regards

-Juha Suomalainen

[Robin Wilson]

Hi Juha,

Thanks for sending such detailed plots and code - that's really helped me to look into this.

To be honest, I'm a bit stuck as to what is going on. I've done some fairly extensive testing of Py6S this morning and I'm sure this isn't a bug in Py6S. I've checked that the inputs are being specified correctly when the sensor altitude is being set, and have checked that the outputs from the 6S output file (which can be viewed with s.outputs.fulltext) are being parsed correctly. I've also used the online version of 6S provided by the 6S authors (http://6s.ltdri.org/pages/run6SV.html) and that gives the same outputs as Py6S, and I can generate the same input file that Py6S generates by answering the questions in the web interface.

So, if it's not a bug with Py6S then it must be the results coming from 6S, and either a) the 6S outputs are wrong, or b) this is expected behaviour.

I would be surprised if there was a bug like this in 6S itself, as 6S has been validated against a number of other radiative transfer models, and has always found good agreement (admittedly, I don't know how much these comparisons involved things like custom altitudes), but I also agree with you that these results look very strange, and I can't really think of a situation in which they're right.

As you probably know, I didn't write the 6S code itself, I just wrote the Python interface. The 6S code is available, and I've spent a while looking through it, but it's rather hard to interpret (it's written in Fortran 77) and it's hard to see exactly how the flow of data through the various calculations works. There seems to be a scaling factor which scales something related to Rayleigh transmittance (I've yet to establish exactly what the variable name means - they're all a maximum of 6 characters long, which makes it a bit difficult!) and is large (nearly 1) when the altitude is near the satellite altitude, and small (near 0.01) when the altitude is near the ground - and if anything this seems like it might be the wrong way around, but I'm not 100% sure. Unfortunately I'm really not experienced enough in the details of the 6S code itself to be able to understand it all.

I've also tested 6SV2.1 (a later version of 6S, which isn't very widely used and isn't yet supported by Py6S) and it gives the same results as 6SV1.1.

All I can really suggest is trying to contact the original 6S authors and see if they have any insight in to what is going on here. I would contact them myself, but I've never had a response from any of them when I've emailed them about issues with 6S - so unfortunately I suspect you may not get a response. However, it's always worth a try and a carefully chosen subject line may attract their attention.

I'm really sorry I can't help more - if you want any more advice (or guidance on where to look in the 6S source code) then I will try and help as much as I can,

Best regards,

Robin

--
You received this message because you are subscribed to the Google Groups "Py6S" group.
To unsubscribe from this group and stop receiving emails from it, send an email to py6s+uns...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/0fe91f9a-8d25-43d5-bcb4-5dbad53f8e0c%40googlegroups.com.

[Juha S]

Thank you for the reply and detailed analysis of the root problem. I ran some simulations more and I think the solution to this is just to ignore the upward and total values on Rayleigh transmittance row.

Proof: here is the full text of simulation with sensor on 0.1km altitude:

*                             downward        upward          total           *

...

*      rayl.  sca. trans. :     0.93267        0.86730        0.80890         *

*      aeros. sca.   "    :     0.89678        0.99507        0.89236         *

*      total  sca.   "    :     0.83040        0.99442        0.82577         *

The total scattering row should be result of multiplication of rayleigh and aerosol rows. We can see here that this is not the case on upward and total columns. The fact that the total scattering row is different (and actually has sensible values in different altitudes!), would indicate that internally the 6S is using some other Rayleigh value than what is shown in output file. Most likely the 6S has just a bug where it retrieves wrong value for "upward Rayleigh transmittance" when writing the output. 

This solves the issue for me, as I actually don't even need the Rayleigh transmittance itself for anything. I just need to have the total upward transmittance which I believe is:

transmittance_total  = result.transmittance_total_scattering.upward *result.transmittance_global_gas.upward

I just got confused as the Rayleigh seemed NOT to be included the in the "transmittance_total_scattering" and thought I had to add it separately.

If you want, you could add a 2 line fix to the Py6S as it is still possible to determine the correct Rayleigh values using the total_scattering values:

result.transmittance_rayleigh_scattering.upward = result.transmittance_total_scattering.upward / result.transmittance_aerosol_scattering.upward
result.transmittance_rayleigh_scattering.total= result.transmittance_total_scattering.total/ result.transmittance_aerosol_scattering.total

Thanks!
-Juha

[Robin Wilson]

Hi,

I agree that the total scattering transmittance should be the result of multiplying the Rayleigh and aerosol scattering transmittances. Interestingly, if you look at the example 6S output files provided by the 6S authors, this holds for all of the examples except when the example is parameterised to have the sensor at a non-ground and non-satellite altitude (specifically, it fails for Example_Out_1.txt but works for numbers 2, 3 and 4).

You're right that if the upward value is wrong then the total value will also be wrong: in the 6S source code, the total value is just the downward multiplied by the upward.

Looking into how 6S calculates these: the total value seems to be calculated completely differently from the Rayleigh and aerosol values: it isn't just a simple multiplication of the two, it's far more complicated than that. In fact, I'm struggling to see how the code relates these at all!

I'm wary of putting a correction in Py6S for this until we get it confirmed as a definite 6S bug - as one of the 'selling points' of Py6S is that all the numbers it provides are exactly the same as you'd get from running 6S manually. However, I'm wondering if I should put in some sort of warning when you access the upward Rayleigh transmittance, which directs to a part of the docs. I'd need to write that docs section first though!

I think it would be very useful to have some input on this from the people who wrote 6S. Are you happy to email them about this? They might take more notice of an email from someone who hasn't pestered them before! Feel free to cc me in to the email if you want, and include your graph and your observation that total should be equal to rayleigh * aerosol. They probably won't be experienced with Py6S code, so you might need to send some example 6S input files instead (you can get these from Py6S by calling s.write_input_file(), which returns a path to a temporary file containing the input file sent to 6S).

Thanks again for your engagement on this,

Robin

--
You received this message because you are subscribed to the Google Groups "Py6S" group.
To unsubscribe from this group and stop receiving emails from it, send an email to py6s+uns...@googlegroups.com.

To view this discussion on the web, visit https://groups.google.com/d/msgid/py6s/846c1817-74b3-4499-9d2a-f80b4fa00b81%40googlegroups.com.