The Image-Based Integrated Method for Determining and Mapping Aerosol Optical Thickness

This chapter focuses on the development of an image-based integrated method for determining and mapping aerosol optical thickness (AOT). Using the radiative transfer (RT) equation, a methodology is developed using nonvariant targets with the revised darkest pixel method to improve the accuracy to calculate AOT over urban areas. This method can be applied to create thematic maps using GIS that present AOT values in a given area. The methodology is applied to Landsat TM/ETM+ satellite images of Limassol, Cyprus over a period of time to visualize and assess the AOT levels over the urban area. An accuracy assessment of the method shows a strong correlation between the in-situ AOT values from the sunphotometers and the AOT values derived from the image-based integrated method.


Introduction
Aerosol optical thickness (AOT) is defined as the degree to which aerosols prevent the transmission of light through absorption or scattering of light. AOT can be retrieved by using radiometers on board satellites [1][2][3], ground-based sunphotometers, and satellite sensors, such as the moderate resolution imaging spectroradiometer(MODIS), which provide AOT measurements on a daily and monthly basis [4,5]. However, the spatial resolution of the MODIS AOT data products is 10 × 10 km, which does not permit the identification of specific AOT distribution over urban areas or complex terrains [5][6][7]. Due to the variability of aerosols, atmospheric aerosol monitoring is difficult. Significant efforts to improve aerosol characterizations have included using in-situ measurements, ground-based remote sensing, and satellite observations [1][2][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 AErosol RObotic NETwork (AERONET, a program established by NASA) or field-based sun photometers [9,10].
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][12][13][14][15][16][17][18][19][20][21][22][23][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][19][20][21][22][23][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][27][28][29]. Several researchers [5,18,26,28,[29][30][31][32][33][34][35][36][37] found that using the radiative transfer and atmospheric modeling in conjunction with field measurements of aerosol optical thickness can yield more accurate atmospheric corrections instead of using simple image-based techniques.
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 thematic maps to display the AOT levels.

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.

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 (L p ) 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][41][42][43]: where L p is the path radiance (integrated band radiance measured in W/m 2 sr), L pr is the atmospheric path radiance due to Rayleigh scattering in W/m 2 sr [44], and L pa is the atmospheric path radiance due to Mie scattering in W/m 2 sr [44].

Calculating path radiance
The path radiance "L p " measured in Wm −2 μm −1 sr −1 , is calculated according to Eq. (2) [42,45]: where L ts is the measured at-satellite radiance, τ(μ) ↑ is the direct (upward) target sensor atmospheric transmittance, ρ tg is the ground target reflectance, and E G 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]: The "E G " is calculated according to Eq. (4) [45]: 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): 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]: where λ C is the central wavelength for each Landsat band.

Calculating Rayleigh scattering (L pr )
To calculate path radiance due to Rayleigh scattering(L pr ), Eq. (7) is used [44]: where P r is Rayleigh scattering phase function and θ v is the sensor viewing angle, as found on the satellite image header file. to 3↑ refers to the transmittance factor due to ozone in the upward direction, and to 3↓ 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].

Calculating Mie scattering
To calculate Mie scattering, Eq. (9) is used [47]: where w a 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/).

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 8: Solve for E G using Eq. (4), using the variables derived from steps 4, 5, and 7. The τ a (AOT) will be unknown.
Step 10: Calculate L ts from the satellite image header file using Eq. (3).
Step 11: Calculate ρ tg from the surface reflectance values.
Step 15: Enter w a 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.

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.

AOT retrieval using the image-based integrated method: the example of Limassol, Cyprus
The Geometric and radiometric correction is performed on the 11 Landsat-5 and Landsat-7 satellite images, as well as atmospheric correction using the modified DP method. The image-based integrated method is applied to the Landsat-5 and Landsat-7 images taken over specific days of 2010-2011. Eqs. (1-9) are used to retrieve AOT from bands 1-4 of the Landsat-5 and Landsat-7 satellite images, following the steps outlined. As the in-situ AOT measurements are at the 500 nm wavelength, Landsat band 1 is used to derive AOT levels in this method.
In this example, the model is run in the ERDAS Imagine and MATLAB software to generate an image that is consists of AOT values. The model is then applied to every surface using the Landsat band 1 satellite image, since the AOT required are in the 500 nm wavelength. The resulting image features the AOT values for every pixel in the image, as featured in Figure 2.
Areas with no data, due to cloud presence or negative AOT values, appear in white. All AOT values can be exported to ASCII files to create a georeferenced AOT dataset, which is then imported into the ArcGIS software.
Using the ArcGIS software, each pixel is converted to a point and each point is associated with the calculated AOT value. The white area of the image contains values of "no data" (Figure 3, left). To create a GIS thematic map showing the AOT distribution in the urban area of Limassol, the interpolation method is used in order to estimate the values from the "no data" regions. The ordinary Kriging interpolation tool is applied, where the unknown values of the "no data" area are interpolated using the weighted average of neighboring samples.   Once the interpolation is completed, an AOT thematic map is created, which is classified according to AOT values. The AOT values are displayed in different colors according to the AOT concentration and a legend of the AOT values is created. This facilitates the ability to visualize the concentration and distribution of AOT values over an urban area. Figure 4 features a map of Limassol following Kriging interpolation, which provides synoptic views of the AOT distribution. As is evident, high AOT levels are present throughout the city, especially in the industrial estates and on busy streets, with lower AOT levels present in parks, stadiums, and outside the city. In the example, the thematic map is overlaid with GIS vector data from the Lands and Surveys Department for Limassol to facilitate the identification of sources of AOT values within the area of interest ( Figure 5). Each polygon in the thematic map is connected to a GIS database, where data such as plot number and area information can be identified. In this way, the GIS database, combined with the AOT data, can create thematic maps to illustrate trends in different areas. In the example, the image-based integrated method is used to produce GIS maps indicating the AOT distribution over Limassol. In order to visualize high and low AOT levels, thematic maps are generated using colors ranging from blue to green to yellow to red for the specified AOT range. A separate thematic map is created for each of the 11 Landsat-5 and Landsat-7 satellite images. Figure 6 features the GIS thematic map using Landsat-7 satellite imagery from 24/6/2010.  Table 1 compares the AOT measurements taken on site using a sunphotometer against the AOT values derived from the GIS maps using the proposed algorithm by location and date, as well as the correlation coefficient (r) between the on-site AOT measurements and algorithmderived AOT values from GIS, for each location. In order to determine the accuracy of the GIS model, the root mean square deviation (RMSD) is also calculated between the on-site AOT measurements and algorithm-derived AOT values from GIS for each location. The results of the RMSD, correlation coefficients, and the coefficient of determination verify that the AOT derived from the GIS map correlate strongly with the on-site AOT values as measured by sunphotometers.

Conclusions
In this chapter, the image-based integrated method is developed to retrieve AOT values from satellite images and display these values in an urban area, using the example of Limassol, Cyprus. The findings show that the proposed integrated image-based method is able to accurately obtain AOT measurements from satellite images. The accuracy assessment shows strong agreement between the in-situ AOT values and the AOT values retrieved through the method described in the chapter. The proposed methodology is able to retrieve AOT using satellite images and visually display AOT over urban areas by using GIS to produce thematic maps. Also, the significance of the image-based integrated method is that the environmental data is not needed, as compared with other methods. This methodology is therefore an alternative to more sophisticated and complex methods of deriving AOT values from satellite images, especially where atmospheric or meteorological data are not accessible. Also, the methodology can be used with archived satellite images where environmental data are not available, thereby providing detailed information regarding spatial aerosol concentration overtime. Also, thematic maps can be generated to illustrate the distribution of AOT overtime in a specific area of interest and enabling trends to be identified.