Tag Archives: Atmospheric Correction

Batch processing atmospheric correction using ARCSI

Previous posts have covered atmospheric correction using the Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI) software. When a large number of scenes require processing then commands to automate steps are desirable and greatly simplify the process.

ARCSI provides a number of commands which are useful for batch processing:

  • arcsisortlandsat.py sorts the landsat ‘tar.gz’ files into individual directories for each sensor (i.e., LS 5 TM) and also builds a standard directory structure for processing.
  • arcsiextractdata.py extracts the contents of \verb|tar| and \verb|tar.gz| archives with each archive being extracted into an individual directory.
  • arcsibuildcmdslist.py builds the ‘arcsi.py’ commands for each scene creating a shell script.

The files used for this tutorial are shown below, but others could be downloaded from from earthexplorer and used instead.

LC80720882013174LGN00.tar.gz	LT52040232011310KIS00.tar.gz
LC82040242013139LGN01.tar.gz	LT52040242011118KIS00.tar.gz
LT52040232011118KIS00.tar.gz	LT52050232011109KIS00.tar.gz

Sort Landsat Scenes

The first command to execute is ‘arcsisortlandsat.py’ which needs to be executed from within the directory containing files:

arcsisortlandsat.py -i . -o .

This will generate the following directory structure and move the files into the appropriate directories.

> ls *
Inputs	Outputs	RAW	tmp

Inputs	Outputs	RAW	tmp

> ls */RAW
LT52040232011118KIS00.tar.gz	LT52040242011118KIS00.tar.gz
LT52040232011310KIS00.tar.gz	LT52050232011109KIS00.tar.gz

LC82040242013139LGN01.tar.gz	LC82060242015047LGN00.tar.gz

Extract Data

To extract the data from the ‘tar.gz’ files the ‘arcsiextractdata.py’ command is used as shown below:

arcsiextractdata.py -i ./LS5/RAW/ -o ./LS5/Inputs/

arcsiextractdata.py -i ./LS8/RAW/ -o ./LS8/Inputs/

Once the files are extracted the directory structure will look like the following:

> ls */RAW
LT52040232011118KIS00.tar.gz	LT52040242011118KIS00.tar.gz
LT52040232011310KIS00.tar.gz	LT52050232011109KIS00.tar.gz

LC80720882013174LGN00.tar.gz	LC82040242013139LGN01.tar.gz

> ls */Inputs
LT52040232011118KIS00	LT52040242011118KIS00
LT52040232011310KIS00	LT52050232011109KIS00

LC82040242013139LGN01	LC82060242015047LGN00

Build ARCSI Commands

To build the ‘arcsi.py’ commands, one for each input file, the following commands are used. Notice that these are very similar to the individual commands that you previously executed but now provide inputs to the ‘arcsibuildcmdslist.py’ command which selects a number of input files and generate a single shell script output.

arcsibuildcmdslist.py -s ls5tm -f KEA --stats -p RAD DOSAOTSGL SREF \
          --outpath ./LS5/Outputs --aeroimg ../WorldAerosolParams.kea \
          --atmosimg ../WorldAtmosphereParams.kea --dem ../UKSRTM_90m.kea \
          --tmpath ./LS5/tmp --minaot 0.05 --maxaot 0.6 --simpledos \
          -i ./LS5/Inputs -e MTL.txt -o LS5ARCSICmds.sh

arcsibuildcmdslist.py -s ls8 -f KEA --stats -p RAD DOSAOTSGL SREF \
        --outpath ./LS8/Outputs --aeroimg ../WorldAerosolParams.kea \
        --atmosimg ../WorldAtmosphereParams.kea --dem ../UKSRTM_90m.kea \
        --tmpath ./LS8/tmp --minaot 0.05 --maxaot 0.6 --simpledos \
        -i ./LS8/Inputs -e MTL.txt -o LS8ARCSICmds.sh

Following the execution of these commands the following two files will have been created ‘LS5ARCSICmds.sh’ and ‘LS8ARCSICmds.sh’. These files contain the ‘arcsi.py’ commands to be executed. Open the files and take a look, you will notice that all the file paths have been convert to absolute paths which means the file can be executed from anywhere on the system as long as the input files are not moved.

Execute ARCSI Commands

To execute the ‘arcsi.py’ commands the easiest methods is to run each in turn using the following command:

sh LS5ARCSICmds.sh

sh LS8ARCSICmds.sh

This will run each of the commands sequentially. However, most computers now have multiple processing cores and to take advantage of those cores we can use the GNU parallel command line tool. Taking advantage of those cores means that processing can be completed much quicker and more efficiently.

cat LS5ARCSICmds.sh LS8ARCSICmds.sh | parallel -j 2 

The switch ‘-j 2’ specifies the number of processing cores which can be used for this processing. If no switches are provided then all the cores will be used. As some parts of the process involve reading / writing large files, depending on the speed of you hard drive and the number of cores you have available you might find limiting the number of cores is needed. Please note that until all processing has completed nothing will be printed to the console. You can check the command are running using the top / htop command.

Once you have completed your processing you should clean up your system to remove any files you don’t need for any later processing steps. In most cases you will only need to keep the ‘tar.gz’ (so you can reprocess the RAW data at a later date if required) and the SREF product. It is recommended that you also retain the scripts you used for processing and a record of the commands you used for a) reference if you need to rerun the data and b) as documentation for the datasets so you know what parameters and options were used.

This post was adapted from notes for a module taught as part of the Masters course in GIS and Remote Sensing at Aberystwyth University, you can read more about the program here.

Inverting for Aerosol Optial Thickness (AOT) using ARCSI

Aerosol Optial Thickness (AOT) as an important parameter in atmospheric correction models such as 6S, which is used by the Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI) software (see previous post for more details). With time and expert knowledge it would be possible to derive an AOT value for each scene manually. In some regions of the world it has been found that a constant can be used (e.g., Australia use an AOT of 0.05; [1]).

However, for most regions a constant is not viable as the atmosphere is too variable (both temporally and spatially) and not correctly parameterising for AOT can lead to large errors in atmospheric correction [2].

ARCSI provides a method of deriving AOT from an image by using a dark object subtraction (DOS) to estimate surface reflectance in the blue channel. Taking this estimated reflectance as input, 6S is then numerically inverted to identify an AOT value which provides a surface reflectance value as close to that estimated from DOS as possible.

To execute this functionality the following command is used:

arcsi.py --sensor ls8 \
    -i Inputs/LC82040242013139LGN01/LC82040242013139LGN01_MTL.txt \
    -o ./OutputsAOTInv \
    --tmpath ./tmp \
    -f KEA --stats \
    --prodlist RAD DOSAOTSGL SREF \
    --aeroimg ../WorldAerosolParams.kea \
    --atmosimg ../WorldAtmosphereParams.kea \
    --dem ../UKSRTM_90m.kea \
    --minaot 0.05 \
    --maxaot 0.6 --simpledos 


--sensor ls8
specifies the sensor (in this example Landsat 8)

-i Inputs/LC82040242013139LGN01/LC82040242013139LGN01_MTL.txt
is the header file of the input image.

-o ./Outputs
specifies the output directory where all output file will be saved to.

--tmpath ./tmp
is a directory where temporary files can be outputted during the processing, these will be deleted afterwards.

--f KEA
specifies that the output image file format should be KEA. Any GDAL supported format can be used.

specifies that the output images should be populate with statistics and pyramids – makes display much faster (only available for KEA and HFA output formats).

--prodlist RAD TOA SREF
specifies the output products to be generated. You only need to specify what you want, e.g., SREF but if other products (e.g., RAD) are required then these will also be produced even if they are not specified.

--aeroimg ../WorldAerosolParams.kea
is an image file which is sampled to identify the aerosol model to be used for the input image (this can be downloaded from ARCSI downloads).

–atmosimg ../WorldAtmosphereParams.kea
is an image file which is sampled to identify the atmosphere model to be used for the input image (this can be downloaded from ARCSI downloads).

–dem ../UKSRTM_90m.kea
is an elevation model covering the scene, in this case the 90 m SRTM product. The higher the resolution of the DEM available the better. SRTM data are available to download from earthexplorer.

--minaot 0.05
is the lower limit of the AOT values tested when estimating the AOT value to be used to correct the scene.

--maxaot 0.6
is the upper limit of the AOT values tested when estimating the AOT value to be used to correct the scene.

specifies that a simple single valued DOS method should be used to provide an estimate of the Blue SREF used in the inversion.

This command will estimate a value for AOT, use this to parameterise and run 6S to build a LUT at elevation steps of 100 m and apply this to produce surface reflectance (multiplied by a scaling factor of 1000 to reduce the file size).

This post was adapted from notes for a module taught as part of the Masters course in GIS and Remote Sensing at Aberystwyth University, you can read more about the program here.

[1] S.S. Gillingham , N. Flood , T.K. Gill. 2013. On determining appropriate aerosol optical depth values for atmospheric correction of satellite imagery for biophysical parameter retrieval: requirements and limitations under Australian conditions; International Journal of Remote Sensing. Vol. 34(6). 2089-2100. http://dx.doi.org/10.1080/01431161.2012.738945
[2] Wilson, R.T., Milton, E.J. & Nield, J.M., 2014. Spatial variability of the atmosphere over southern England, and its effect on scene-based atmospheric corrections. International Journal of Remote Sensing. 35(13). 5198-5218. http://dx.doi.org/10.1080/01431161.2014.939781 (‘Behind the paper’ and pre-print link).

Atmospheric Correction of Landsat Data using ARCSI

The Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI) software was developed at Aberystwyth University, under contract from the Norwegian Space Centre and allows automatic correction of Landsat data to Top Of Atmosphere (TOA) or Surface Reflectance (SREF). The conversion to Surface Reflectance uses the 6S codes, accessed through the Py6S Python interface.

Getting ARCSI

If you’re running OS X or Linux, you can get ARCSI through conda. For Windows users, unless you’re very familiar with compiling software from source it is recommended you set up a virtual machine with Linux to run ARCSI.

For the latest installation instructions see the software page.

Download and untar data

This tutorial uses Landsat 8 data (although ARCSI supports other sensors). If you don’t already have a Landsat 8 scene, you can download one from http://earthexplorer.usgs.gov/.

The data is downloaded as a tarred and compressed file. After you have downloaded a scene, create a directory for the scene and copy the tar.gz file to it. Untar the file using:

tar -xf LC80430342013291LGN00.tar.gz

There are a number of tiff files (one for each each band) and associated header information in the MTL file.

After you have downloaded and untared the data create a new directory to store the output files (I’m using one called ‘OutputImages’).

Conversion to Radiance

Landsat data are provided as digital number (DN) images, the first step is to calibrate to radiance, which is the amount of energy received by the sensor per second (W) per steradian (sr) per square metre (m^{2}). This is written as W sr^{-1}m^{-2} or W sr^{-1}m^{-2}nm^{-1} for a given wavelength (in nm). This can also be defined as “The power passing through a unit area in a unit solid angle about the normal to the area (per unit spectral interval)”.

For Landsat 8, a linear equation is used to convert from DN to radiance (L):

L = (\text{DN} \times G) + I.

The gain (G) and intercept (I) are different for each band and are provided in the Landsat header file.

To convert each band to radiance in ARCSI, and stack all bands together, the following command is used:

arcsi.py -s ls8 -f KEA --stats -p RAD -o ./OutputImages \
  -i LC80430342013291LGN00/LC80430342013291LGN00_MTL.txt

Where ‘-o’ is the output directory and ‘-i’ the input header (MTL) file. Statistics and overviews (‘–stats’) are calculated for fast display.

Once this has run open the output file (ending in ‘_rad.kea’) in TuiView.

Top of Atmosphere Reflectance

When spectral radiance values are compared within and between sensors, variations frequently occur because of different sun-target-sensor configurations. Hence, these data are often converted to Top of Atmosphere (TOA) or planetary reflectance or at sensor reflectance as this eliminates the effect of different solar zenith angles at the time of image acquisition, minimises the solar irradiance band differences occurring between different sensors accounts for differences in the Earth-Sun distance. Conversion to TOA reflectance is achieved using:

\rho_{\lambda} = \frac{\pi L_{\lambda}d^2}{\text{ESUN}_{\lambda}cos(\theta_s)},

where \lambda is the wavelength, \rho_{\lambda} is the spectral (planetary or TOA) reflectance for wavelength \lambda, L_{\lambda} is the spectral radiance, d is the Earth-Sun distance (in astronomical units), ESUN is the mean solar exoatmospheric irradiance and \theta_s is the solar zenith angle.

TOA reflectance does not take into account any atmospheric effects and is simply a ratio of the at-sensor radiance with the incoming radiance from the sun, where the distance and angle of the sensor to the target are considered. The approach is easy to apply, as no external information is required.

To calibrate Landsat 8 data to TOA reflectance in ARCSI the following command is used:

arcsi.py -s ls8 -f KEA --stats -p RAD TOA -o ./OutputImages \
   -i LC80430342013291LGN00/LC80430342013291LGN00_MTL.txt

Once this has run, open the output file (ending in ‘_toa.kea’) in TuiView and compare with the radiance image (in a separate window). Note the output reflectance, which has values 0 – 1, is multiplied by a scaling factor of 1000 so the results can be stored as a 16 bit integer.

Surface Reflectance

To get to surface reflectance, an atmospheric model can be used to account for the effects of scattering and absorption in the atmosphere. Within ARCSI the 6S radiative transfer model is used. A number of parameters are required for 6S:

  • Atmospheric Profile
  • Aerosol Profile + Aerosol Optical Depth (AOD)
  • Ground Reflectance
  • Geometry
  • Altitude
  • Wavelength

Many of these parameters are calculated within ARCSI or derived from the header file (e.g., geometry, altitude and wavelength). There are a number of preset options for atmospheric profile, aerosol profile and and ground reflectance in ARCSI, to show all available options type:

arcsi.py -h

As an example we’ll run using the standard input parameters:

arcsi.py -s ls8 -f KEA --stats -p RAD SREF --aeropro Maritime \
  --atmospro MidlatitudeSummer --aot 0.25  -o ./OutputImages \
  -i LC80430342013291LGN00/LC80430342013291LGN00_MTL.txt

This will take a bit longer to run through. Once it’s finished open the top of atmosphere and surface reflectance images in linked TuiView displays and compare the profiles.

Comparison of Landsat data corrected to top of atmosphere (left) and surface reflectance (right).

Additional tasks

  1. Try changing the values of AOT (–aot) between 0 and 1, remembering to create separate output directories and observe the changes in reflectance over the range of wavelengths.
  2. You can investigate the different settings available in 6S, without applying to an image through Py6S, the quick start guide is a good place to start.

Getting help with ARCSI

The documentation for ARCSI is available at http://rsgislib.org/arcsi/

If you have any questions about ARCSI which are not covered by the documentation the RSGISLib mailing list on Google Groups is used for ARCSI support (rsgislib-support@googlegroups.com). Note you need to be a member of the group to post.

This post was adapted from notes for a module taught as part of the Masters course in GIS and Remote Sensing at Aberystwyth University, you can read more about the program here.