Atmospheric correction software for Landsat5 Thematic Mapper data set

Eric Vermote

Draft Outline

TABLE OF CONTENTS

I Summary 5

II Cloud screening procedure 5

III Atmospheric correction procedure 7

III.1 Intrinsic reflectance and transmission effects 7

III.1.1 Theoretical background 7

III.1.2 Software Implementation 8

III.2 Correction for target non-uniformity 9

III.2.2 software implementation 10

IV ATMOSPHERIC PARAMETERS 25

IV.1 Pressure 25

IV.2 Ozone 25

IV.3 Water Vapor 25

IV.4 Aerosol 26

V OUTPUT PRODUCTS 27

V.1 corrected reflectances 27

V.2. QC 27

V.3 log file 28

V.3.1 Input parameters 28

V.3.2 output parameters 29

VI FUNCTIONAL DESCRIPTION 31

VI.1 Check inputs 31

VI.2. Clouds, shadow mask aerosol statistics 31

VI.3 Optical depth inversion 31

VI.4 interpolation of aerosol grid 32

VI.5 atmos. corr. (step 1) grid creation, adjacency effect filter creation 33

VI.6 atmos. corr. (step 1) 33

VI.7 atmos. corr. (step 2) 33

VI.8 generation of process report 34

VII DETAILED CODE DESCRIPTION 36

VII.1 general processing chart 36

VII.2 Global variables, definition 37

VII.2.1 definitions 37

VII.2.2 structures 40

VII.2.2 global variables 41

VII.3 main program 43

VII.4 C routines 45

VII.4.1 void compsums(int k,int tmpint, u_char *qcbuf[],int block) 45

VII.4.2 void fillgaot(float aot[],float sigmaaot[]) 46

VII.4.3 void interaot(float aot[]) 47

VII.4.4 void atm_cor(u_char *buf,u_char *qc,int ib,int bufsize,int firstline, … 48

VII.4.5 int prepwb(u_char *aebuf1,u_char *aebuf2,u_char *aebuf3,u_…. 49

VII.4.6 int ae_cor(u_char *aebuf1,u_char *aebuf2,u_char *aebuf3,int i… 59

VII.4.7 float getpressure(u_char *dem,double dlat,double dlon) 51

VII.4.8 float getwv(u_char *wvcontent,double dlat,double dlon,short month) 51'

VII.4.9 int initializepars(char *file,params_t *inpars) 52

VII.4.10 int rdparamsline(FILE *fd,int *par,int *ind,char *value) 53

VII.4.11 int nbofvalues(char *str) 54

VII.4.12 void print_pars(FILE *fd,params_t *inpars) 55

VII.4.13 void process_error(int level,char *msg) 55

VII.4.14 void logmsg(FILE *fdlog,char *msg) 56

VII.4.15 void cloud_mask(u_char *b3,u_char *b4,u_char *b6,… 56

VII.4.16 float comp_dsol(short jday) 58

VII.4.17 short julian_day(short year,short month,short day) 59

VII.4.18 int getGP(int i,int j,int bh,int bw) 60

VII.4.19 void initcoefsH2O(float *aH2O,float *bH2O,float *cH2O) 60

VII.4.20 void initcoefsO3(float *aO3) 61

VII.4.21 initcoefsCO2(float *aCO2,float *bCO2) 61

VII.4.22 initcoefsO2(float *aO2,float *bO2) 62

VII.4.23 void initcoefsNO2(float *aNO2,float *bNO2) 63

VII.4.24 initcoefsCH4(float *aCH4,float *bCH4) 63VII.5 fortran routines 64

VII.5.1 subroutine aeroso (iaer,co,xmud) 64

VII.5.2 subroutine bdm 66VII.5.3 subroutine caluoz(uo3,jour,xlat) 67

VII.5.4 subroutine chand (xphi,xmuv,xmus,xtau,xrray) 67

VII.5.5 subroutine compgrac(ttaer,tts,ttv,tphi,uwv,uoz,pres,ib,ttgog,ttg…68

VII.5.6 subroutine compgrao(tmref,tsref,smref,ssref,ts,tv,phi,uwv,uoz,pres…. 70

VII.5.7 subroutine csalbr(xtau,xalb) 72

VII.5.8 real function fintexp3(xtau) 73

VII.5.9 real function fintexp1(xtau) 74

VII.5.10 subroutine discre(ta,ha,tr,hr,it,nt,yy,dd,ppp2,ppp1,zx) 74

VII.5.11 subroutine dust 75VII.5.12 subroutine gauss(x1,x2,x,w,n) 76

VII.5.13 subroutine iso(tamoy,trmoy,pizmoy,tamoyp,trmoyp,palt,nt,mu,rm,gb,xf) 77

VII.5.14 subroutine kernel(is,mu,rm,xpl,psl,bp) 79VII.5.15 subroutine ocea 80

VII.5.16 subroutine os(tamoy,trmoy,pizmoy,tamoyp,trmoyp,palt, 81phirad,nt,mu,np,rm,gb,rp,xl) 81VII.5.17 subroutine print_error(tex) 83

VII.5.18 subroutine psprdf(pres,ttu,tdu,ib,dl,mp,alpha,nt) 84

VII.5.19 subroutine soot 85

VII.5.20 subroutine trunca(coeff) 86

VII.5.21 subroutine wate 87

VIII Glossary 88

I Summary

This document describes the software and theoritical background for atmospheric correction of Thematic Mapper data (band 1-5,7). Atmospheric correction includes correcting for Rayleigh scattering, gaseous absorption (ozone, oxygen, water vapor) and aerosol scattering. The software corrects both for intrinsic reflectance and transmittance effect as well as adjacency effects which are quite important for TM with its 30m pixel size. The products of the correction are the surface reflectances (lambertian targets are assumed) as if they were measured at ground level without atmosphere. The correction is expected to perform the best for homogeneous aerosol distribution (50km), zero cloudiness and for aeras where some vegetation is present. In all other cases, the algorithm will encounter problems because of the difficulties in assessing correctly the aerosol effect or correcting for adjacency effects produced by clouds.

The document presents first the adopted cloud screening procedure (Section II), then describes the two step atmospheric correction and implementation (Section III), the next part (Section IV) describes how the ancillary data necessary for atmospheric correction are generated. Section V describes the output image products, quality control information. and a description of the log file to be transmitted to the validation scientist. In section VI, we give a functional description of the code, the description is detailled at the routine level in section VII. A glossary summarizing the abreviations and symbols used in this document is given in section VIII.

II Cloud screening procedure

The cloud screening procedure stems from the work of Stowe et al (1991) simplified and adapted to the case of Landsat Thematic Mapper. It uses three types of considerations for the detection of clouds, their ability to reflect the radiation in the visible and near infrared (Reflectance threshold), their relative whiteness (Simple ratio, visible/near infrared test) and the that fact they are above the surface and therefore cooler than surface (Thermal threshold test). If a pixel succeeds in all three tests it is considered as a clear pixel, if it fails all three tests it is considered as a cloudy pixel, in all others cases the pixel is considered as mixed. The distinction is quite important, because no atmospheric correction is conducted on cloudy pixels but is applied to all other cases.

The following table summarize the procedure:

	Cloudy	Clear
Reflectance Threshold test	r(TM3) =>44% C1	r(TM3) <44% C'1
Thermal threshold test (1) with T(1) ECMWF (2) default	T(TM6)<T(ECMWF)-5 T(TM6)<249K C2	T(TM6)>T(ECMWF)-2 T(TM6)>293K C'2
Reflectance ratio	C3	C'3

cloudy=C1.and.C2.and C3

clear=C'1,and.C'2.and.C'3

Based on the cloud mask and T(TM6) a cloud shadow mask is generated and included in the QC mask. It is used by the aerosol retrieval module (shadow, mixed and cloudy pixels are discarded for aerosol thickness inversion).

references

Stowe, L. L., McClain, E. P., Carey, R., Pellegrino, P., Gutman, G. G., Davis, P., Long, C. and Hart, S. (1991), Global Distribution of Cloud Cover derived from NOAA/AVHRR operational Satellite Data, Adv. Space Res. 11(3):51-54.

III Atmospheric correction procedure

III.1 Intrinsic reflectance and transmission effects

III.1.1 Theoritical background

The general equation for atmospheric correction assumes in this first step that the target is lambertian and uniform and that gaseous absorption and particle scattering can be simply decoupled, that is :

(1)

where:

(2a)

(2b)

with , a,b,c depend on the wavelength (see Fig. 1 to 6 for spectral response used in each channel at the end of this section) and the gaseous absorption (see Table 1 and Figure 7,12 at the end of this section ), is the concentration of water vapor in [g/cm2] and is the concentration in ozone in [cm.atm]. The content in H2O and O3 are read from gridded ancillary data. For all other gases the concentration is equal to,P [atm], the pressure in units of atmosphere, this pressure, P [atm], is determined from the altitude, zt, given by the gridded ETOPO5 elevation model if ECMWF pressure is not available, assuming:

P=P0(1-) (3)

is computed by subroutines OS.f

(tR(P0) is given in table 1 for each channel)

and,

is computed from the subroutine ISO.f

and is also computed from subroutine ISO.f

In all the expressions above and following, tA is computed from the gridded ancillary data.(tA in TM channel 1,2 and 3). To account for the spectral variation of the aerosol optical depth, we adjust the optical depth in channel 4,5,7 using a law in l^-1 upon the central wavelength.

In future implementation we plan to refine that step by using the actual dependence measured from channel 1,2,3.

(5)

with

(6)

III.1.2 Software Implementation

From equation (1) it can be seen that inversion of the surface reflectance from at-sensor measured radiance reduce to a two step operation:

(7)

(8)

where d2 is the earth-sun distance constant for the whole TM scene and is computed from the julian day as:

(9)

Es is taken from table 1. Ccal and Coffset are calibration coefficients used for converting digital counts to at sensor radiance [W/m2/str/mm].

The computation of each term of equations (7) and (8) which has to be performed for each band is not done for each pixels of the images but only for 16 grided points corresponding to the spatial resolution of the aerosol products. The values for each pixel are interpolated using bi-linear interpolation techniques.

III.2 Correction for target non-uniformity

An additional step has been developed for adjacency effect correction. The correction procedure stems from the modelling work by Tanré et al (1981) and is simplified in this code for operational application. When taking into account the adjacency effect, the signal at the top at the atmosphere can be rewritten by decoupling the photons coming directly from the target () from those coming from areas adjacent to the target and then scattered to the sensor (td(mv)) :

(10)

where rs is the pixel reflectance and <rs> is the contribution of the pixel background to the top of the atmosphere signal that is computed as:

(11)

where x,y denote the coordinate to a local reference centered on the target, and f(r) is the environment function (see definition later). In our case the analytical integral will become:

(12)

where i and j refer to the coordinates of an array of pixels centered on the target pixels.

In the previous step, we obtained a corrected reflectance, that we will now call, , that stands for surface reflectance with adjacency effect. We can see from (1) that if we consider that: , then rs and are related through the following equation:

(13)

therefore, we can correct the reflectance obtained in the previous step, for adjacency effects using:

(14)

with , with t=tA+tR

III.2.2 software implementation

In practice, <rs> is computed from a sub-zone of the original image centered on the pixel to be corrected using:

(15)

with r(i,j) representing the distance between the pixels (i,j) and the center of the zone that is approximated by:

:

(16)

where Dx and Dy are taken to be equal to 28.50m for rorg which considered the short range adjacency effect (pixels are taken from the full resolution image), and Dx, Dy equal to 280m in the ravg terms which representents the remaining contribution (long range), for which an average contribution is considered (the pixels are averaged by box of 10x10).

Finally, knowing r(i,j), f(r) (with r in [km]) can be computed using:

(18)

with

(19a)

(19b)

and

(20a)

(20b)

references

Tanré, D., Herman, M. and Deschamps, P. Y. (1981), Influence of the background contribution upon space measurements of ground reflectance, Applied Optics 20:3673-3684.

Putsay, M. (1992), A simple atmospheric correction method for the short wave satellite image, Int. J. Remote Sensing 13(8):1549-1558.

Table 1: Constants for each TM channel used in correction software

Figure 1: Spectral response of TM band 1 used to generate table 1 constants.

Figure 2: Spectral response of TM band 2 used to generate table 1 constants.

Figure 3: Spectral response of TM band 3 used to generate table 1 constants.

Figure 4: Spectral response of TM band 4 used to generate table 1 constants.

Figure 5: Spectral response of TM band 5 used to generate table 1 constants.

Figure 6: Spectral response of TM band 7 used to generate table 1 constants.

Figure 7: Comparison of approximation (Eq 2a) and exact transmission (6S) computation for water vapor for each TM bands under various conditions.

Figure 8: Comparison of approximation (Eq 2b) and exact transmission (6S) computation for ozone for each TM bands under various conditions.

Figure 9: Comparison of approximation (Eq 2a) and exact transmission (6S) computation for carbon dyoxide for each TM bands under various conditions.

Figure 10: Comparison of approximation (Eq 2c) and exact transmission (6S) computation for oxygen for each TM bands under various conditions.

Figure 11: Comparison of approximation (Eq 2c) and exact transmission (6S) computation for nitrous oxyde for TM band 7 under various conditions.

Figure 12: Comparison of approximation (Eq 2c) and exact transmission (6S) computation for methane for each TM bands under various conditions.

IV ATMOSPHERIC PARAMETERS

IV.1 Pressure

The pressure P is fixed for the entire scene and is either read from ECMWF data (the closest measurement) or computed from elevation model ETOPO5 (resolution reduced to 1°x1°) assuming a pressure at sea level of 1013mb . The pressure should be in unit of millibars.

IV.2 Ozone

The ozone content should be in cm.atm, (Dobson/1000). It is expected to come from TOMS database. In case, no data are available a negative value should be entered to the correction program, the correction program will in this case use climatologie of ozone content based on London et al (1976)

IV.3 Water Vapor

A profile of Pressure,Temperature and relative humidity is needed as input for the computation of integrated water content. The ECMWF profile closest to the TM scene is used for that purpose. The pressure is expected to be in millibars, the temperature in Kelvin and the relative humidity, Rh, is a unit-less number in % between 0% (no water vapor, for example above 10km altitude) and 100% (saturated in water vapor).

The integrated water vapor content, , in g/cm2 can be computed from specific humidity Uspec (g/g or kg/kg), which represents the ratio between the mass of water vapor at that level and the mass of air, and pressure (Millibar) by:

The specific humidity can be computed from the relative humidity, Rh, through the dew point temperature, Tdew point. The dew point temperature is obtained at each level from

with a=0.01, b=2.858487 and ti values are given in the table below:

and the specific humidity is then computed by

IV.4 Aerosol

The aerosol product is computed for 16 grid points equally space on the original image (4 lines, 4 collums) in three spectral bands (1,2,3). Based on the detection of dark targets in each grid cells. If in some of the grid cells dark targets are absents, the values derived in others grid cells are used to generate the missing data. In the extreme case where no dark targets can be found the aerosol optical depth are set to zero and therefore only a rayleigh correction is applied.

references

London, J., Bojkov, D. R., Oltmans, S. and Kelly, J. L., 1976, Atlas of the Global Distribution of the Total Ozone July 1957-June 1967 NCAR/TN/113+STR, NCAR,Boulder,Colorado,USA.

V OUTPUT PRODUCTS

V.1 corrected reflectances

The output corrected reflectances (between 0 and 1) are scaled accordingg to the dynamic range of the instrument to be encoded on 1 byte. Valid range is between 1 and 254, each value lower than 1 is set to 0 , each value greater than 254 is set to 255 and in both cases the bad data flag is activated (see QC). The scaling factor for each band are provided in the log file.

V.2. QC

In addition to the sixs band of surface reflectance a Quality Control band (encoded on 1 byte) is produced that describes the corrections that have been done on the data and their integrity. The QC carry the cloudy/non cloudy information.

V.3 log file

In order to provide quality check, the aerosol product, the water vapor, ozone and some statistics about clouds have to be gathered in an ascii file and automatically sent with the input parameter file to the validation scientist (eric@kratmos,gsfc.nasa.gov). If inconsistencies are detected by the validation scientist original data should be transfered for further investigations.

V.3.1 Input parameters

Calibration coefficients for each band: slope, offset to convert from digital count to [W/m2/str/mm].

example for band 1:

CAL1=1.04325E-03,1.82985

Date: Year-Month -day

example June, 3, 1994

date=940603

III) ozone in [cm.atm] or dobson/1000

UOZ=0.280

IV) ECMWF water vapor, pressure and surface temperature

water vapor content [g.cm-2]

UWV=2.4

surface pressure [millibar]

ECMWFSP=1013

surface temperature [Kelvin]

ECMWFST=298

V) Image size

NUMBER OF LINES=5965

NUMBER OF PIXELS=6967

VI) Grid point data (16), for each point line,,collumn latitude,longitude, solar zenith angle, solar azimuth angle, view zenith angle, view azimuth angle.

example the grid point near the upper left corrner of the image at line number 745, collumn 870, corresponding to a latitude of 45.00 deg north and longitude of 10.00 deg west with sun zenith angle of 32.63 deg, sun azimuth of 270.00 deg, view zenith angle of 7.38 deg and view azimuth of 90.00 deg.

GP1=745,870,45.00,-10.00,32.63,270.20,7.38,90.00

......etc

V.3.2 output parameters

pressure, ozone and water vapor used for correction:

PSURF=1013

UOZ=0.280

UH2O=2.4

cloud statistics, percent clear,cloudy,mixed, shadow, bad data [%]

CLOUD=10

MIXED=5

SHADOW=1

BADDATA=0.2

C) aerosol optical thickness grid points, in three bands for each grid point+standart deviation in TM band 1 (if 0 then grid point was interpolated).

AOT1=0.65,0.40,0.35,0.01

AOT2=0.65,0.40,0.35,0

.... etc

VI FUNCTIONAL DESCRIPTION

VI.1 Check inputs

Check the presence in the input parameters files of all informations needed: calibration coefficients, Date,Grid of lat-lon,angles, surface temperature, surface pressure water vapor content, ozone content.

If calibration coefficients or date or grid points are missing, the process stops.

if any of the others parameters are missing, default values are provided through climatogoly.

Input: Input parameter file (see section 5.)

Output: valid parameters or status report for job failure.

VI.2. Clouds, shadow mask aerosol statistics

Apply cloud masking algorithm based on channel 3,4 reflectance and 6 temperature, generate cloud shadow mask based on cloud apparent temperature (and background temperature if surface temperature is missing). (see Section 1.). Generate mean and standard deviation of TM channel 1,2,3,7 for each pixels classified as DDV for each grid-cells (centered on gridpoints). Also compute statistics on cloud,shadow, bad pixels percentages.

Input: TM band 1,2,3,4,6,7

Output: aerosol input statistics, temporary cloud+shadow mask,cloud+shadow+qc stats

VI.3 Optical depth inversion

Compute aerosol optical depth at three wavelength based on comparison between computed top of the atmosphere reflectance in channel 1,2,3 (using channel 7 as predictor of surface reflectance) and observed top of the atmosphere reflectance.

Input: statistics for aerosol generated by (2), grid points geometrical conditions

Ouput: optical depth in 3 channels.

subroutines used:

aerosol.f: compute aerosol characteristics for a given model (phase function,single scattering albedo).

trunca.f: decompose the aerosol phase function in series of Legendre polynomial used in OS.f and ISO.f and compute truncation coefficient, f, to modify aerosol optical thickness, t , and single scattering albedo, w0 according to legendre transformation.

gauss.f: Compute for a given n, the gaussian quadrature (the n gaussian angles and the their respective weights). The gaussian quadrature is used in numerical integration involving the cosine of emergent or incident direction zenith angle (used by iso.f,os.f)

kernel.f: compute for the value of Legendre polynomials used in the successive order of scattering method.

discre.f: Decompose the atmosphere in a finite number of layers. For each layer, it provides the optical thickness, the proportion of molecules and aerosols assuming an exponential distribution. This layer decomposition is used in the integrating scheme of iso.f and os.f.

iso.f: compute atmospheric transmission,T.

os.f: compute atmospheric intrinsic reflectance.

VI.4 interpolation of aerosol grid

The aerosol grid produce by 2 may contains 0. in case no good pixels were found in the grid cells. It is therefore interpolated by the bilinear interpolation method starting with a distance of 1 grid point to the missing value, and increasing the distance till all the values are filled. In case, not even one aerosol value was different than zero, the aerosol optical depth is set to zero for the whole grid.

input: raw aerosol optical depth grid

output: interpolated aerosol optical depth grid

VI.5 atmos. corr. (step 1) grid creation, adjacency effect filter creation

Based on grid point geometrical conditions, aerosol optical depth, surface pressure, water vapor, ozone content, atmospheric correction parameters (see Eq. (1) section 2) are computed for each grid points. Subroutines iso.f and os.f as well as all the other mentionned in process 3. are used in this process.

Based on average aerosol optical depth over the scene and assuming a nadir view, a precalculated filter function for adjacency effect is computed.

input: interpolated aerosol optical depth grid, ozone,water vapor content,pressure,geometric conditions grid

output: atmospheric parameters grid, filter function for adjacency effect.

VI.6 atmos. corr. (step 1)

For each pixels, atmospheric correction parameters are interpolated from values computed at the grid points. A buffer of corrected values to be used by adjacency effect correction module (7) that works on an array of data is created.

input: TM band 1,2,3,4,5,7, atmospheric correction grid,calibration coefficients, geometrical grid.

ouput: buffer of corrected values

VI.7 atmos. corr. (step 2)

adjacency effect correction (section 2.2) is performed on buffer of step 1 corrected values. Step 2 corrected values are scaled (x400) converted to one byte and written to the ouput files.

input: buffer of corrected values (step 1).

ouput: output file of scaled corrected values (step 2)

VI.8 generation of process report

Statistics and selected buffers produced during the process are gathered and written to a parameters output file.

input: valid parameters for atmospheric correction, aerosol grid,atmospheric correction grid,cloud/shadow/bad data statistics.

ouput: ascii output file of input values.

disks file:

(1) climatology, (2) inputs provided by ESRIN, (3) output, (4) temporary.

buffer in memory.

Processes:

(1) quick, (2) I/O consuming (3) CPU consuming (4) I/O,CPU consuming

VII DETAILED CODE DESCRIPTION

VII.1 general processing chart

To be completed soon: will give the tree structure (using gprof utility).

VII.2 Global variables, definition

VII.2.1 definitions

GENERALS

u_char: is a redefinition of type unsigned char

M_PI: redefine if not defined already to pi (3.1415etc)

DGTORAD: conversion factor from degree to radians

LEAPYR(y): define function of y to check if y is bissectile year

MARCH: index of month of march (to compute julian day)

ERRLEVEL_FATAL: fatal error

ERRLEVEL_NORMAL: error (non fatal)

ERRLEVEL_WARNING: warning

GRID POINT/CELL

GPS_COL: number of grid point/cell accross the scene

GPS_LINE:number of grid point/cell along the scene

NB_GPS: total number of grid point/cell (=GPS_COL*GPS_LINE)

PROCESSING BUFFER SIZE

QCBUF_SIZE: number of lines for input and output buffers in

step 1 processing.

AEBUF_SIZE: number of lines for input and output buffers in

step 2 processing.

PLANK LAW CONVERSION TO TEMPERATURE

K1: coefficient used for conversion of the radiance in TM

band 6 to temperature (Plank's law)

K2: coefficient used for conversion of the radiance in TM

band 6 to temperature (Plank's law)

TMIN6: minimum temperature (in K) used for scaling of band 6

temperature (option -noatmos only)

TMAX6: maximum temperature (in K) used for scaling of band 6

temperature (option -noatmos only)

WATER VAPOR CLIMATOLOGICAL FILE/ DEFAULT VALUE

WVFILE: name of the file

WV_NBLAT: number of cell in latitude

WV_DLAT: size of a cell in latitude

WV_LATMIN: first latitude

WV_LATMAX: last latitude

WV_NBLON: number of cell in longitude

WV_DLON: size of a cell in longitude

WV_LONMIN: first longitude

WV_LONMAX: last longitude

WV_DFTVALUE: default value for water vapor (in g/cm2).

DIGITAL ELEVATION MODEL FILE/ PRESSURE DEFAULT VALUE

DEMFILE: name of the file

DEM_NBLAT: number of cell in latitude

DEM_DLAT: size of a cell in latitude

DEM_LATMIN: first latitude

DEM_LATMAX: last latitude

DEM_NBLON: number of cell in longitude

DEM_DLON: size of a cell in longitude

DEM_LONMIN: first longitude

DEM_LONMAX: last longitude

P_DFTVALUE: default value for pressure (in millibars).

AEROSOL OPTICAL DEFAULT VALUES

AOPTB1_DFTVALUE: aerosol optical depth in TM band 1, default

value.

AOPTB2_DFTVALUE: aerosol optical depth in TM band 2, default

value.

AOPTB3_DFTVALUE: aerosol optical depth in TM band 3, default

value.

ADJACENCY EFFECT CORRECTION SIZING

PSPMAT_RAD: radius of the correction in pixels.

MPSP_NBELT: number of elts of the correction mattrix (2*PSPMAT_RAD+1)*(

2*PSPMAT_RAD+1)

TWOPSMAT_RAD: twice the correction radius

PSPMAT_NLC: number of line/collumn of the correction mattrix.

CALIBRATION CONSTANT

float LowerRadLimit[7]: offset for conversion input data to radiance

float UpperRadLimit[7]: highest radiance for scaling input data

float revgain[7]: revised gain derived from SCAR-A experiment

static float LowerRadLimit[7]={-1.5, -2.8, -1.2, -1.5, -0.37, 0.0, -0.15};

static float UpperRadLimit[7]={152.1,296.8,204.3,206.2,27.19,0.0,14.38};

static float revgain[7]={0.07737,0.1459,0.1101,0.09372,0.013486,0.0055159,0.006

3857};

VII.2.2 structures

TMhdr_t

structure to handle information contained in the TM header.

Contains:

short year,month,day,sat: the year,month,day and satellite

of acquisition

float gain[7],bias[7], gain and biases for conversion of

digital count to radiance for each of the seven bands

int usgsproj,usgszone,nbpixs,nblines: USGS projection

identifier, USGS zone identifier, number of pixels per line,

number of lines.

double usgsparams[15]: USGS projection parameters.

double ul_lon,ul_lat,ur_lon,ur_lat,ll_lon,ll_lat,lr_lon,lr_lat:

coordinates in latitude and longitude of the four corner of the

image (ul:Upper Left corner, ll: Lower Left corner, lr: Lower

Right corner, UR: Upper Right corner)

float ul_easting,ul_northing,ur_easting,ur_northing;

float ll_easting,ll_northing,lr_easting,lr_northing;

coordinate in Easting and Northing of the four corner of the

image.(ul:Upper Left corner, ll: Lower Left corner, lr: Lower

Right corner, UR: Upper Right corner)

float ts,fs: Solar zenith (ts) and azimuth (fs) angle

at the center of the scene.

double center_lon,center_lat: Latitude and longitude

of the center of the scene.

float center_easting,center_northing: Easting and Northing

of the center of the scene.

int center_x,center_y,center_hoffs: Pixels coordinate of

the center of the scene, offset of the center of the scene

to the center of the WRS zone (in pixels).

VII.2.2 global variables

GRID CELLS DARK TARGET STATISTICS

float avg_b1[NB_GPS]: table (the index being the grid cell)

containing the average count observed in TM band 1 for dark

targets pixels.

float avg_b2[NB_GPS]: table (the index being the grid cell)

containing the average count observed in TM band 2 for dark

targets pixels.

float avg_b3[NB_GPS]: table (the index being the grid cell)

containing the average count observed in TM band 3 for dark

targets pixels.

float avg_b7[NB_GPS]: table (the index being the grid cell)

containing the average count observed in TM band 7 for dark

targets pixels.

float std_b1[NB_GPS]: table (the index being the grid cell)

containing the standard deviation count observed in TM band 1 for dark

targets pixels.

float std_b2[NB_GPS]: table (the index being the grid cell)

containing the standard deviation count observed in TM band 2 for dark

targets pixels.

float std_b3[NB_GPS]: table (the index being the grid cell)

containing the standard deviation count observed in TM band 3 for dark

targets pixels.

float std_b7[NB_GPS]: table (the index being the grid cell)

containing the standard deviation count observed in TM band 7 for dark

targets pixels.

double sum_b1[NB_GPS]: table (the index being the grid cell)

containing the sum count observed in TM band 1 for dark

targets pixels.

double sum_b2[NB_GPS]: table (the index being the grid cell)

containing the sum count observed in TM band 2 for dark

targets pixels.

double sum_b3[NB_GPS]: table (the index being the grid cell)

containing the sum count observed in TM band 3 for dark

targets pixels.

double sum_b7[NB_GPS]: table (the index being the grid cell)

containing the sum count observed in TM band 7 for dark

targets pixels.

double sum2_b1[NB_GPS]: table (the index being the grid cell) containing the sum of square count observed in TM band 1 for dark targets pixels.

double sum2_b2[NB_GPS]: table (the index being the grid cell)

containing the sum of square count observed in TM band 2 for dark

targets pixels.

double sum2_b3[NB_GPS]: table (the index being the grid cell)

containing the sum of square count observed in TM band 3 for dark

targets pixels.

double sum2_b7[NB_GPS]: table (the index being the grid cell)

containing the sum of square count observed in TM band 7 for dark

targets pixels.

int samples[NB_GPS]: table (the index being the grid cell)

containing the number of dark target pixels.

int blockheight: the number of line of a cell (all cells

have the same dimensions).

int blockwidth: the number of pixel in a cell line (all cells

have the same dimensions).

int nbofblock: the number of buffers that should be read in the first

step of the processing (cloud screening, aerosol inversion). This

number of buffer is computed by dividing the total number of line in a

image by the size of the buffer given by QCBUF_SIZE a constant described

earlier.

VII.3 main program

!Description: ingest fast format data and a file containing

processing parameters. In regular mode it operates in two steps:

(1) reads band 1,2,3,4,6,7 and generate a Quality control Mask (see

cloud mask in TM_utils.c) as well as an aerosol optical

depth grid corresponding to the scene;

(2) reads band 1,2,3,4,5,7 the QC mask and performs atmospheric

correction for each non cloudy pixels, the correction include

correcting for adjacency effects.

The processing runs in three differents mode:

TOA mode: generate top of the atmosphere reflectance (no atmospheric

correction)

AC1 mode: apply atmospheric correction but no adjacency effect

correction

AC2: apply in addition to regular atmospheric correction, adjacency

effect correction.

!Input parameters:

int *argc: the number of parameter on the command line.

char **argv: table containing the command line parameter

argv[0]: contains the command name (process_TM-V1)

argv[1]: input prefix to generate name of the seven

TM bands used in the processing. File name for band "$x"

is $prefix.band$x

argv[2]: contains name of the processing parameter

file.

argv[3]: is reserved for options, that

could be: -notatmos (no atmospheric correction)

-noadja (no adjacency correction)

argv[4]: is reserved for options, that

could be: -notatmos (no atmospheric correction)

-noadja (no adjacency correction)

!Output Parameters:

int main: status -1 (process aborted)

0 (OK)

VII.4 C routines

VII.4.1 void compsums(int k,int tmpint, u_char *qcbuf[],int block)

!Description: routine used to computed summ and summ of squares of the digital

counts of channel 1,2,3,7 that are "darks" in channel 7 (0.015 < ref7 < 0.05 an

d not water)

!Input parameters:

int tmpint: the pointer in input buffer of the current pixel

int k: the grid point reference for that pixels

u_char *qcbuf[]: cloud/shadow mask buffer

int block: block index for circular buffer cb1,cb2,cb3,cb7

!Global variable

INPUTS

slope[k][0]: conversion slope from count to reflectance for band 1 at gridcell k

offset[0]: conversion offset from count to reflectance for band 1

slope[k][1]: conversion slope from count to reflectance for band 2 at gridcell k

offset[1]: conversion offset from count to reflectance for band 2

cb1[block][tmpint]: digital count in channel 1

cb2[block][tmpint]: digital count in channel 2

cb3[block][tmpint]: digital count in channel 3

cb7[block][tmpint]: digital count in channel 4

OUPUTS

each is for gridcell "k"

sum_b1[k]: summ of digital count in channel 1 "dark target" pixels

sum2_b1[k] summ of square digital count in channel 1 "dark target pixels

sum_b2[k]: summ of digital count in channel 2 "dark target" pixels

sum2_b2[k] summ of square digital count in channel 2 "dark target" pixels

sum_b3[k]: summ of digital count in channel 3 "dark target" pixels

sum2_b3[k] summ of square digital count in channel 3 "dark target" pixels

sum_b7[k]: summ of digital count in channel 7 "dark target" pixels

sum2_b7[k] summ of square digital count in channel 7 "dark target" pixels

samples[k]: number of "dark" target for gridcell "k"

!Output Parameters:

NONE

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 15:46:05 EST 1995

VII.4.2 void fillgaot(float aot[],float sigmaaot[])

!Description: fill in missing values in the aerosol grid band 1,2,3 based

on existing values (spatial interpolation). Missing values have been previously

set to -1. A filling can be distincted from a original value by looking at the

standart deviation of the optical depth which is set to -1 for a filling.

!Input parameters:

float aot[]: table containing average optical depths in band 1,2,3 for each gridcell

float sigmaaot[]: same as aot but contains

!Output Parameters:

float aot[]: table containing average optical depths in band 1,2,3 for each gridcell

float sigmaaot[]: same as aot but contains

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 15:46:05 EST 1995

VII.4.3 void interaot(float aot[])

!Description: interpolate aerosol thickness from band 1,2,3 to band 4,5,7.

as a function of wavelength (spectral interpolation).

!Input parameters:

float aot[]: table containing average optical depths in band 1,2,3 for each gridcell

!Output Parameters:

float aot[]: table containing average optical depths in band

1,2,3,4,5,7 for each gridcell

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 15:46:05 EST 1995

!Design Notes: For now a aerosol model similar to the continental

is assumed for spectral interpolation (Lambda^-1) in the future we

expect to use 1,2,3 to determine the aerosol model and its spectral

dependence.

VII.4.4 void atm_cor(u_char *buf,u_char *qc,int ib,int bufsize,int firstline,

params_t *inpars)

!Description: perform atmospheric correction ACP on a buffer of pixels.

Atmospheric correction consists in correcting atmospheric intrinsic

reflectance, transmission and multiple reflection effects. The

parameters for correction are interpolated at each pixel location from

a set of grid points for which the correction parameters are exactly

computed.

!Input parameters:

u_char *buf: buffer containing the pixels to be corrected (digital

counts)

u_char *qc: buffer containing cloud/shadow mask informations

int ib: band number to be corrected.

int bufsize: number of line in inputs buffer

int firstline: line number in the whole scene of the first line

of the input buffer.

params_t *inpars: structure containing processing parameters

!Output Parameters:

u_char *buf: buffer containing the corrected pixels (reflectances

scaled between rmin and rmax for that particular band)

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 15:46:05 EST 1995

VII.4.5 int prepwb(u_char *aebuf1,u_char *aebuf2,u_char *aebuf3,u_char *wbuff,int bufsize,params_t *inpars)

!Description: fill in the value into the working buffer used for adjacency

effect correction. The value in the working buffer consists in the

buffer being corrected and a number of line before and after the current

buffer that are needed for adjacency effect correction.

!Input parameters:

u_char *aebuf1: buffer containing values before the current one

u_char *aebuf2: buffer containing the current values (to be corrected)

u_char *aebuf3: buffer containing values after the current one

int bufsize: number of lines of buffer 1,2,3.

params_t *inpars: structure containing the number of pixels per

line of input buffer.

!CONSTANT USED

ADJACENCY EFFECT CORRECTION CONSTANTS

PSPMAT_RAD: radius of the correction in pixels.

MPSP_NBELT: number of elts of the correction mattrix (2*PSPMAT_RAD+1)*(2*PSPMAT_RAD+1)

TWOPSMAT_RAD: twice the correction radius

PSPMAT_NLC: number of line/collumn of the correction mattrix.

!Output Parameters:

u_char *wbuff: buffer for adjacency effect correction (containing

edges of the area to be corrected)

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 16:51:17 EST 1995

VII.4.6 int ae_cor(u_char *aebuf1,u_char *aebuf2,u_char *aebuf3,int ib, int bufsize,int AdjaCorrect,params_t *inpars)

!Description: performs correction for adjacency effect.

!Input parameters:

u_char *aebuf1: buffer containing values before the current one

u_char *aebuf2: buffer containing the current values (to be corrected)

u_char *aebuf3: buffer containing values after the current one

int bufsize: number of lines of buffer 1,2,3.

int AdjaCorrect: flag for adjacency effects correction (0 don't

apply, 1 apply)

params_t *inpars: structure containing the number of pixels per

line of input buffer.

!Output Parameters:

u_char *aebuf2: buffer containing the corrected values

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 17:00:24 EST 1995

VII.4.7 float getpressure(u_char *dem,double dlat,double dlon)

!Description: supply default value for pressure value in case the value

is missing from the processing parameters. The value is computed from

a Digital Elevation Model (ETOPO5) asuming a sea-level pressure of one

atmosphere.

!Input parameters:

u_char *dem: table containing altitude on 1degx1deg basis

double dlat: latitude of the location

double dlon: longitude of the location

!Output Parameters:

float getpressure: pressure at the location in millibar

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 17:00:24 EST 1995

VII.4.8 float getwv(u_char *wvcontent,double dlat,double dlon,short month)

!Description: supply default value for water vapor content, in case the value

is missing from the processing parameters. The value is read from the

OORT climatology for water vapor.

!Input parameters:

u_char *wvcontent: table containing water vapor climatogie 1degx1deg basis,

on each month.

double dlat: latitude of the location

double dlon: longitude of the location

short month: month of the Year (1=january,12=december)

!Output Parameters:

float getwv: integrated water vapor content in g/cm2.

!Revision History:

$LOG: process_TM.c,v1.0$

Original version: Thu Dec 7 17:00:24 EST 1995

VII.4.9 int initializepars(char *file,params_t *inpars)

!Description: decode the data in the parameter files and initialize

associated variables.

!Input parameters:

char *file: the name of the parameter file.

!Output Parameters:

params_t *inpars inpars: a structure of type params_t containing the

parameters.the structure params_t is described in details in file

TM_struc.h

int initializepars: status -1 (problem with opening parameter file)

0 (OK)

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

!Design Notes: The decoding process use a parser for decoding the parameter file.

Default parameters values are set and parameters present in the parameter file

are updated. If the parameter file can not be opened for read the routine returns

the status -1 if it returns 0.

VII.4.10 int rdparamsline(FILE *fd,int *par,int *ind,char *value)

!Description: parser to decode the parameter file, called by initializepars

function.

!Input parameters:

FILE *file: file handler pointing to the parameter file.

!Output Parameters:

int *par: the current parameter type decoded in the current line

of the parameter file. (from 1 to 13, 99 is the invalid type)

int *ind: a subtype of parameter for the special case (GP: grid point,

CAL: calibration coefficients ).

char *value: the actual parameter list

int rdparamsline: status 0 (OK),

-1 blank parameter line

1 unvalid parameter line

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.11 int nbofvalues(char *str)

!Description: count the number of parameters contained in str. A parameter

is separated by a other parameter using a comma (,).

!Input parameters:

char *str: string containing a list of parameters.

!Output Parameters:

int nbofvalues: the number of parameters

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.12 void print_pars(FILE *fd,params_t *inpars)

!Description: write to a file the processing parameter contained in inpars,

a structure of params_t type (for a description of params_t structure look into the

TM_struc.h).

!Input parameters:

FILE *fd: file handler to which to write the processing parameters.

params_t *inpars: variable containing the processing parameters.

!Output Parameters:

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.13 void process_error(int level,char *msg)

!Description: error messaging and handling routine, a message is always write

to the standart error output (stderr), and in case of fatal error an exit(-1) is

performed.

!Input parameters:

int level: the error level (ERRLEVEL_FATAL,ERRLEVEL_NORMAL,ERRLEVEL_WARNING).

char *msg: variable containing the error message.

!Output Parameters:

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.14 void logmsg(FILE *fdlog,char *msg)

!Description: write a message to a log file.

!Input parameters:

FILE *fdlog: File handler to the log file.

char *msg: message to be writen to log file.

!Output Parameters:

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.15 void cloud_mask(u_char *b3,u_char *b4,u_char *b6, float tsurf,float thrshld1,float thrshld2,float fk1,float fk2, u_char *qcbuf[], int imageline, int buf_size,params_t *inpars)

!Description: generate a cloud/shadow mask for a given buffer of input data.

!Input parameters:

u_char *b3: buffer containing the digital count observed in TM-channel3.

u_char *b4: buffer containing the digital count observed in TM-channel4.

u_char *b6: buffer containing the digital count observed in TM-channel6.

float tsurf: the assumed surface temperature (processing parameter)

float thrshld1: the maximum temperature of a cloudy pixel

float thrshld2: the minimum temperature of a clear pixel

float fk1: value used to convert radiance to temperature (Planck function)

float fk2: value used to convert radiance to temperature (Planck function)

int imageline: line number in the full image of the first line of the buffer.

int buf_size: number of line of the input buffer

params_t *inpars: structure containing the processing parameter (see TM_struc.h

for details)

!Global variable (input):

extern int blockheight: number of line of the input buffer

extern int blockwidth: number of pixels of the input buffer

extern float slope[7][NB_GPS]: table containing conversion coefficient (slope)

from count to reflectance (band 1,2,3,4,5,7) for each grid-cell.

extern float offset[7]: table containing offset for conversion count to reflectance.

(band 1,2,3,4,5,7) or count to radiance (band 6).

extern float slopecf[7]: table containing conversion coefficient (slope)

from count to radiance for band 1 to 7.

!Output Parameters:

u_char *qcbuf[]: the cloud/shadow mask also called QC buffer or quality

control buffer. it classifies the pixels into severals categories, CLOUDY,

CLEAR, MIXED, SHADOW, LAND, WATER. The CLEAR,MIXED,CLOUDY categories are

based on three test using threshold on the expected reflectance in channel 3,

the apparent temperature in channel 6, the ratio between channel 3 and 4.

The SHADOW mask is generated from the CLOUDY pixels using apparent temperature

of the pixels, expected surface temperature tsurf, and the sun-sensor geometry.

The LAND/WATER information is based on the value of the ratio between the

reflectance in channel 3 and 4 (CH4/CH3 < 0.9 -> WATER; CH4/CH3 >1.1 LAND)

The mask are as followed: bit 0 (2^0): not used here

bit 1 : not used here

bit 2 : not used here

bit 3 : 1 WATER, 0 LAND

bit 4-5 : 00 CLEAR, 01 MIXED, 10 CLOUDY, 11 SHADOW

BIT 6 : 0 VALID DATA , 1 INVALID DATA

BIT 7 : not used here

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.16 float comp_dsol(short jday)

!Description: compute the variation in solar constant due to variation of

sun-earth distance as a function of the number of day since the 1rst January

of the year of the acquisition date.

!Input parameters:

short jday:number of day since the 1rst January

of the year of the acquisition date.

!Output Parameters:

float comp_dsol: multiplicative correction coefficient for

the solar constant.

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

!References and Credits: the routine came originally from a FORTRAN subroutine

for the 6S code.

VII.4.17 short julian_day(short year,short month,short day)

!Description: compute the number of day between the 1rst January

of the year of the acquisition date and the acquisition date.

!Input parameters:

short year: year of acquisition

short month: month of acquisition between 1 and 12

short day: calendar day in the month

!Output Parameters:

short julian_day: the number of day between the 1rst January

of the year of the acquisition date and the acquisition date

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.18 int getGP(int i,int j,int bh,int bw)

!Description: function used to compute the closest grid point coordinate

to a given pixel (reference in line collumn within the image).

!Input parameters:

int i: line coordinate of the pixels in the image

int j: collumn coordinate of the pixels in the image

int bh: Grid cell number of line

grid cell: Grid cell number of collumn

!Output Parameters:

int getGP: grid point index.

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.19 void initcoefsH2O(float *aH2O,float *bH2O,float *cH2O)

!Description: initialize coefficient for water vapor transmission

functions (in TM bands 1,2,3,4,5,7, 6 must be ignored).

!Input parameters:

!Output Parameters:

float *aH2O: table containing coefficients "a" for water vapor transmission

float *bH2O: table containing coefficients "b" for water vapor transmission

float *cH2O: table containing coefficients "c" for water vapor transmission

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.20 void initcoefsO3(float *aO3)

!Description: initialize coefficient for ozone transmission

functions (in TM bands 1,2,3,4,5,7, 6 must be ignored).

!Input parameters:

!Output Parameters:

float *aO3: table containing ozone equivalent optical depth.

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.21 initcoefsCO2(float *aCO2,float *bCO2)

!Description: initialize coefficient for carbon dyoxide transmission

functions (in TM bands 1,2,3,4,5,7, 6 must be ignored).

!Input parameters:

!Output Parameters:

float *aCO2: table containing coefficients "a" for carbon dyoxide transmission

float *bCO2: table containing coefficients "b" for carbon dyoxide transmission

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

!References and Credits: coefficient were derived by running 6S for the TM bands,

original expression for transmission was proposed by D. Tanre in: Tanre, D., Holben, B. N. and Kaufman, Y. J. (1992), Atmospheric Correction algorithm for NOAA-AVHRR Products: Theory and Application, IEEE Trans. Geosci. Rem. Sens. 30(2):231-248.

VII.4.22 initcoefsO2(float *aO2,float *bO2)

!Description: initialize coefficient for oxygen transmission

functions (in TM bands 1,2,3,4,5,7, 6 must be ignored).

!Input parameters:

!Output Parameters:

float *aO2: table containing coefficients "a" for oxygen transmission

float *bO2: table containing coefficients "b" for oxygen transmission

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.23 void initcoefsNO2(float *aNO2,float *bNO2)

!Description: initialize coefficient for nitrous oxyde transmission

functions (in TM bands 1,2,3,4,5,7, 6 must be ignored).

!Input parameters:

!Output Parameters:

float *aNO2: table containing coefficients "a" for nitrous oxyde transmission

float *bNO2: table containing coefficients "b" for nitrous oxyde transmission

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.4.24 initcoefsCH4(float *aCH4,float *bCH4)

!Description: initialize coefficient for methame transmission

functions (in TM bands 1,2,3,4,5,7, 6 must be ignored).

!Input parameters:

!Output Parameters:

float *aCH4: table containing coefficients "a" for methane transmission

float *bCH4: table containing coefficients "b" for methane transmission

!Revision History:

$LOG: TM_utils.c,v1.0$

Original version: Sat Oct 28 11:54:01 EDT 1995

VII.5 fortran routines

VII.5.1 subroutine aeroso (iaer,co,xmud)

!Description: routine uses to initialize aerosol phase function, extinction and

!scattering coefficient at ten wavelengths.

!Input parameters:

! int iaer: the aerosol model identifier

! iaer=1 continental aerosol

! iaer=2 maritime aerosol

! iaer=3 urban

! iaer=4 user defined see co(4) in that case

! iaer=5 desert background model

! float co(4): composition of the aerosol in case of user defined

! aerosol type (iaer=4). co(1) is the fraction of dust like particle

! co(2) the fraction of water soluble particle, co(3) the fraction of

! oceanic type particles, co(4) is the fraction of soot particles.

! (co(1)+co(2)+co(3)+co(4)=1)

! float xmud: is the cosine of the scattering angle

!common variables:

! INPUT

! common /sixs_sos/phasel(10,83),cgaus(83)

! real cgaus(83): the cosine of the 83 angle values for which the

! phase function is computed. only used to compute phase(10) by

! interpolation.

! OUTPUT

! common /sixs_aer/ ext(10),ome(10),gasym(10),phase(10)

! real ext(10): the extinction coefficient at ten wavelenghts

! real ome(10): the single scattering albedo at ten wavelenghts

! real gasym(10): the assymetry parameters at ten wavelenghts

! real phase(10): the value of the phase function at ten wavelengths

! for scattering angle of cosine xmud

! common /sixs_sos/phasel(10,83),cgaus(83)

! real phasel(10,83): the aerosol phase function at 10 wavelenghts

! and 83 angles.

! INPUT

! common /sixs_sos/phasel(10,83),cgaus(83)

! real cgaus(83): the cosine of the 83 angle values for which the

! phase function is computed.

!Output Parameters:

!Revision History:

!$LOG: AEROSO.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.2 subroutine bdm

!Description: routine uses to initialize aerosol phase function at ten

!wavelengths in the case of desert background aerosol (called by AEROSO.f)

!Input parameters:

!common variables:

! OUTPUT

! common /sixs_aerbas/ ph(10,83)

! real ph(10,83): the aerosol phase function at 10 wavelenghts

! and 83 angles.

!Output Parameters:

!Revision History:

!$LOG: BDM.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.3 subroutine caluoz(uo3,jour,xlat)

!Description: routine used to interpolate the ozone climatological content

!from the London et al. climatology.

!Input parameters:

! short jour: the julian day (between 1 and 366)

! float xlat: the latitude in deg (between -90 and 90)

!Output Parameters:

! float uo3: the amount of ozone at this day of year and

! location in unit of cm.atm

!Revision History:

!$LOG: CALUOZ.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

VII.5.4 subroutine chand (xphi,xmuv,xmus,xtau,xrray)

!Description: routine used to compute the intrinsic reflectance of a pure

! molecular atmosphere.

!Input parameters:

! float xphi: relative azimuth angle: difference between the sun azimuth

! and the sensor azimuth (in degree)

! float xmuv: cosine of the view zenith angle (unitless).

! float xmus: cosine of the sun zenith angle (unitless).

! float xtau: optical depth of the molecules (unitless).

!Output Parameters:

! float xrray: intrinsic atmospheric reflectance in case of

! pure molecular scattering.

!Revision History:

!$LOG: CHAND.f,v1.0$

!Design Notes: this is a 6S routine.

VII.5.5 subroutine compgrac(ttaer,tts,ttv,tphi,uwv,uoz,pres,ib,ttgog,ttgh2o,

tttd,tttu,trray,troatm,ttdd,ttdu,tsatm)

!Description: routine used to compute atmospheric transmission and reflectance

! parameters for the grid points (16) of the TM image. This routine interfaces the

!main C processing code with the 6S routines. This routine is called for each

!TM bands at the exception of the thermal band.

!Input parameters:

! float ttaer: table containing the aerosol optical depth at wavelenght

! corresponding to band number ib for each grid point (unitless)

! float tts: table containing the solar zenith angle for each grid point

! (degree)

! float ttv: table containing the view zenith angle for each grid point

! (degree)

! float tphi: table containing the relative azimuth angle

! (difference between the sun azimuth angle and the view azimuth

! angle) for each grid point (degree)

! float uwv: value of the integrated water vapor content that is

! constant for the whole scene (g/cm2).

! float uoz: value of the integrated ozone content that is

! constant for the whole scene (cm.atm).

! float pres: value of the surface pressure that is

! constant for the whole scene (Millibars).

! int ib: TM band number ib=(1,2,3,4,5,6) for bands (1,2,3,4,5,7)

!Output Parameters:

! float ttog: table containing for each grid point the gaseous

! transmission for all gases except water vapor (unitless).

! float ttgh2o: table containing for each grid point the water vapor

! gaseous transmission water vapor (unitless).

! float tttd: table containing for each grid point the total downward

! scattering transmission (molecules+aerosols) (unitless)

! float tttu: table containing for each grid point the total upward

! scattering transmission (molecules+aerosols) (unitless)

! float ttray: table containing for each grid point the intrinsic

! reflectance of the molecules only (unitless)

! float ttroatm: table containing for each grid point the intrinsic

! atmospheric reflectance (molecules+aerosols) (unitless)

! float ttdd: table containing for each grid point the diffuse downward

! scattering transmission (molecules+aerosols) (unitless)

! float ttdu: table containing for each grid point the diffuse upward

! scattering transmission (molecules+aerosols) (unitless)

! float ttdu: table containing for each grid point the atmosphere

! spherical albedo (molecules+aerosols) (unitless)

!Revision History:

!$LOG: COMPGRAC.f.f,v1.0$

!Original version: Wed Dec 13 09:33:38 EST 1995

VII.5.6 subroutine compgrao(tmref,tsref,smref,ssref,ts,tv,phi,uwv,uoz,pres,ib, tmaer,tsaer)

!Description: routine used to compute aerosol optical depth for band ib

! It uses as input the mean and standart deviation of the reflectance measured

! in channel ib (tmref,smref) for dark target pixels,

! and the predicted signal of the surface in that band (smfref,ssref) based

! on channel 7 measured reflectance.

!Input parameters:

! float tmref: the measured mean reflectance of

! dark target pixels in channel ib for each gridell.

! float tsref: table containing the standart deviation of from

! mean reflectance (tmref) of dark target pixels in channel ib

! for each gridell.

! float smref: table containing the measured mean surface reflectance of

! dark target pixels in channel ib for each gridell (based on channel 7)

! float ssref: table containing the measured mean surface reflectance of

! dark target pixels in channel ib for each gridell (based on channel 7)

! float ts: the sun zenith angle corresponding to that observation

! (degree)

! float tv: the view zenith angle corresponding to that observation

! (degree)

! float phi: the relative azimuth angle

! (difference between the sun azimuth angle and the view azimuth

! angle) corresponding to that observation (degree)

! float uwv: value of the integrated water vapor content

! corresponding to that observation (g/cm2).

! float uoz: value of the integrated ozone content that is

! corresponding to that observation (cm.atm).

! float pres: value of the surface pressure

! corresponding to that observation (Millibars).

! int ib: TM band number ib=(1,2,3,4,5,6) for TM bands (1,2,3,4,5,7)

! NB: this modules is only called with ib=1,2,3 now.

!Output Parameters:

! float tmaer: the retrieved mean aerosol optical depth

! at wavelength corresponding to TM band ib central wavelenght.

! float tsaer: the retrieved standart deviation of aerosol optical depth

! at wavelength corresponding to TM band ib central wavelenght.

! This computation is based on values of tsref and ssref.

!Revision History:

!$LOG: COMPGRAO.f,v1.0$

!Original version: Wed Dec 13 09:33:38 EST 1995

VII.5.7 subroutine csalbr(xtau,xalb)

!Description: routine used to compute the spherical albedo of a pure

!molecular atmosphere of optical depth xtau

!Input parameters:

! float xtau: the optical depth of the pure molecular

! atmosphere (unitless).

!Output Parameters:

! float xalb: the spherical albedo (unitless).

!Revision History:

!$LOG: CSLABR.f,v1.0$

!Original version: Wed Dec 13 09:33:38 EST 1995

!Design Notes: This is a 6S routine.

VII.5.8 real function fintexp3(xtau)

!Description: this function returns the exponential integral of the third

! order for value xtau. It is used by CSALBR routine.

!Input parameters:

! float xtau: argument (unitless).

!Output Parameters:

! float fintexp3: The exponential integral of order 3.

!Revision History:

!$LOG: CSLABR.f,v1.0$

!Original version: Wed Dec 13 09:33:38 EST 1995

Design Notes: This is a 6S routine.

VII.5.9 real function fintexp1(xtau)

!Description: this function returns the exponential integral of the third

! order for value xtau. It is used by CSALBR routine.

!Input parameters:

! float xtau: argument (unitless).

!Output Parameters:

! float fintexp1: The exponential integral of order 1 for xtau.

!Revision History:

!$LOG: CSLABR.f,v1.0$

!Original version: Wed Dec 13 09:33:38 EST 1995

Design Notes: This is a 6S routine.

VII.5.10 subroutine discre(ta,ha,tr,hr,it,nt,yy,dd,ppp2,ppp1,zx)

!Description: routine used to compute the altitude of the bottom of

!the atmospheric layer number=it based on the relative proportion

!of rayleigh and molecules in the atmosphere. The routine returns a altitude

! that minimize makes each layer about the same thickness.

!Input parameters:

! float ta: the aerosol optical depth (unitless).

! float ha: the aerosol scale height (km)

! float tr: the molecules optical depth (unitless)

! float hr: the molecules scale height (km)=8km

! int it: the number of the current layer (0=Top ,nt=bottom)

! int nt: the total number of layer.

! float yy: the total optical depth from layer 0 to it-1.

! float dd: the optical depth of layer (it-1)

! float ppp2: altitude of the top of the atmosphere (km)=300km

! float ppp1: altitude of the bottom of the atmosphere (km)=0km

!Output Parameters:

! float zx: the altitude of the top of layer number=it.

!Revision History:

!$LOG: DISCRE.f,v1.0$

!Original version: Wed Dec 13 09:33:38 EST 1995

!Design Notes: This is a 6S routine.

VII.5.11 subroutine dust

!Description: routine uses to initialize aerosol phase function at ten

!wavelengths in the case of dust-like particles (called by AEROSO.f)

!Input parameters:

!common variables:

! OUTPUT

! common /sixs_aerbas/ ph(10,83)

! real ph(10,83): the aerosol phase function at 10 wavelenghts

! and 83 angles.

!Output Parameters:

!Revision History:

!$LOG: DUST.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: This is a 6S routine.

VII.5.12 subroutine gauss(x1,x2,x,w,n)

!Description: this subroutine generates a gauss quadrature (angle+weight)

!for a give range (x1,x2) and a chosen number of angle. This quadrature

!is used in the successive order of scattering code (ISO.f,OS.f) to compute

!signals at those specific angles,and performs integration. This routine

! is called by COMPGRAO.f and COMPGRAC.f.

!Input parameters:

! float x1: minimum angle of the quadrature (radian)

! float x2: maximum angle of the quadrature (radian)

! int n: desired number of angle in the quadrature

!Output Parameters:

! float x(n): table containing the cosine of the n gauss angles (unitless)

! float w(n): table containing the weighting factor to be

! used for integration.

!Revision History:

!$LOG: GAUSS.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.13 subroutine iso(tamoy,trmoy,pizmoy,tamoyp,trmoyp,palt,nt,mu,rm,gb,xf)

!Description: this subroutine computes the transmission and spherical

!albedo of an atmosphere consisting of a mixing of molecules and aerosols.

!Input parameters:

! float tamoy: aerosol optical depth (unitless)

! float trmoy: molecules optical depth (unitless)

! float pizmoy: aerosol single scattering albedo

! float tamoyp: aerosol optical depth of the portion of the

! atmosphere under the observer (unitless)

! float trmoyp: molecules optical depth of the portion of the

! atmosphere under the observer (unitless)

! float palt: altitude of the observer (km)

! int nt: number of atmosphere layers to be used in the computation.

! int mu: indicator of the number of angles to be used

! in the computation (number of gauss angles= 2*mu-1)

! float rm(-mu:mu): table containing the cosine of the angles to be

! used in the computation.

! rm(-mu+1:-1),rm(1:mu-1) are the gauss angles

! rm(mu): contains the cosine of the view angle.

! rm(-mu): contains the -cosine of the view angle.

! rm(0): contains the cosine of the view angle.

! float gb(-mu:mu): table containing the weight to be

! used in the integrals computation. Only gb(-mu+1:mu-1) with

! the exception of gb(0) are significant.

!Output Parameters:

! float xf(-1:1): table containing the result of the computation.

! xf(-1): partial atmospheric transmission upward till observer level

! xf(0): atmosphere spherical albedo

! xf(1):total atmospheric transmission upward

!NB this routine can be used to compute downward transmission by using

!the reciprocity principle and inverting view and sun zenith angle.

!Revision History:

!$LOG: ISO.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.14 subroutine kernel(is,mu,rm,xpl,psl,bp)

!Description: This routine compute the legendre polynomials (xpl, at order is, for

! the needed angle rm(mu). That are used by OS.f and ISO.f to compute

! the radiation field at each atmospheric level.

!Input parameters:

! int is: the s parameter in legendre series.

! int mu: number of angles in table rm.

! float rm(-mu:mu): table containing the cosine of the angle

! for which legendre polynomial need to be evaluated.

!Output Parameters:

! float psl(-1:81,-mu:mu): table containing the legendre polynomial

! of order -1 to 81, for rm(-mu:mu) angle values.

! float xpl(-mu:mu): table containing the the legendre polynomial

! of order=is (if is=0,1,2), for rm(-mu:mu) angles values.

! float spl(0,-mu:mu): table containing the product

! Betal(is)Psl(is,mu)Psl(is,mu') where mu' is the incidence angle

! and the coefficient of the fourier decomposition (order is) of

! the phase function.

!Revision History:

!$LOG: KERNEL.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.15 subroutine ocea

!Description: routine uses to initialize aerosol phase function at ten

!wavelengths in the case of oceanic particles (called by AEROSO.f)

!Input parameters:

!common variables:

! OUTPUT

! common /sixs_aerbas/ ph(10,83)

! real ph(10,83): the aerosol phase function at 10 wavelenghts

! and 83 angles.

!Output Parameters:

!Revision History:

!$LOG: OCEA.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.16 subroutine os(tamoy,trmoy,pizmoy,tamoyp,trmoyp,palt,

phirad,nt,mu,np,rm,gb,rp,xl)

!Description: this subroutine computes the intrinsic atmospheric reflectance

!of an atmosphere consisting of a mixing of molecules and aerosols.

!Input parameters:

! float tamoy: aerosol optical depth (unitless)

! float trmoy: molecules optical depth (unitless)

! float pizmoy: aerosol single scattering albedo

! float tamoyp: aerosol optical depth of the portion of the

! atmosphere under the observer (unitless)

! float trmoyp: molecules optical depth of the portion of the

! atmosphere under the observer (unitless)

! float palt: altitude of the observer (km)

! float phirad: relative azimuth of the observation (azimuth

! of the sun minus the azimuth of the observer) (radian)

! int nt: number of atmosphere layers to be used in the computation.

! int mu: indicator of the number of zenith angles to be used

! in the computation (number of gauss angles= 2*mu-1)

! int np: indicator of the number of azimuth angles to be used

! in the computation (number of gauss angles= 2*np-1)

! float rm(-mu:mu): table containing the cosine of the zenith angles to be

! used in the computation.

! rm(-mu+1:-1),rm(1:mu-1) are the gauss angles

! rm(mu): contains the cosine of the view zenith angle.

! rm(-mu): contains the -cosine of the view zenith angle.

! rm(0): contains the cosine of the sun zenith angle.

! float gb(-mu:mu): table containing the weight to be

! used in the integrals computation. Only gb(-mu+1:mu-1) with

! the exception of gb(0) are significant.

! float rp(np): table containing the cosine of the azimuth angles to be

! used in the computation.

!Output Parameters:

! float xl(-mu:mu,np): table containing the result of the computation.

! that is:

! xl(-mu+1:-1,np): bottom of the atmosphere (downward) radiance field.

! xl(1:mu-1,np): top of the atmosphere (upward) radiance field.

! xl(-mu,1): apparent atmosphere intrinsic reflectance for the given

! geometrical conditions (phirad,xmus,xmuv) at the observer altitude.

! xl(mu,1): atmosphere intrinsic radiance for the given

! geometrical conditions (phirad,xmus,xmuv)

!NB The results in xl(-mu:mu,np) are given in radiance. It is in this case

!the normalized radiance (incident flux=pi). To convert those radiances

!to reflectances just divide by the cosine of the sun zenith angle.

!Revision History:

!$LOG: OS.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.17 subroutine print_error(tex)

!Description: simple error handler routine used by the 6S routines. Simply

! send message to the iwr logical unit.

!Input parameters:

! char *tex: error message to be printed.

!Output Parameters:

!NONE

!Revision History:

!$LOG: PRINT_ERROR.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.18 subroutine psprdf(pres,ttu,tdu,ib,dl,mp,alpha,nt)

!Description: This routine computes the atmospheric point spread mattrix

!(-nt:nt:-nt:nt) pixels to be used for adjacency effect correction. The computation

! is done using the average value of aerosol optical depth over the whole

! TM scene.

!Input parameters:

! float pres: the surface pressure (millibar)

! float ttu: the average total upward transmission (unitless)

! float tdu: the average diffuse upward transmission (unitless)

! int ib: the TM band number ib=(1,2,3,4,5,6)

! for TM band=(1,2,3,4,5,7)

! float dl: the size of a pixel in meters.

! int nt: sizing parameters for the point spread function

! mattrix in number of pixels corresponding to the point spread

! function radius. The size of the mattrix is (2*nt+1)*(2*nt+1)

!Output Parameters:

! float alpha: ratio of total to diffuse transmission to be

! used in the adjacency correction parts.

! float mp(-nt:nt,-nt:nt): point spread function mattrix.

!Revision History:

!$LOG: PSPRDFUNC.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

VII.5.19 subroutine soot

!Description: routine uses to initialize aerosol phase function at ten

!wavelengths in the case of soot particles (called by AEROSO.f)

!Input parameters:

!common variables:

! OUTPUT

! common /sixs_aerbas/ ph(10,83)

! real ph(10,83): the aerosol phase function at 10 wavelenghts

! and 83 angles.

!Output Parameters:

!Revision History:

!$LOG: SOOT.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.20 subroutine trunca(coeff)

!Description: This routine performs the decomposition in fourier series

! of the aerosol phase function. It returns through the common variable

! zone the coefficient of decomposition Betal(is), and as an output argument

! the truncation coefficient that represents the difference between the

! original phase function and its fourier transform.

!global variable

! INPUT

! common /sixs_trunc/pha(1:83),betal(0:80)

! float pha(1:83): table containing the aerosol phase function at

! 83 angles.

! OUPUT

! common /sixs_trunc/pha(1:83),betal(0:80)

! float betal(1:83): table containing the coefficient of fourier

! decomposition till order 81.

!Input parameters:

!NONE

!Output Parameters:

! float xcoef: truncation coefficient equal to the integral

! over two pi of the phase function.

!Revision History:

!$LOG: TRUNCA.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VII.5.21 subroutine wate

!Description: routine uses to initialize aerosol phase function at ten

!wavelengths in the case of water soluble particles (called by AEROSO.f)

!Input parameters:

!common variables:

! OUTPUT

! common /sixs_aerbas/ ph(10,83)

! real ph(10,83): the aerosol phase function at 10 wavelenghts

! and 83 angles.

!Output Parameters:

!Revision History:

!$LOG: WATE.f,v1.0$

!Original version: Tue Dec 12 15:56:22 EST 1995

!Design Notes: this routine is from the 6S code.

VIII Glossary

C Ci Condition to determine presence of cloud in a pixel (i=1,2,3)

C'i Condition to determine absence of cloud in a pixel (i=1,2,3)

d2 sun-earth distance variation factor [unitless]

Es Solar constant (W.m-2.str-1.mm-1)

f azimuth angle with respect to the north

[degree]

fs azimuth solar angle

fv azimuth view angle

l wavelength [mm]

m air mass [unitless]

m cosine of zenithal angle [unitless]

ms solar zenith angle

mv view zenith angle

n Angstrom exponent determines the spectral variation of optical depth

(l-n)

P surface pressure [mb]

P0 default surface pressure at sea-level (1013 mb)

q zenithal angle with respect of local normal (90°-elevation)

[degree]

qs solar zenith angle

qv view zenith angle

r Reflectance [unitless]

r(TMi) Top of the atmosphere reflectance in TM band i (i,1..6,7)

rR+A Reflectance due to molecules (Rayleigh) scattering and aerosols.

rTOA Top of the atmosphere reflectance (after calibration)

rint Intermediate reflectance (before multiple reflection correction)

Intermediate reflectance (before adjacency effect corrections)

<rs> correction term for adjacency effect

rs Final reflectance product

S Spherical albedo [unitless] controls multiple reflection process.

SR Spherical albedo due to molecules (Rayleigh)

SR+A Spherical albedo due to aerosol+molecules

T Temperature [Kelvin]

T(TM6) Temperature from TM band 6

TA Atmospheric transmission for aerosols

TA(mv) Upward transmission for aerosols

TA(ms) Downward transmission for aerosols

TR Atmospheric transmission for molecules (Rayleigh)

TR(mv) Upward transmission for molecules

TR(ms) Downward transmission for molecules

TR+A Atmospheric transmission for molecules and aerosols.

TR+A(mv) Upward transmission for molecules and aerosols.

TR+A(ms) Downward transmission for molecules and aerosols.

t Optical depth [unitless]

tA Aerosol optical depth

tR Molecules (Rayleigh) optical depth