Robin's Blog

Can I atmospherically-correct my images with Py6S?

I’m a big fan of Matt Might’s blog, and thought I’d implement one of his tips for blogging as an academic – namely Reply to Public. I’ve had a number of emails from Py6S users asking me questions about how to atmospherically-correct entire satellite images with Py6S – so I thought ‘d respond online, so that in the future people either find this response without needing to email me (ideally), or at the very least, I can point them to this blog post when they do email me.

Unfortunately, the simple answer to the question is: Py6S cannot (yet!) atmospherically correct satellite images – at least, not in a sensible manner – so I’d suggest using other atmospheric correction software such as ATCOR or FLAASH. The longer answer is below…

If you read the Py6S documentation you’ll find that there is an atmospheric correction option which you can use to take an at-sensor radiance (often called a Top of Atmosphere radiance when you’re dealing with satellite sensors) and atmospherically-correct it to get a corrected surface reflectance (or radiance) value. You could write some code to do that – it might look a bit like this:

from Py6S import *
# Create a SixS object called s
# (used as the standard name by convention)
s = SixS()

# Set the atmospheric conditions as usual
s.aero_profile = AeroProfile.PredefinedType(AeroProfile.Maritime)
s.aot550 = 0.05
s.atmos_profile = AtmosProfile.UserWaterAndOzone(2, 0.318)

# Set the wavelength
s.wavelength=Wavelength(PredefinedWavelengths.MODIS_B6)

# Set the altitudes
s.altitudes.set_target_sea_level()
s.altitudes.set_sensor_satellite_level()

# Set the geometry
s.geometry = Geometry.User()
s.geometry.solar_z= 35.8
s.geometry.solar_a= 149.4
s.geometry.view_z= 5.1
s.geometry.view_a= 106.8
s.geometry.month=06
s.geometry.date=28

# Turn on atmospheric correction mode and set it to do the
# correction assuming a Lambertian surface with a TOA
# radiance of 137.5 W/m^2
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromRadiance(137.5)   

# Run the model
s.run()

# Print the result of the atmospheric correction
# (assuming Lambertian reflectance)
# This is the ground-reflectance for the given radiance,
# under the given atmospheric conditions
print s.outputs.atmos_corrected_reflectance_lambertian

This works fine, and so you might think that all you need to do to correct a satellite image is to take the code above, and run it for each pixel of the image using a for loop, as in the pseudo-code below:

# Load in the image (eg. using GDAL)

for pixel in image:
    reflectance = run_py6s(various, parameters, here)

    # Store reflectance in the output image

# Save output image

Unfortunately, that won’t work very well. Or rather, it will work – it’ll just take a long time! Let’s do some maths to work out roughly how long:

• As a rough approximation, it takes Py6S around two seconds to atmospherically correct a single pixel in a single band (the exact time depends on how complicated the parameterisation is – for example, a full manual AERONET-based aerosol parameterisation will take longer than a standard aerosol model)
• A full Landsat scene is roughly 5500 x 6000 pixels (170km x 185km), which is around 33 million pixels
• 33 million pixels x 2 seconds per pixel = 66 million seconds, or 763 days!
• If I started running this now – on one band of a Landsat scene – it’d be finished around April 2015. Even if I used a simple parameterisation that only took one second per pixel, it’d still take over a year!

Now, I suspect you’re now thinking the Py6S is really awful because it takes such a long time. Well, unfortunately it’s not down to Py6S (if it was, I could try and improve it!). In fact, the Python bit of Py6S adds very little time – the vast majority of the time is taken by 6S itself.

So, maybe it’s 6S that is rubbish. Well, again, unfortunately not, or we could switch to something that is better. All other Radiative Transfer Models, such as MODTRAN, SCIATRAN etc, suffer from the same problem.

So, how on earth do people manage to correct satellite images? Obviously people have managed to do it – and there is commercial software available to do it. Well, they use something called a Lookup Table (LUT). This involves running the Radiative Transfer Model many times to produce a table of results, and then looking up the result you want in this table. Creating the table takes a long time (but nowhere near as long as 763 days!) and it often takes a lot of space to store the resulting table (for example, the Lookup Table used in ATCOR is a multi-Gigabyte file), but once you’ve got the table you can correct a pixel in a tiny fraction of a second – rather than one or two seconds – which means the correction of an image is a lot quicker.

A lookup table is a good example of two standard trade-offs in computing:

• Setup time vs Run time: Often algorithms can either be designed to have a very short (or no) setup time, but then take a long time to run; or they can be designed to have a long setup phase, and then run very quickly from then onwards. Often the decision of which one to focus on depends on how often you’re going to run your process – if you want to run it multiple times (like an atmospheric correction algorithm) then it is normally better to have a long setup time (which you only have to do once!) and then a short run time.
• Memory vs Time: Again, algorithms often trade speed against memory – a fast algorithm normally requires more memory, and vice-versa. In this case, a Lookup Table trades speed for memory – we have to store the Lookup Table in memory, but once we have access to it, we can do corrections very fast.

I haven’t written any code to use Py6S to create a lookup table and then use that to correct an image (although I will almost certainly write some code to do this sometime during my PhD) – but it wouldn’t be too difficult to do. The procedure would be something like the following:

Setup Phase:

1. Come up with a list of parameter values for all of the important parameters – choosing a number of parameters across their realistic range (there is a trade-off between accuracy and speed here, as always). For example maybe: AOT = [0.1, 0.2, 0.3, 0.4, 0.5]; solar_z = [10, 20, 30, 40, 50]; radiance = [10, 20, 30, 40, 50, 60, 70] etc.
2. Write some code to loop over all of these lists and run every possible combination of them
3. Store the result of all of these runs of Py6S (ie. the atmospherically-corrected reflectance) and the parameters used to get the result in a big table of some sort.

Correction Phase:

Loop through every pixel in the image, and for each pixel:

1. Get the parameters (AOT, angles, radiance etc) for that pixel
2. Choose the nearest parameters that were used to create results in the table (for ALL of the parameters you used)
3. Interpolate between these parameters as needed (eg. for an AOT of 0.15, get the results for an AOT of 0.1 and and AOT of 0.2 and interpolate) to get an interpolated ground reflectance.
4. Store this reflectance in another image.

This basic procedure has been used in almost all atmospheric correction tools for the last thirty years, and is described in various papers including Fraser et al. (1992) and Liang et al. (1997). Another way to speed-up the computation even more is by introducing parallel computation. If we can assume that each pixel’s atmospheric correction is entirely independent from every other pixel’s (which may not be the case in the real-world, but is generally assumed by most atmospheric correction algorithms) then we can split the pixels between a number of processors and thus correct many pixels in parallel – which gets an even greater speedup as we can do this both when generating the LUT, and when correcting the image.

Unfortunately, actually implementing this is a bit more complicated than I’ve explained here and various decisions have to be taken (including what ranges of parameters to use, how to interpolate, how best to store the multi-dimensional lookup table and more…) – but it is definitely on my ‘to do’ list as part of my PhD, and as part of Py6S.

So, in conclusion: Py6S can’t do atmospheric correction of satellite imagery in a sensible time at the moment, but should be able to within the next few years. In the meantime, I suggest using other atmospheric correction software (most of which, unfortunately, is commercial software). You may think that not being able to do this makes Py6S useless…but, as my next post will show, Py6S is still very useful for remote-sensing research.

References

Fallah-Adl, Hassan, et al. "Fast algorithms for removing atmospheric effects from satellite images." Computational Science & Engineering, IEEE 3.2 (1996): 66-77.

Fraser, R. S., Ferrare, R. A., Kaufman, Y. J., Markham, B. L., & Mattoo, S. (1992). Algorithm for atmospheric corrections of aircraft and satellite imagery.International Journal of Remote Sensing, 13(3), 541-557.

Liang, Shunlin, et al. "An operational atmospheric correction algorithm for Landsat Thematic Mapper imagery over the land." Journal of Geophysical Research 102.D14 (1997): 17173-17.

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, Programming, Py6S, Python, Remote Sensing