Atmospheric Correction of Multispectral VNIR Remote Sensing Data: Algorithm and Inter‐sensor Comparison of Aerosol and Surface Reflectance Products

Optical imaging satellites, such as SPOT and Cartosat‐2S, provide visible/near infrared (VNIR) multispectral data at very high spatial resolution. The applications of these data sets are associated with precise mapping, monitoring, and change detection of Earth's surface, given that the measurements can be compensated for atmospheric effects. Existing atmospheric correction (AC) algorithms use visible and shortwave infrared channels and therefore cannot be used for AC of data from VNIR sensors. This article describes an algorithm for aerosol optical depth (AOD) retrieval and AC of VNIR imaging data. The AOD algorithm relies on the fact that for vegetated surfaces there exists a visible/NIR surface reflectance relationship due to the absorption of solar radiation by photosynthetic pigments in visible bands, while high reflectance in NIR bands governed by structural discontinuities in the leaves of healthy vegetation. We then describe how retrieved AOD is used to derive surface reflectance. To test the algorithm, the aerosol and surface reflectance products generated from 106 Cartosat‐2S data sets are compared with MODIS‐terra products. The algorithm significantly removes the haze from the images making surface feature visible. The comparison of Cartosat‐2S and MODIS‐terra AOD involving >1,500 data points shows good correlation of 0.95 with a relative difference of ≤25%. Similarly, the comparison of surface reflectance involving >4,500 data points shows good correlation ranging from 0.75 to 0.86 with a relative difference ranging from 24% to 37%. The normalized difference vegetation index shows a correlation of 0.89, with a relative difference of ≤18%. Results show that the given algorithm may be useful for AC of data from VNIR sensors.


Introduction
The space-borne multispectral visible/shortwave infrared (VSWIR) imaging spectrometers such as Moderate Resolution Optical Imaging Spectrometer-terra/aqua (MODIS-terra/aqua) (Salomonson et al., 1989), Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi National Polar-Orbiting Partnership (Cao et al., 2014), and Advanced Baseline Imager (ABI) onboard GOES-16 (Schmit et al., 2005) map radiance in multiple narrow spectral bands of solar reflective spectrum (400-2,500 μm). Surface reflectance features in the solar reflective spectrum are a good source of information about the chemistry and composition of Earth's terrestrial domains (Asner et al., 2017;Jetz et al., 2016;Ustin et al., 2004). However, since space-borne spectrometer measures signal at the top of atmosphere (TOA) that contains coupled information from both atmosphere and surface, therefore, atmospheric correction (AC) of TOA measurement is necessary to derive the surface reflectance before using it for surface studies. In general, if atmospheric state parameters like aerosol and water vapor concentration are known, the AC of any optical data set can be achieved with the aid of radiative transfer (RT) calculations (Berk et al., 2014;Mayer & Kylling, 2005;Ricchiazzi et al., 1998;Tanre et al., 1986;Vermote et al., 1997). For land surface, the existing ready-to-use AC software and algorithms such as Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) and Atmospheric Removal Program (ATREM) (Kruse, 2004;Perkins et al., 2012;Richter & Schläpfer, 2002) include a method for retrieving the aerosol amount using a dark pixel visible and SWIR reflectance ratio method (Kaufman et al., 1997;Levy et al., 2007). Thus with these models, the aerosol retrieval and therefore AC are only possible when the image contains bands in visible as well as in SWIR wavelengths. While for water surface, the measurements in two near-infrared bands are used to retrieve aerosol information (Ahmad et al., 2010;Antoine & Morel, 1999;Gordon, 1997;Sahay et al., 2014;Sarangi et al., 2015). Thus, the AC of data sets from VSWIR sensors (such as MODIS, VIIRS, and ABI) is within the reach of the present technology; consequently, these sensors have a rich history of applications in various fields, such as mineral exploration, vegetation, atmospheric, and ocean color studies. However, due to coarser spatial resolution (0.25-1.0 km), the data from these sensors have limited utility in certain fields such as urban developments, precision agriculture, bathymetry, coastal zone mapping, and other strategic applications where high spatial resolution is required.
Other kinds of sensors like SPOT-7 (Satellite Pour l-Observation de la Terre) and Cartosat-2S that provide images covering only visible and near-infrared (VNIR) part of the solar reflective spectrum contain a significant amount of spatial and spectral information about the Earth's surface. The disadvantage of these sensors is the lesser number of spectral bands (typically 3-5) in the VNIR region and no SWIR bands, making AC a challenging task with existing technology, but these sensors have the advantage of lower costs and a very high spatial resolution of the order of sub-meters to few meters. The high aerosol loading and adverse imaging condition (large sun or sensor zenith angle) magnify the atmospheric perturbation in TOA radiance, thus making AC a necessary task. For example, Figure 1 shows the Cartosat-2S natural color composite (NCC) image of a part of New Delhi, India. Due to significant aerosol scattering in the atmosphere, the NCC image ( Figure 1) appears hazy making the Earth's surface features unrecognizable in some parts of the image. Figure 2 shows the TOA reflectance spectrums for some vegetated pixels taken from data shown in Figure 1 along with typical vegetation spectrum convolved for Cartosat-2S bands. It is clear that due to atmospheric effects the shapes of TOA reflectance spectrums are completely different from the typical vegetation spectrum. Thus, to get useful information (qualitative and quantitative) about Earth's surface, AC of these data sets is required. Note that the typical spectrum shown in Figure 2 is derived using in situ surface reflectance spectrums acquired elsewhere, and magnitude wise may not represent the shown region. To the best of our knowledge, AC researches mainly remained focused on VSWIR sensors while the problem of AC of data acquired by VNIR sensors (SPOT-7 and Cartosat-2S) over land remained unaddressed. High spatial resolution data from VNIR sensors, if corrected for atmospheric effects, will be very useful for various applications involving precise mapping, monitoring, and change detection of Earth surface.
In this paper, we address the problem of AC of data acquired from VNIR sensors such as Cartosat-2S and SPOT-7. Using Cartosat-2S data sets, we first described an algorithm for aerosol retrieval and then described  Table 1. AC to derive "surface reflectance" assuming horizontal surface having Lambertian reflectance. After brief introduction in Section 1, Section 2 describes the features of VNIR instrument onboard Cartosat-2S. Section 3 describes the data sets and the study area used in this study. Section 4 describes the aerosol and AC algorithm for aerosol optical depth (AOD) and the surface reflectance retrieval, respectively. Section 5 presents the aerosol and atmospherically corrected Cartosat-2S products and shows their comparison with MODIS-terra products. Finally, Section 6 concludes the salient features of work.

Instrument Description
Cartosat-2S is an Indian satellite meant for applications like mapping, monitoring, and change detection of Earth's surface (Cartosat-2 brochure, n.d; Sharma et al., 2016). It contains one multispectral (MX) and one panchromatic (PAN) sensor having a high signal-to-noise ratio. The spacecraft is in a polar sun-synchronous orbit (orbital altitude of about 500 km) with approximately 9:30 hr of equatorial cross over time. The spatial resolution of MX and PAN data is around 2.0 and 0.65 m, respectively, with a swath of about 10 km. The salient features of Cartosat-2S MX payload are provided in Table 1. Figure 3 shows the spectral response functions (SRF) for MX bands and PAN band along with extraterrestrial irradiance at mean Earth-Sun distance.

Data Used and Study Area
The Cartosat-2S MX data sets acquired over different parts of the world covering parts of the Indian subcontinent (India, Nepal, Pakistan, and Bangladesh), China, and the United States are used in the present study. A total of 106 Cartosat-2S data sets acquired in January, July, November, and December months in the years 2016 and 2017 are used. Figure 4 shows the scene-center location of all 106 data sets. Each Cartosat-2S data set contains radiometrically calibrated radiance for all four bands. The Cartosat-2S data are routinely calibrated using ground measurements for reference targets by the data processing team (Sharma et al., 2016). The metadata file of each data set contains information about acquisition date/time, solar zenith angle, sensor zenith angle, azimuth angles, satellite altitude, mean elevation of the scene, and scene center latitude/longitude coordinates. As an example, Table 2 shows metadata information of two Cartosat-2S data sets acquired on 01 November 2016 and 03 October 2016 over the parts of New Delhi and Ahmadabad region, respectively. These two data sets are used to display the visual effect of AC on the TOA apparent reflectance image. The reason for choosing these two data sets among 106 data sets is due to significant aerosol loading in the atmosphere on these days of observation. To test the validity of our algorithm and the accuracy of Cartosat-2S products, the Cartosat-2S aerosol and surface reflectance products at all 106 globally distributed locations ( Figure 4) are compared with MODIS-terra products collocated in space and time domain. The MODIS-terra AOD products (DT aerosol retrievals, Collection 6.1), namely, "MOD04_3K MODIS/terra Aerosol 5-Min L2 Swath 3 km" at spatial resolution of 3 km, encompassing all 106 Cartosat-2S data set locations are used for comparison with Cartosat-2S AOD. For comparison of surface reflectance product, MODIS-terra surface reflectance product, namely, "MOD09GA MODIS/terra Surface Reflectance   (Kaufman et al., 1997;Levy & Hsu, 2015;Levy et al., 2007;Vermote & Wolfe, 2015). The equatorial passing time of Cartosat-2S and MODIS-terra are not the same but fall within 1-hr time window. Moreover, due to the agile platform of Cartosat-2S, the observation time difference may vary within a time window of few minutes to hours. In the present study, after careful analysis of acquisition times of all 106 data sets, it is found that the time difference between Cartosat-2S and MODIS-terra observations ranges from 30 min to 2 hr, out of which >90% observation falls within time difference of 1.5 hr and >80% of observation fall within 1-hr range. Due to the difference in time of observation, the precise spatiotemporal collocation of MODIS-terra and Cartosat-2S aerosol retrieval is not possible which may introduce error in comparison results particularly when the difference in observation time is very large. However, as mentioned earlier, most of the observation (>80% of total observations) fall within 1-hr time difference; therefore, no significant diurnal variation of AOD is expected in normal atmospheric conditions. The uncertainty in MODIS-terra AOD varies location wise because of location-dependent validity of visible/SWIR surface reflectance relationship used in dark-target aerosol algorithm (Kaufman et al., 1997;Levy et al., 2007). However, global validation studies using in situ AOD from Aerosol Robotic Network (Holben et al., 1998) have shown that about 67% of MODIS 3 km AOD retrievals fall within an error envelope (EE) of ±0.05 ± 0.20 (Remer et al., 2013) while about 66% of MODIS 10 km AOD retrievals fall within an EE of ±0.05 ± 0.15 (Levy et al., 2010). On average, the overall uncertainty in MODIS AOD is up to 20-30% (Levy et al., 2007) for the AOD value of 0.3 to 0.5. Thus in the absence of in situ AOD observations collocated with Cartosat-2S AOD retrievals, the comparison with MODIS-terra AOD product will give a good assessment of the performance of the presented aerosol algorithm.

Radiative Transfer
Based on RT theory of coupled planetary atmosphere-surface system (Chandrasekhar, 1960;Lenoble, 1985), an approximate but reasonably accurate vector RT code called Simulation of the Satellite Signal in the Solar Spectrum (6SV) has been developed in recent past (Kotchenova et al., 2006;Vermote et al., 1997). The relation expressing the TOA apparent reflectance for Lambertian surface is (Vermote et al., 1997) Here s , v , Δ , and refer to solar zenith, sensor zenith, relative azimuth angles, and wavelength, respectively. The various quantities in equation (1) are expressed in terms of equivalent reflectance defined by relation = L∕ s F 0 , where L refers to measured radiance, F 0 is the extraterrestrial irradiance at TOA, and s is the cosine of solar zenith angle. In equation (1), ra ( s , v , Δ , ) refer to the reflectance due to molecule and aerosols, and T( s , ) and T( v , ) refer to the downward and upward transmittance of the atmosphere, respectively. Symbol S( ) and T g refer to the spherical albedo of the atmosphere and gaseous transmittance (downward and upward direction both), respectively. Rearrangement of equation (1) gives where Given the satellite measurement, equation (2) can be solved for surface reflectance s ( ) if state of the atmosphere (gaseous transmittance and aerosol parameters) is known. Thus, AC of data from VNIR sensors such as Cartosat-2S requires accurate gas information and most essentially an algorithm to derive aerosol information from the data itself.

Gaseous Absorption
Out of many atmospheric gases, water vapor H 2 O, ozone O 3 , and oxygen O 2 are the main gases that show observable features in VNIR wavelengths under typical atmospheric conditions. Figure 5 shows atmospheric transmittance simulated for H 2 O = 10 mm and O 3 and O 2 gases content taken from standard tropical atmosphere. It is seen that Cartosat-2S MX bands are mainly affected by absorption due to O 3 and H 2 O while the PAN band is additionally affected by O 2 . In the present work, we used the standard tropical atmospheric model to calculate O 2 transmittance. Since the atmospheric column content of O 3 changes slowly with latitude and season, therefore, we used climatology data to calculate O 3 transmittance.

Aerosol Optical Depth Retrieval Algorithm 4.3.1. Physical Basis
It is a well-known fact that in the visible spectrum the vegetated surface shows low reflectance particularly at blue and red wavelengths due to absorption by photosynthetic pigments (mainly chlorophyll and carotenoids); therefore, reflectance in red and blue channels is generally less than 15% for vegetated surfaces (Peñuelas & Filella, 1998). However, from Figure 5, it is clear that for moderate aerosol loading the apparent reflectance due to aerosols at TOA may be up to 20% in blue wavelengths. For a very thick aerosol layer, it can be more than 20% depending on the absorption and scattering nature of aerosols. The major challenge in any aerosol retrieval algorithm is to decouple the TOA reflectance into the contribution from the surface and atmospheric constituents and then separating the molecular scattering contribution from 10.1029/2019EA000710 atmospheric contribution giving the aerosol contribution. For water surfaces due to the absorption of solar radiation, the surface reflectance in red and NIR wavelengths are nearly zero, therefore, assuming the entire signal in these wavelengths as the atmospheric contribution is a good approximation for deriving aerosol information. However, land surfaces show significant non-zero reflectance depending on the surface types. MODIS-aqua/terra dark-target aerosol algorithm estimates the visible bands (blue and red wavelength) surface reflectance using the measurement in SWIR band (Kaufman et al., 1997;Levy & Hsu, 2015). However, due to the unavailability of SWIR measurements in VNIR sensors like Cartosat-2S, the dark-target aerosol algorithm cannot be used. Observing that the effect of aerosol scattering is very small (excluding the case of dust events) on TOA measurements in NIR wavelengths (see Figure 5), we propose the use of NIR channel to characterize surface reflectance in visible wavelengths for reference target types while the process of aerosol retrieval. For vegetated targets, in visible channel (blue and red wavelengths) reflectance is low due to absorption from photosynthetic pigments, while in NIR, where there are no significant absorption features by pigments, the magnitude of reflectance is governed by structural discontinuities encountered in the leaves of healthy vegetation (Peñuelas & Filella, 1998). This forms the physical basis of visible/NIR surface reflectance relationship over the green vegetated target. Moreover, Lorenzen and Jensen (1988), by conducting the field measurements of canopy reflectance of the vegetation in the blue, green, red, and NIR wavebands, have shown that the surface reflectance in these wavebands is a strong function of green biomass. In fact, the surface reflectance in blue and red wavelengths decreases with increasing green biomass, while for NIR wavelength the nature of relation gets inverted (Lorenzen & Jensen, 1988). Based on this discussion, the physical bases of our aerosol retrieval algorithm are as follows: (a) Solar radiation in NIR wavelength is very less sensitive to aerosols, and (b) for vegetated target there exists a visible/NIR surface reflectance relationship that depends on greenness of the target defined by quasi-normalized difference vegetation index (NDVI ′ ): where RCR refers to Rayleigh corrected reflectance. This method is based on the vegetated land surface modeling used in the NASA Deep Blue aerosol retrieval algorithm (Hsu et al., , 2017(Hsu et al., , 2019.

Visible/NIR Surface Reflectance Relationship
To determine the spectral surface reflectance relationships between visible and NIR channels of Cartosat-2S, the best strategy is to use the atmospherically corrected Cartosat-2S surface reflectance data generated using in situ aerosol measurements. However, in situ aerosol observations collocated with Cartosat-2S data were not available; therefore, hyperspectral data from Airborne Visible/Infrared Imaging Spectrometer Next Generation (AVIRIS-NG) acquired during late 2015 and early 2016 under Indian campaign are used (Mishra et al., 2019;Thompson et al., 2018Thompson et al., , 2019. The in situ AODs were measured using sun-photometers during the AVIRIS-NG flights. The AVIRIS-NG data sets were corrected for atmospheric effects using measured AOD, and then some of the pure vegetated and soil spectra were chosen. We selected six different types of dense vegetation spectrum and six different soil type spectra. By mixing these soil and vegetation spectra, we synthesized a total of 630 spectra (Figure 6a) with different values of NDVI (normalized difference vegetation index). Then by convolution of these spectra with Cartosat-2S MX SRF (Figure 6b), we estimated the TOA reflectance for NIR and red/blue channels using RT calculations (i.e., using equation (1)) assuming two values of AOD at 0.55 μm equal to 0.1 and 0.5 and continental aerosol type. Based on these results, we developed an empirical relationship that estimates the surface reflectance at visible channels (0.485 and 0.649 μm) as a function of the reflectance in the NIR band (i.e., 0.808 μm) and NDVI ′ values. The empirically developed relationships are given by Here RCR refers to Rayleigh corrected reflectance. The model coefficients a, b, c, d, e, and f are calculated according to two NDVI ′ classes. In addition, aerosol retrieval is only done if NDVI ′ ⟩0.1 and ESR 0.485 ⟨0.15. Table 3 shows the resulting coefficients for each class.
To assess the performance of relationships in equation (5), we used four AVIRIS data sets acquired over Kalaburgi and Muddur in Karnataka province, Patna in Bihar province, and Zawar in Rajasthan province of India. The AVIRIS-NG data at 5.0-nm spectral resolution are convolved to Cartosat-2S MX SRF. Figure 7 shows the NCC image of all four data sets. Kalaburgi and Muddur are mainly urban and mixed-agriculture sites, respectively. Patna is a city situated at the bank of Ganga River, and the region shown in the image is an urbanized area with some portion of the Ganga River catchment area. Zawar is a geological site mainly of mineralogical interest. Figure 7 shows the scatter plots between atmospherically corrected surface reflectance (SR 0.485 and SR 0.649 ) and the estimated surface reflectance (ESR 0.485 and ESR 0.649 ) for blue and red channel using NIR measurements (RCR 0.808 ) in visible/NIR relationship (equation (5)). The linear regression statistics of comparison are provided in Table 4. High correlation coefficients, r, are found at all four stations with r ranging between 0.88 to 0.96 and 0.98 to 0.99 for the blue and red bands, respectively. The slope of the linear regression lines shows variation about unity by 21% and 24%, respectively, for blue and red bands. The relative error of estimates for all four sites ranges from 10% to 20% and 7% to 16% for the blue and red bands, respectively. The maximum root mean square error observed is 0.012 and 0.019 for the blue and red bands, receptively. An error in surface reflectance of the order of 0.01 leads to an error of up to 20% in inverted AOD (Mishra, 2018).

Aerosol Look-Up Tables
The inversion of atmospheric contribution to AOD requires RT calculations of various radiative quantities, namely, path reflectance, aerosol transmittance, and spherical albedo, for different aerosol types and aerosol loading. Since on-line RT simulations are computationally very time consuming, therefore to avoid on-line RT simulations, pre-computed look-up tables (LUTs) of radiative quantities are used in operational aerosol retrieval code. For RT calculations, the optical properties of seven aerosol models are used. Out of these, five aerosol models are generated by external mixing of three basic aerosol types, namely, water-soluble aerosol, dust-like aerosol, and shoot aerosol, assuming number density as given in Table 5. The parameters required in Mie calculation to derive optical properties of the basic aerosol types are provided in Table 6 with an assumption of log-normal particle-size distribution function: Here r is the particle radius, r m is the mean radius, and is the width of distribution function. It is to be noted that the physical properties of the basic aerosol components described in Table 6 are from WMO-WCP112   (1986); however, updated physical properties of these components as well as of many more aerosol models can be referred from Hess et al. (1998), to generate improved LUTs which may closely represent real aerosols. The parameters for the rest two models, namely, biomass burning and desert dust, are taken from Dubovik et al. (2002) (Dubovik et al., 2002). In the present version of the algorithm, 6SV RT code is used for generating LUTs. For every aerosol model, LUTs of path reflectance, transmittance, and spherical albedo corresponding to AOD = 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0 at 0.55 μm are generated. LUTs for path reflectance are generated for discrete values of solar zenith, sensor zenith, and relative azimuth, while LUTs for transmittance are generated for discrete values of solar and sensor zenith. The spherical albedo for every value of AOD is also tabulated. Four sets of these LUTs corresponding to each Cartosat-2S waveband are generated and saved for operational use.

Retrieval Method
For aerosol retrieval, we first divide the whole image into 25 × 25 pixel boxes (for Cartosat-2S MX data it is approximately 50 m by 50 m box). For each box, we first classify pixels based on NDVI ′ , and then we retrieve aerosol properties by taking advantage of the spectral relationship (equation (5)) of surface reflectance between visible bands centered at 0.485 and 0.649 μm and NIR band centered at 0.808 μm to account for the effect of the ever-changing dynamics of vegetation phenology on the surface reflectance. Very bright and thick cloud contaminated pixels are avoided using a cloud mask generated using the threshold test in visible band reflectance. The remaining pixels are sorted based on blue channel reflectance, and 50% brightest  and 20% darkest values are also masked to avoid cloud contamination. It is to be noted that the pixels with NDVI ′ ⟨0.1 are also excluded to avoid pure soil or any inland water body pixel. The remaining pixels are used in the aerosol retrieval algorithm to derive AOD value for the 25 × 25 pixel box. After estimating the surface reflectance at 0.485 and 0.649 μm for the remaining pixels using relations in equation (5) for vegetated pixels in each box, we estimate mean surface reflectance for each box. Using these surface reflectance values and pre-computed LUTs, modeled TOA reflectance at 0.485 and 0.649 μm are computed for given sun-sensor geometry and for different values of indexed AOD, and single scattering albedo. These 0.485 and 0.649 μm modeled TOA reflectance are then compared to Cartosat-2S MX measured TOA reflectance to match the appropriate values of AOD at 0.55 μm and aerosol model (i.e., single scattering albedo) for each box. Since AOD at 0.55 μm is retrieved at a spatial resolution of 1/25 of the original resolution, the cubic interpolation is performed to re-sample derived aerosol optical thickness map at original spatial resolution.

Atmospheric Correction
The method for deriving surface reflectance from Cartosat-2S MX TOA radiance data is shown in Figure 8, which consists of following major steps: (a) One combined transmittance spectrum for O 3 and O 2 is calculated using sun-sensor geometry and climatological values for O 3 and O 3 content. One H 2 O transmittance spectrum is generated using total column H 2 O content from forecast reanalysis data. Both transmittance spectra are convolved for four MX bands. (b) The TOA 11-bit digital count image for each MX band is converted to TOA apparent reflectance using relation where DN, s , and F 0 refer to digital count, cosine of solar zenith, and extraterrestrial irradiance corrected for earth-sun distance, respectively. (c) The TOA reflectance is corrected for gaseous absorption by dividing TOA of all bands with gaseous transmittance computed in step (a). (d) The entire image is divided into 25 × 25 pixel boxes. AOD for each 25 × 25 box of the scene is retrieved by the method described in Section 4.3, and re-sampling of AOD map to original resolution using cubic interpolation is done. (e) Using estimated AOD and aerosol model, the radiative quantities, namely, the atmospheric reflectance, spherical albedo, and the total upward and downward transmissions, are estimated with the aid of pre-computed LUTs for all MX bands. (f) Finally, the surface reflectance for each MX band is derived based on equation (2). In our computer program for AC, the first three steps need to be performed only once at the beginning of the program execution, and the last three steps are performed for each pixel. The output consists of the AOD map at 0.55 μm and surface reflectance for all MX bands. It is important to note that the two-step approach for AC (first estimate AOD and then use that to estimate surface properties) described above introduces inconsistency in that the eventual AC surface reflectance may not be consistent with the surface reflectance assumed in the AOD retrieval part even if they are done at the same resolution.

Visual Improvement in Images After Atmospheric Correction
Using aerosol retrieval and AC algorithm, the AOD and surface reflectance products are generated for all 106 Cartosat-2S acquired over locations shown in Figure 4. Figures 9-14 show the results (aerosol and surface reflectance) obtained from the Cartosat-2S data sets acquired over Delhi and Ahmedabad (see Table 1). Figure 9b shows the retrieved AOD map from Cartosat-2S MX data acquired on 01 November 2016 over a part of Delhi. It is clear from the AOD map that the left part of the image shows a high AOD (>1.2) relative to the right part of the image. On average, whole image shows heavy aerosol loading with mean AOD > 1.2. Heavy  aerosol loading on 01 November 2016 may have occurred due to the transportation of dust aerosols from the regions of Afghanistan-Pakistan border over which a heavy dust storm event was detected in MODIS NCC image (Figure 9a) on 01 November 2016. In addition to this, the practice of agriculture residue burning in November month in adjoining northern states (Haryana and Punjab) of Delhi also has introduced a large amount of carbonaceous aerosol as shown by MODIS-aqua hot spots (red dots) in Figure 9a. Similar qualitative information is seen from Cartosat-2S MX NCC image shown in Figure 10a, where appreciable haziness, especially in the left part of the image, is apparent. Figures 10a and 10b show TOA apparent reflectance and retrieved surface reflectance NCC images, respectively, derived from Cartosat-2S MX data using the present algorithm. Figures 10c and 10d show TOA and surface reflectance, respectively, for the Subset 1 (red square in Figure 10a), while Figures 10e and 10f show TOA and surface reflectance, respectively, for the Subset 2 (red square in Figure 10a). From Figure 10, it is clear that derived surface reflectance (right panels) is free from haze as compared to TOA NCC images (left panels) and therefore the surface features especially in the left region that was obscured in TOA reflectance image became visible after AC. Figure 11 shows a quantitative comparison of TOA apparent reflectance and derived surface reflectance from Cartosat-2S MX data acquired on 01 November 2016 (shown in Figure 10) over some of the vegetated and soil type pixels. The standard vegetated and soil spectra are also shown in Figure 11 to compare the shapes of the derived spectra. From Figures 11a and 11b, it is clear that over soil and vegetated targets, respectively, the shapes of derived surface reflectance spectra are in good agreement with standard spectra. While TOA apparent reflectance spectra show very high reflectance values and shapes are completely different from the standard spectra. Figure 12 shows the histogram analysis of TOA reflectances and the derived surface reflectances for Bands 1, 2, and 3 of Cartosat-2S MX data acquired on 01 November 2016. It is clear from the figure that after AC the mean values (standard deviation values) of derived surface reflectance are significantly lower (larger) than that of TOA reflectance for all three bands. It is clear that AC has appreciably increased (by a factor of >2.4) the standard deviation (increased contrast) compared to TOA reflectance which shows low standard deviation due to homogeneity introduced by heavy atmospheric aerosol layer. The maximum change in mean reflectance value is observed for Band 1 (blue waveband), which is expected as blue wavelengths are most sensitive to aerosol scattering.  Similarly, we processed the Data Set 2 acquired by Cartosat-2S on 03 October 2017 that covers part of Ahmedabad. Figure 13 shows a visual comparison of TOA reflectance and derived surface reflectance NCC images. From Figure 13a we can see that significant region of the image is contaminated by bright haze and thick clouds and therefore the surface features underneath the thick haze layer are almost invisible (see region encompassed by redlined boxes and adjoining region in Figure 13a). On the other hand, in derived surface reflectance NCC image ( Figure 13b) the haze is significantly cleared and the surface features became visible. Figures 13c and 13d, and 13e and 13f, show the effect of AC for very hazy regions represented by Subsets 1 and 2, respectively. It is clear from these figures that except the regions with a very thick cloud where no information about surface reaches TOA, the surface features become visible and clear from the haze. It is well known that the effect of high aerosol scattering on TOA data is to decrease the heterogeneity of the surface features and generally increase the brightness of the scene. Figure 14 shows the histogram analysis of TOA reflectances and the derived surface reflectances for Bands 1, 2, and 3 of Cartosat-2S MX data acquired on 03 October 2017. Similar to the previous case, here also after AC, the mean values (standard deviation values) of derived surface reflectance are significantly lower (larger) than that of TOA reflectance for all three bands. The AC has appreciably increased (by a factor of >1.7) the standard deviation compared to TOA reflectance which shows low standard deviation due to homogeneity introduced by heavy atmospheric aerosol layer.
Here also the maximum change in mean reflectance value is observed for Band 1 (blue waveband), which is expected. Therefore, these statistical parameters (decreased mean and increased standard deviation after atmospheric correction) and visuals presented in Figures 9-14 prove that our algorithm effectively removes atmospheric effects (haze due to aerosols), thereby increasing the contrast of the surface features.

Inter-sensor Comparison of Aerosol and Surface Reflectance Product
The qualitative analysis in the previous subsection clearly shows the improvement in Carosat-2S MX image quality, the shape of the spectrum, and the statistical parameters (mean and standard deviation of reflectance), after AC. However, the validation of aerosol and surface reflectance products with other sources is needed to establish the validity of Cartosat-2S products and the presented algorithms. Due to lack of in situ measurements of AOD and surface reflectance collocated in time and space with Cartosat-2S observations, in the present work, MODIS-terra aerosol product at the spatial resolution of 3 km and surface reflectance product at a spatial resolution of 500 m are used to compare Cartosat-2S products. For more information on MODIS-terra products used in this work, see Section 3.

Comparison Methodology
The Cartosat-2S AOD and surface reflectance products are at different and very fine spatial resolution of 50 and 2 m, respectively, as compared to the MODIS-terra product. Moreover, the projection of Cartosat-2S products is in UTM, while MODIS-terra aerosol and surface reflectance products are provided with per pixel latitude and longitude information as separate layers for a given tile. Due to all these differences, proper collocation, down-sampling, and data selection criteria of Cartosat-2S data must be defined for faithful comparison with MODIS-terra product. The reason for selecting the MODIS-terra product is because Cartosat-2S equatorial pass time is very near to that of MODIS-terra. As the first step, MODIS-terra AOD is geo-referenced to geographic Lat/Lon projection using the latitude and longitude layers provided with each data set. It is to be noted that before geo-referencing the AOD layer, the pixels with quality flag equal to 2 or 3 (i.e., very good or best quality retrievals) are only retained; all other pixels are masked with a fill value. In the second step, the Cartosat-2S data are converted from UTM to geographic Lat/Lon projection. For every MODIS Lat/Lon grid box, the corresponding Cartosat-2S AOD pixels are identified based on latitude and longitude information. The pixels having valid AOD retrieval among all of the identified pixels are counted. If valid AOD pixel count is over ≥50% of total pixels, the average value of valid pixel values is filled. Doing so the 10.1029/2019EA000710 down-sampled Cartosat AOD at 3 km and surface reflectance at 500 m spatial resolution is obtained which are perfectly collocated with MODIS-terra geo-reference AOD and surface reflectance products, respectively. Finally, pixel by pixel values of AOD (at 3-km resolution) and surface reflectance (at 500-m resolution) from 106 down-sampled Cartosat-2S products and corresponding values from MODIS-terra geo-referenced product are tabulated, and linear regression statistics are estimated. The normalized difference vegetation index (NDVI) from collocated surface reflectance of Cartosat-2S and MODIS-terra is also calculated for comparison.

Comparison of Aerosol Optical Depth
The comparison of Cartosat-2S AOD ( c ) with MODIS-terra AOD ( m ) is plotted in Figure 15. The linear regressions statistics of these validations are provided in Table 7. High correlation coefficients, r, are found with a value of r ≥ 0.95. Slopes of the linear regression lines show small variation about 1 by 3-6%, probably caused by the use of slightly different aerosol optical properties in the Cartosat-2S aerosol algorithm as compared to MODIS-terra aerosol algorithm where a dynamic mixture of fine and coarse mode aerosols is used. For a similar reason or due to uncertainty in empirical surface reflectance (due to the BRDF effect) model, the Cartosat-2S AOD is slightly negatively biased; however, the observed bias is within the expected uncertainty of 25-30% for such algorithms (Levy et al., 2007). The slope of regression curve and negative bias of Cartosat-2S AOD relative to MODIS-terra AOD retrieval may be improved by using the location-based aerosol model selection scheme with dynamic mixing of course and fine mode aerosols with different scattering and absorbing properties. The relative difference of Cartosat-2S AOD estimates ( ∕⟨P⟩ m ), that is, the standard error of regression ( ) compared to the average observed AOD from MODIS-terra (⟨P⟩ m ), is found to be 25%. The result shown in Figure 15 and in Table 7 established that c have an average deviation of around 25% with respect to m at 3-km spatial resolution; however, from scatter plots in Figure 15, it is evident that the deviation of c relative to m is a function of m itself. Therefore, scatter plot in Figure 15 is somewhat heteroscedastic in nature which may lead to bias in linear regression fit parameters shown in Figure 15 and in Table 7. This makes necessary to derive uncertainty in c given by variable Δ c as a function of m . For this all 1,500 data points are first sorted in increasing order of m and then are divided into 15 bins with equal number of data point (approximately 100 data points in each bin) entries. For each bin, the modulus of difference of AODs | c − m | is calculated. It is found that | c − m | shows near-normal distribution; therefore, the first standard deviation of this value (i.e., the value of | c − m | below which 68% of retrieval falls) is defined as the uncertainty value Δ c for the corresponding bin . Doing so 15 values of Δ c are estimated corresponding to average value of m for each bin. The results are plotted in Figure 16, where from the linear regressions analysis shows that Δ c = 0.03 + 0.14 m , which means the uncertainty in Δ c increases with increasing value of m . The extensive global validation  of MODIS-terra aerosol retrievals with in situ AOD ( i ) shows an uncertainty of Δ m = ±(0.05 + 0.15 i ) in m (Levy et al., 2010;Remer et al., 2013). Using the equation for Δ m and Δ c , it is found that the error in retrieved AOD at low values ( i = 0.25) may be up to 65% of i . However, the error in retrieved AOD decreases to 30% for moderate AOD values ( i = 0.5) and then again starts increasing with an error of around 40% for high AOD values ( i = 1.0). These observations show that the Cartosat AOD retrievals may have high errors when actual AOD is too low (<0.25) and too high (>1.0). It is to be noted that these error estimates hold valid when Δ m and Δ c are of the same sign; however, if Δ m and Δ c are of opposite sign, the actual error will be less.

Comparison of Surface Reflectance
The Cartosat-2S surface reflectance (SR c,blue , SR c,green , SR c,red , SR c,nir ) comparison with MODIS-terra (SR m,blue , SR m,green , SR m,red , SR m,nir ) at spatial resolution of 500 m are plotted in Figure 16. The linear regressions statistics of these validations are provided in Table 7. High correlation coefficients, r, for all bands are found with values in range 0.75 to 0.86. The red band shows maximum correlation with r equal to 0.86, which is expected as red channels of both sensors have similar center wavelengths (for SRF and wavelengths of MODIS-terra bands, see Xiong et al., 2006). Variation of slopes of the linear regression lines about unity by 19-28% may be associated with the difference in full width half maximum (FWHM), center wavelengths, and the SRF of Cartosat-2S and MODIS-terra bands. The relative error of estimates ( ∕⟨P⟩ m ), that is, the standard error of regression ( ) compared to the average observed surface reflectance from MODIS-terra (⟨P⟩ m ), is found to be in range 24% to 37%. The high relative differences in MODIS and Cartosat-2S surface reflectance are also associated with the difference in FWHM and the shape of SRFs of the two sensors. Moreover, this discrepancy will be surface-type dependent. The comparison of NDVI obtained from Cartosat-2S and MODIS-terra surface reflectance is plotted in Figure 17. The linear regressions statistics of NDVI comparison are provided in Table 7. Very good correlation r with value 0.89 and slope near unity is found. The relative error of estimates for NDVI is found to be ≤18%, which is due to the large difference in the center wavelength of NIR channels.

Summary
Ability of aerosol retrieval and accurate AC of data acquired from low cost and high spatial resolution ( 1-2 m) optical VNIR sensors such as Cartosat, World-view-4, and SPOT-7 will keep these sensors relevant in remote sensing of Earth's surface and atmosphere. In fact, it will make high-resolution VNIR sensors complementary to advanced but coarser VSWIR multi/hyperspectral instruments. Here we described an aerosol retrieval and AC algorithm for multispectral VNIR data acquired over land in fewer spectral bands (three to five bands). The AOD, surface reflectance, and NDVI obtained by processing Cartosat-2S MX data at 106 globally distributed locations have been compared with corresponding MODIS-terra products collocated in time and space to test the performance of our algorithms. The comparison of Cartosat-2S and Figure 17. (a-d) The comparison of Cartosat-2S surface reflectance (SR c,blue , SR c,green , SR c,red , SR c,nir ) with MODIS-terra (SR m,blue , SR m,green , SR m,red , SR m,nir ) for blue, green, red, and NIR bands, respectively, at spatial resolution of 500 m. (e) Comparison of NDVI derived from Cartosat-2S and MODIS-terra surface reflectances at 500-m spatial resolution.