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…]; 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.

Py6S: Run a radiative transfer model with a user-defined input spectrum

Version 1.2 of Py6S has recently been released (including a couple of bug fix releases, taking the most recent version to v1.2.2), and the major new feature in this version is the ability to use any spectrum as the ground reflectance for a model run. Previously, users were restricted to using the built-in 6S ground spectra (vegetation, clear water, lake water and sand) or doing some complicated parameterisation to set the user-defined reflectances up in the correct manner – so complicated that it was almost never worth doing!

Luckily, that has all changed in this version – read on to find out more…

Setting the Ground Reflectance

Setting the ground reflectance parameterisation in Py6S is very simple – you just use one of the methods of the GroundReflectance class. For example, to create a 6S model and set a constant reflectance of 0.5 for all wavelengths, simply run:

from Py6S import *
s = SixS()
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.5)
s.run()

In previous versions the parameter given to HomogeneousLambertian (or HeterogeneousLambertian) could be any of:

• A single reflectance value as a float, which is used as a uniform reflectance for all wavelengths: 

  s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.7)

• A constant stating a built-in 6S reflectance value:

  s.ground_reflectance = GroundReflectance.HomogeneousLambertian(GroundReflectance.GreenVegetation)

• An array giving a user-defined spectrum given in steps of 2.5nm for the wavelength range used for that simulation (the last bit is important!)

  s.ground_reflectance = GroundReflectance.HomogeneousLambertian([0.6, 0.8, 0.34, 0.453])

But now there is one more option:

• A two-dimensional array containing wavelengths (column 0, in micrometres) and reflectance values (column 1, as a fraction). This array will then be taken by Py6S, resampled to 2.5nm, subset to only the range of wavelengths being used for the current simulation, and used with the underlying 6S model.

  s.ground_reflectance = GroundReflectance.HomogeneousLambertian(spectrum_array)

This saves a lot of time and extra code, compared to doing all of the resampling yourself and using the user-defined spectrum option which was always available in Py6S. However, there are more benefits than this… Py6S will store the whole spectrum that you give, and just re-sample it to the correct wavelength range every time you run the model, which means you can use the lovely helper methods like run_vnir and run_landsat_tm to run simulations for a spectrum for a number of wavelengths with very little code. For example:

from Py6S import *
s = SixS()
# Somehow generate a 2D array called spectrum_array here
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(spectrum_array)
wavelengths, results = SixSHelpers.Wavelengths.run_vnir(s, output_name="apparent_radiance")

Getting spectral data

Of course, that example won’t work until we add in some code to generate this 2D array for us. There are a two ways we could get this array:

• Generate it within our code somehow – for example, using a simple model, or interpolation between some hard-coded values.
• Load it from some sort of external datafile, for example from a CSV file using np.loadtxt

Both of those are ways that you might want to use, but I’ve added functions to Py6S to simplify a couple of ways of doing this that may be useful to people.

Firstly, the functions in the Spectra module make it easy to import spectra from the USGS Spectral Library and the ASTER Spectral Library. These two libraries seem to be the most commonly used – as far as I am aware, at least – and have a very wide range of spectra in them. The two functions are very simple to use – you simply pass the URL or file path to the spectrum data file, and it returns the right sort of array. Thus you can do:

spectrum_array = Spectra.import_from_usgs("http://speclab.cr.usgs.gov/spectral.lib06/ds231/ASCII/V/russianolive.dw92-4.30728.asc")
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(spectrum_array)

Or

s.ground_reflectance = GroundReflectance.HomogeneousLambertian(Spectra.import_from_usgs("http://speclab.cr.usgs.gov/spectral.lib06/ds231/ASCII/V/russianolive.dw92-4.30728.asc"))

(Either of the URLs could be replaced with a simple local file path instead, which is particularly useful if you have downloaded the USGS Spectral Library archive file, or obtained the CD of the ASTER Spectral Library).

The other way to get data is to use a model. Many models can produce spectral data, and most will output in a format like CSV that can be imported into a NumPy array really easily – but to make things even easier I have also released a Python interface to the ProSAIL called – you guessed it – PyProSAIL. Again, the methods have been carefully designed to return the right sort of arrays, so you can simply do:

# Make sure you have both PyProSAIL and Py6S installed
import pyprosail
from Py6S import *

spectrum = pyprosail.run(1.5, 40, 8, 0, 0.01, 0.009, 1, 3, 0.01, 30, 0, 10, 0, pyprosail.Planophile)
s = SixS()
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(spectrum)
s.run()

So, hopefully that has been a useful introduction to how to use user-defined spectra with Py6S. I’ll leave you with an example of how to put all of this together to do something useful – that is, to run a Py6S simulation with the ground reflectance set to the spectra from every USGS spectral library file in a directory, and store the results:

from Py6S import *
import glob
import numpy as np

# Create the SixS model and set some parameters
s = SixS()
s.aero_profile = AeroProfile.PredefinedType(AeroProfile.NoAerosols)
s.altitudes.set_sensor_satellite_level()
s.altitudes.set_target_sea_level()

# Get a list of files that have the extension .asc
# (the extension used by the USGS spectral library files)
files = glob.glob("./USGSSpecLib/Selected/*.asc")

results = []

# Run the model using each filename as the ground reflectance spectrum
# using the Landsat ETM+ bands
for filename in files:
	s.ground_reflectance = GroundReflectance.HomogeneousLambertian(Spectra.import_from_usgs(filename))
	wvs, res = SixSHelpers.Wavelengths.run_landsat_etm(s, output_name='apparent_radiance')
	results.append(res)

# Stack all of the results into an array
results_arr = np.vstack(results)

# Write the results out to a CSV file for further analysis
np.savetxt("results.csv", results_arr, delimeter=",")