**2. Methodology**

The innovation of this methodology is the retrieval of AOT by using an image-based algorithm based on the radiative transfer equation to derive AOT using the reflectance values from the darkest pixel method of atmospheric correction. In the image-based integrated method, the RT equation retrieves AOT values from the satellite images with as much precision as possible. Also, the integrated use of field spectroscopy, GIS, and remote sensing for AOT retrieval is an important contribution to the field of remote sensing.

The first step of the methodology is the preprocessing of the satellite images, including geometric and radiometric correction. Also, the satellite images need to be atmospherically corrected. The darkest pixel atmospheric correction method can be applied, as it provides reasonable atmospheric correction, especially in cloud-free skies [28, 38, 39]. If there is a known nonvariant target or overpass *in-situ* measurements available in the selected site where the true ground reflectance values of the darkest target are known, then the modified DP method can be applied to result in more accurate reflectance values to the corrected satellite image [40]. The images are then processed using the RT equations in order to retrieve the AOT levels. Specific parameters from the satellite images, sun photometer measurements, and the reflectance values from the atmospherically corrected images are incorporated. In this methodology, the AOT values are calculated at the 500 nm wavelength, so Landsat band 1 will be used to derive AOT levels through the algorithm. Also, GIS analysis is done to produce thematic maps showing the AOT levels in the area of interest. In this chapter, an accuracy assessment is performed in the example of Limassol, Cyprus, where the AOT values retrieved from the GIS analysis are compared with the *in-situ* sunphotometer measurements from the sunphotometers. **Figure 1** shows the methodology for the AOT retrieval.

characterizations have included using *in-situ* measurements, ground-based remote sensing, and satellite observations [1–3]. AOT derived from satellite images often requires further validation [8]. The accuracy of satellite-derived AOT is frequently assessed by comparing satellite-based AOT with **AE**rosol **RO**botic **NET**work (AERONET, a program established by

Research indicates that the determined AOT from satellite image data can be used as a tool to assess air pollution as well as identify the sources of local emissions [2, 7, 11–24]. AOT values can be used as a way of measuring air quality; determining AOT in large-scale pollution areas provides a synoptic, cost-effective means to further assess the air quality in such areas [16, 18– 24]. Indeed, the AOT values derived from the atmospheric path radiance can be utilized to assess and monitor air quality and atmospheric pollution [8]. The image-based integrated method presented in this chapter can accurately calculate AOT values retrieved from satellite imagery by using radiative transfer (RT) equations and GIS. The method can be used to visually display AOT levels using thematic maps in order to identify concentrations of AOT over an urban area [11]. The image-based method can also be used with archived satellite images, thereby providing detailed information regarding spatial aerosol concentration overtime [11]. AOT is directly related to the atmospheric aerosol load, which is the main variable describing the effects of aerosols on radiative transfer in the Earth's atmosphere. According to Guanter et al. [25], modeling atmospheric constituents and surface reflectance involves modeling the radiative transfer across the atmosphere. The key parameter for assessing atmospheric pollution is the aerosol optical thickness, which is also the most important unknown of every atmospheric correction algorithm for solving the radiative transfer equation and removing atmospheric effects from remotely sensed satellite images [8, 26–29]. Several researchers [5, 18, 26, 28, 29–37] found that using the radiative transfer and atmospheric modeling in conjunction with field measurements of aerosol optical thickness can yield more accurate atmospheric

The image-based integrated method discussed in this chapter combines the RT equation, satellite imagery, the modified darkest pixel (DP) method of atmospheric correction, and GIS to derive AOT measures. An example including 11 Landsat satellite images with *insitu* measurements over a specific period of time is used to assess the AOT values based on the 30 × 30 m spatial resolution of Landsat over the city of Limassol, Cyprus and create

The innovation of this methodology is the retrieval of AOT by using an image-based algorithm based on the radiative transfer equation to derive AOT using the reflectance values from the darkest pixel method of atmospheric correction. In the image-based integrated method, the RT equation retrieves AOT values from the satellite images with as much precision as possible. Also, the integrated use of field spectroscopy, GIS, and remote sensing for AOT retrieval is an

**NASA**) or field-based sun photometers [9, 10].

142 Aerosols - Science and Case Studies

corrections instead of using simple image-based techniques.

thematic maps to display the AOT levels.

important contribution to the field of remote sensing.

**2. Methodology**

### **2.1. Detailed description of the image-based integrated method for determining AOT**

According to Hadjimitsis and Clayton [40], RT equations can be used to retrieve AOT values from the path radiance. Therefore, the algorithm used in the image-based integrated method is formulated to solve for AOT, which is included in the equation for path radiance (*Lp*) and Mie scattering [11]. To calculate AOT, Eq. (1) is used.

Path radiance is the sum of the atmospheric path radiance due to Rayleigh (*L*pr) and Mie scattering (*L*pa) [11, 29, 40–43]:

$$L\_{\rho} = L\_{\rm pa} + L\_{\rm pa} \tag{1}$$

where *Lp* is the path radiance (integrated band radiance measured in W/m2 sr), *L*pr is the atmospheric path radiance due to Rayleigh scattering in W/m2 sr [44], and *L*pa is the atmospheric path radiance due to Mie scattering in W/m2 sr [44].

#### *2.1.1. Calculating path radiance*

The path radiance "*Lp*" measured in Wm−2 μm−1 sr−1, is calculated according to Eq. (2) [42, 45]:

$$L\_{\rho} = L\_{\rm ts} - \frac{\pi \left(\mu\right) \uparrow \rho\_{\rm u} E\_G}{\pi},\tag{2}$$

where *Lts* is the measured at-satellite radiance, *τ*(*μ*) ↑ is the direct (upward) target sensor atmospheric transmittance, *ρ*tg is the ground target reflectance, and *EG* is the global irradiance reaching ground.

Eq. (3) can be used to calculate the at-satellite radiance (*L*ts) for Landsat satellite images that contain gain and calibration offset values in the satellite image's header file, as described in the Landsat-7 Handbook [46]:

$$L\_{\rm ts} = \mathbf{gain} \times \mathbf{DN} + \text{offset} \tag{3}$$

The "*EG*" is calculated according to Eq. (4) [45]:

$$E\_G = E\_{0\lambda} \mu exp \left[ -\left(\frac{1}{2\tau\_r} + \frac{1}{6\tau\_a}\right) / \mu \right] \tag{4}$$

where *E*0λ is the solar irradiance at the top of the atmosphere in Wm−2, found in the calibration file of the satellite image, *μ* is the cosine of the solar zenith angle (cos*θ*0), found on the satellite image header file, *τr* is the Rayleigh optical thickness, and *τa* is the AOT (which is unknown).

The "*τ*(*μ*) ↑" [42] is calculated in Eq. (5):

$$\text{tr}(\mu)\uparrow = \exp\left[-\left(\tau\_r + \tau\_a\right)/\mu\right] \tag{5}$$

Although several equations to calculate *τr* have been developed, most radiative transfer equations used the equation provided by Moller [47], as indicated in Eq. (6), which is considered to be the most appropriate [34]:

The Image-Based Integrated Method for Determining and Mapping Aerosol Optical Thickness http://dx.doi.org/10.5772/65279 145

$$
\pi\_r \left( \mathcal{J}\_\varepsilon \right) = 0.00879 \mathcal{J}\_\varepsilon^{-4.09} \tag{6}
$$

where *λC* is the central wavelength for each Landsat band.

#### *2.1.2. Calculating Rayleigh scattering (Lpr)*

*LLL <sup>p</sup>* = + pr pa (1)

= - (2)

ts *L* = ´+ gain DN offset (3)

) = é- + ù exp / ( *r a* ) ë û (5)

sr), *L*pr is the

(4)

sr [44], and *L*pa is the atmospheric

where *Lp* is the path radiance (integrated band radiance measured in W/m2

sr [44].

The path radiance "*Lp*" measured in Wm−2 μm−1 sr−1, is calculated according to Eq. (2) [42, 45]:

( ) tg ts , *<sup>G</sup>*

p

where *Lts* is the measured at-satellite radiance, *τ*(*μ*) ↑ is the direct (upward) target sensor atmospheric transmittance, *ρ*tg is the ground target reflectance, and *EG* is the global irradiance

Eq. (3) can be used to calculate the at-satellite radiance (*L*ts) for Landsat satellite images that contain gain and calibration offset values in the satellite image's header file, as described in the

where *E*0λ is the solar irradiance at the top of the atmosphere in Wm−2, found in the calibration file of the satellite image, *μ* is the cosine of the solar zenith angle (cos*θ*0), found on the satellite image header file, *τr* is the Rayleigh optical thickness, and *τa* is the AOT (which is unknown).

> tt

Although several equations to calculate *τr* have been developed, most radiative transfer equations used the equation provided by Moller [47], as indicated in Eq. (6), which is consid-

 m

 r

tm

*E*

atmospheric path radiance due to Rayleigh scattering in W/m2

*p*

*L L*

path radiance due to Mie scattering in W/m2

*2.1.1. Calculating path radiance*

144 Aerosols - Science and Case Studies

reaching ground.

Landsat-7 Handbook [46]:

The "*EG*" is calculated according to Eq. (4) [45]:

The "*τ*(*μ*) ↑" [42] is calculated in Eq. (5):

ered to be the most appropriate [34]:

*t*(m To calculate path radiance due to Rayleigh scattering(*L*pr), Eq. (7) is used [44]:

$$L\_{pr} = \left\{\frac{E\_{0\lambda}\cos\theta\_0 \,^\circ P\_r}{4\pi \left(\cos\theta\_0 + \cos\theta\_\nu\right)}\right\} \left\{1 - \exp\left[-\tau\_r \left(1/\_{\cos\theta\_0} + 1/\_{\cos\theta\_\nu}\right)\right]\right\} \left\{\imath o\_{3\uparrow} \,\imath o\_{3\downarrow}\tag{7}$$

where *Pr* is Rayleigh scattering phase function and *θv* is the sensor viewing angle, as found on the satellite image header file. *to*<sup>3</sup>↑ refers to the transmittance factor due to ozone in the upward direction, and *to*<sup>3</sup><sup>↓</sup> refers to the transmittance factor due to ozone in the downward direction, both of which are considered negligible and therefore have a value of 1 [44].

The Rayleigh scattering phase function (*Pr*) is calculated using Eq. (8) [44]:

$$P\_r = \frac{3}{4} \left[ + \cos^2 \left( 180 - \theta\_0 \right) \right] \tag{8}$$

#### *2.1.3. Calculating Mie scattering*

To calculate Mie scattering, Eq. (9) is used [47]:

$$L\_{pa} = \left\{\frac{E\_{0A}\cos\theta\_0\omega\_0 P'\_a}{4\pi\cos\theta\_0 + \cos\theta\_v}\right\} \left\{1 - \exp\left[-\tau\_a \left(\left(\frac{1}{\cos\theta\_0}\right) + \left(\frac{1}{\cos\theta\_v}\right)\right)\right] \right\} \tag{9}$$

$$\exp\left[-\tau\_r \left(\left(\frac{1}{\cos\theta\_0}\right) + \left(\frac{1}{\cos\theta\_v}\right)\right)\right]$$

where *wa* is the aerosol single scattering albedo as calculated by Waggoner et al. [48] and *P*′*a* is the phase function, as calculated by Gilabert et al. [44], Forster [29], and Turner et al. [45], provided in graph form or downloaded from AERONET stations taken with the Cimel sunphotometer (http://aeronet.gsfc.nasa.gov/).

#### *2.1.4. Steps for performing the image-based integrated method*

In order to retrieve AOT using the above equation and solve for *τa*, the below steps should be followed:

**Step 1:** Download the required Landsat image from the USGS website (glovis.usgs.gov) by selecting the area of interest and the date and time of the satellite overpass.

**Step 2:** Retrieve the solar zenith angle (*θ*0) at the overpass time, according to the header file on the satellite image.

**Step 3:** Retrieve the sensor viewing angle (*θv*) at the overpass time, according to the header file on the satellite image.

**Step 4:** Retrieve *E*0 according to satellite image calibration file.

**Step 5:** Calculate *μ*, which is the cos(*θ*0) using the variable derived from step 2.

**Step 6:** Retrieve *λC*, using satellite image calibration file.

**Step 7:** Calculate *τr* using Eq. (6), using the variable determined in step 6.

**Step 8:** Solve for *EG* using Eq. (4), using the variables derived from steps 4, 5, and 7. The *τa* (AOT) will be unknown.

**Step 9:**Find *t*(*μ*)↑ using Eq. (5), using the variables derived from steps 5 and 7. The *τa* (AOT) will be unknown.

**Step 10:** Calculate *Lts* from the satellite image header file using Eq. (3).

**Step 11:** Calculate *ρtg* from the surface reflectance values.

**Step 12:** Calculate *Lp* using Eq. (2), using variables from steps 8–11.

**Step 13:** Calculate *Pr* using Eq. (8), using variables from steps 2 and 5.

**Step 14:** Solve for *Lpr* using Eq. (7), using variables from steps 2, 3, 4, 7, and 13 and using 1 for *to*<sup>3</sup>↑ and *to*<sup>3</sup>↓.

**Step 15:** Enter *wa* which is between 0.73 and 0.87 in urban areas and 0.89–0.95 in remote areas.

**Step 16:** Enter *P*′*a* using phase function graph or from AERONET.

**Step 17:** Calculate *Lpa* using Eq. (9), using variables from steps 2, 3, 4, 7, 15, and 16.

**Step 18:** Solve for *τa* (AOT) using Eq. (1), using variables from steps 14 and 17.

#### **2.2. Image-based integrated method for GIS analysis**

In order to solve for AOT and conduct a GIS analysis, it may be necessary to simplify the equation due to the complexity of the RT equations and their logarithmic components. The equations for path radiance, Rayleigh scattering, and Mie scattering can be simplified and then combined into one equation, to be solved in an image-based software. After the AOT values are calculated for every pixel in the satellite image, the values can be exported into a GIS geospatial database, with the AOT value and coordinates. In order to display AOT on a GIS thematic map, the Kriging method can be used to interpolate the AOT values and then overlay the data into GIS vector maps or satellite images of the area of interest to display a spatial distribution of AOT. The maps can also be overlaid with vector data including roads, building lots, parks, green areas, stadiums, industrial areas, and municipal boundaries to identify AOT distribution over the area of interest. *In-situ* AOT values, if available within the area of interest, can be correlated with the AOT values from the GIS maps to perform an accuracy assessment for a specific location.

**Step 2:** Retrieve the solar zenith angle (*θ*0) at the overpass time, according to the header file on

**Step 3:** Retrieve the sensor viewing angle (*θv*) at the overpass time, according to the header

**Step 8:** Solve for *EG* using Eq. (4), using the variables derived from steps 4, 5, and 7. The *τa*

**Step 9:**Find *t*(*μ*)↑ using Eq. (5), using the variables derived from steps 5 and 7. The *τa* (AOT)

**Step 14:** Solve for *Lpr* using Eq. (7), using variables from steps 2, 3, 4, 7, and 13 and using 1 for

**Step 15:** Enter *wa* which is between 0.73 and 0.87 in urban areas and 0.89–0.95 in remote areas.

In order to solve for AOT and conduct a GIS analysis, it may be necessary to simplify the equation due to the complexity of the RT equations and their logarithmic components. The equations for path radiance, Rayleigh scattering, and Mie scattering can be simplified and then combined into one equation, to be solved in an image-based software. After the AOT values are calculated for every pixel in the satellite image, the values can be exported into a GIS geospatial database, with the AOT value and coordinates. In order to display AOT on a GIS thematic map, the Kriging method can be used to interpolate the AOT values and then overlay the data into GIS vector maps or satellite images of the area of interest to display a spatial distribution of AOT. The maps can also be overlaid with vector data including roads, building lots, parks, green areas, stadiums, industrial areas, and municipal boundaries to identify AOT distribution over the area of interest. *In-situ* AOT values, if available within the area of interest,

**Step 17:** Calculate *Lpa* using Eq. (9), using variables from steps 2, 3, 4, 7, 15, and 16.

**Step 18:** Solve for *τa* (AOT) using Eq. (1), using variables from steps 14 and 17.

**Step 4:** Retrieve *E*0 according to satellite image calibration file.

**Step 6:** Retrieve *λC*, using satellite image calibration file.

**Step 11:** Calculate *ρtg* from the surface reflectance values.

**Step 5:** Calculate *μ*, which is the cos(*θ*0) using the variable derived from step 2.

**Step 7:** Calculate *τr* using Eq. (6), using the variable determined in step 6.

**Step 10:** Calculate *Lts* from the satellite image header file using Eq. (3).

**Step 12:** Calculate *Lp* using Eq. (2), using variables from steps 8–11.

**Step 16:** Enter *P*′*a* using phase function graph or from AERONET.

**2.2. Image-based integrated method for GIS analysis**

**Step 13:** Calculate *Pr* using Eq. (8), using variables from steps 2 and 5.

the satellite image.

146 Aerosols - Science and Case Studies

file on the satellite image.

(AOT) will be unknown.

will be unknown.

*to*<sup>3</sup>↑ and *to*<sup>3</sup>↓.
