robintw / Py6S

A Python interface to the 6S Radiative Transfer Model
GNU Lesser General Public License v3.0
186 stars 105 forks source link

set_target_pressure has no effect #72

Closed simonrp84 closed 3 years ago

simonrp84 commented 3 years ago

Hello, I'm trying to use the SixS.altitudes.set_target_pressure command but it doesn't seem to have any effect on the simulation results. For example:

# Initialise
s = SixS()

# Sensor is TOA
s.altitudes.set_sensor_satellite_level()

# Set atmospheric profile with WV and Ozone
s.atmos_profile = AtmosProfile.UserWaterAndOzone(3.6, 0.9)

# Set aerosol profile
s.aero_profile = AeroProfile.PredefinedType(AeroProfile.Maritime)

# Set wavelength range
s.wavelength = Wavelength(0.5)

for pres in np.arange(400, 1000, 50):
    s.altitudes.set_target_pressure(pres)
    s.run()
    print(pres,
          s.outputs.atmospheric_intrinsic_reflectance,
          s.outputs.total_gaseous_transmittance)

Returns:

400 0.094 0.942
450 0.094 0.942
500 0.094 0.942
550 0.094 0.942
600 0.094 0.942
650 0.094 0.942
700 0.094 0.942
750 0.094 0.942
800 0.094 0.942
850 0.094 0.942
900 0.094 0.942
950 0.094 0.942

i.e: The reflectance and transmittance do not appear to change.

If I instead iterate over altitude using:

for alt in np.arange(10, 0, -1):
    s.altitudes.set_target_custom_altitude(alt)

Then as expected I do get changes in the above values:

10 0.054 0.946
9 0.056 0.945
8 0.059 0.945
7 0.062 0.945
6 0.065 0.944
5 0.069 0.944
4 0.073 0.943
3 0.078 0.943
2 0.083 0.943
1 0.088 0.942

Do you have any ideas why set_target_pressure might not be working?

robintw commented 3 years ago

Thanks for reporting this issue.

I've looked into it, and it seems that setting a target pressure may not actually be supported by the underlying 6S model, which is very strange.

If you go to the 6S website you have an option to run 6S interactively. You can go through a series of questions to set the various parameters, and then you can see the 6S input file it generated, plus the output of running that through 6S.

If you go to the Target and Sensor Altitude page, there is an option for setting pressure in millibars. You can set this, and see that the value is shown in the resulting input file as a positive number in millibars. However, if you work through the web interface setting two different values for pressure, then you'll find that the output from 6S is the same. In fact, if you look earlier in the output file you'll see that for two different pressures (I used 500mb and 1000mb) you get the same output in the part of the output file that shows your inputs:

*                       target elevation description                          *
*                       ----------------------------                          *
*           ground pressure  [mb] 1018.00                                     *
*           ground altitude  [km] 0.000                                       *

I've looked into the 6S code to try and find where it processes the specified pressure value. The same input value is used for both setting a pressure (as a positive value, in mb) and setting an altitude (as a negative value, in km). This parameter is called xps in the 6S Fortran code. In the 6S code that processes this, it says:

if (xps.ge.0.) then
        xps=0.
        uwus=1.424
        uo3us=0.344

That is, if xps > 0 then it just sets xps to 0 to represent ground level. There is also a comment above this line of code that says xps >=0. means the target is at the sea level.

So, it seems that the 6S web interface and 6S manual both say that you can set the target altitude as a pressure - but the actual 6S code doesn't support this.

I'm sorry that Py6S provided an option that doesn't work - but I think you can probably see how that happened. I'll make sure this option is removed in the next release.

Hope this helps

simonrp84 commented 3 years ago

Thanks for looking into this. I also took a look into the 6S source today and noticed the same thing (hadn't checked git so didn't see your msg!)

Rather than removing this from Py6S, it'd be better to fix 6S itself as I imagine this is a common use case. Do you host the 6S executable yourself? I see it comes from the sixs package but don't understand enough about conda to know if that's yours or something else...

robintw commented 3 years ago

I'm afraid I'm not the original author of 6S, and I would rather not modify the 6S source code in any way, for various reasons: the original authors aren't responsive to emails, there is no license saying whether the code can be modified and re-released legally, and one of Py6S's 'selling points' is that it produces exactly the same results as using 6S manually. On top of that, I'm extremely inexperienced at writing Fortran, so I'm definitely not the right person to do this.

It's possibly best to try and contact the original 6S authors and explain the mismatch between the web UI/documentation and the code, and see if they want to do anything about it.

Thanks

simonrp84 commented 3 years ago

Coming back to this issue, I contacted the 6S authors and the code downloadable from the 6S website now correctly accounts for surface pressure. I've verified that the code works as expected, both when running the updated 6s manually and also by supplying py6s with the path to this executable.

Is it possible for you to now update the sixs package to reflect this code update? It is the case that py6s no longer produces the same results as running 6S manually!

robintw commented 3 years ago

Ah that's good to know. I'm impressed you got a response from the 6S authors - was this already fixed in the latest version, or have they fixed it once you reported it?

Unfortunately it will actually be quite a lot of work to update Py6S to work properly with 6S v2.1. There is no proper documentation on what has changed, and it will take a lot of looking through the Fortran code to understand where things need to be changed (eg. parameter setting values have changed in some places, sometimes in incompatible ways). In some situations Py6S works with the new version, but it is fragile and may fail in various situations.

I think there would probably need to be separate versions of Py6S for the two major 6S versions - as there are definitely people who want to stay using v1.1. There is also no way to tell which version of 6S a user has installed without running it (a --version flag or similar), which doesn't help.

So, basically, it should be possible, but would be quite a lot of work. I'm now working freelance, and am struggling to find time to devote to major extension work like this. Ideally I would be able to find some funding for this work, but I think that is unlikely. So, I will try and get to this when I can, but I can't guarantee it will be very soon.