the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Improved mean field estimates from the Geostationary Environment Monitoring Spectrometer (GEMS) Level3 aerosol optical depth (L3 AOD) product: using spatiotemporal variability
Sooyon Kim
Yeseul Cho
Hanjeong Ki
Seyoung Park
Dagun Oh
Seungjun Lee
Yeonghye Cho
Jhoon Kim
Wonjin Lee
Jaewoo Park
Ick Hoon Jin
This study presents advancements in the processing of satellite remote sensing data, focusing mainly on aerosol optical depth (AOD) retrievals from the Geostationary Environment Monitoring Spectrometer (GEMS). The transformation of Level2 (L2) data, which includes atmosphericstate retrievals, into higherquality Level3 (L3) data is crucial in remote sensing. Our contributions lie in two novel improvements to the processing algorithm. First, we improve the inversedistanceweighting algorithm by incorporating quality flag information into the weight calculation. By assigning weights that are inversely proportional to the number of unreliable grids, the method can provide more accurate L3 products. We validate this approach through simulation studies and apply it to GEMS AOD data across various regions and wavelengths. The use of quality flags in the algorithm can provide a more accurate analysis of remote sensing. Second, we employ a spatiotemporal merging method to address both spatial and temporal variability in AOD data, a departure from previous approaches that solely focused on spatial variability. Our method considers temporal variations spanning previous time intervals. Furthermore, the computed mean fields show similar spatiotemporal patterns to previous studies, confirming their ability to capture realworld phenomena. Lastly, utilizing this procedure, we compute the mean field estimates for GEMS AOD data, which can provide a deeper understanding of the impact of aerosols on climate change and public health.
 Article
(23829 KB)  Fulltext XML
 BibTeX
 EndNote
In satellite remote sensing missions, observed data are processed at different levels. Using retrievals from the atmosphericstate (Level 2; L2), L2 aerosol optical depth (AOD) products are regridded to Level3 (L3) products through the processes of gap filling and filtering out noises (Cressie, 2018). We first introduce the theoretical background of the mean field algorithm used to generate L3 data applied to AOD retrievals from the Geostationary Environment Monitoring Spectrometer (GEMS) instrument. We consider an oversampling method for generating L3 AOD data, an inversedistanceweighting (IDW) algorithm, and a modified mean field algorithm while accounting for spatiotemporal variability in AOD data in the algorithm.
Aerosols play a critical role in radiative forcing, climate change, and air quality (Brauer et al., 2015; Charlson et al., 1992; Stocker, 2013; Kaufman et al., 2002). They directly change the planetary albedo by reflecting solar radiation and absorbing terrestrial radiation, affecting the radiation balance. Indirectly, as cloud condensation nuclei, aerosols modify cloud properties and increase cloud droplet concentration, impacting solar radiation and cloud albedo (Alexander et al., 2013). Aerosols affect human health and air quality, especially in regions affected by longrange transport or regions with heavy aerosol emissions due to rapid industrialization and a high population density. These effects are linked to cardiovascular, respiratory, and allergic diseases, as well as increased mortality rates (Pöschl, 2006; Tager, 2004).
Additionally, high aerosol concentrations can severely reduce visibility, leading to hazardous weather conditions, such as haze, smog, and dust storms (Charlson, 1969; Chen and Tsai, 2001). Thus, understanding the multifaceted impacts of aerosols is crucial for addressing issues concerning climate change, public health, and environmental visibility. The distribution of aerosols is characterized by its complexity, leading to increased uncertainty in determining the radiativeforcing effects of aerosols (Chen et al., 2022). Analyzing the spatiotemporal distribution of aerosols remains crucial for developing air pollution control policies and understanding the climate impacts of aerosols. Although accurate aerosol optical properties (AOPs) and their vertical profiles can be obtained from groundbased measurements with a high temporal resolution, the AOPs only represent localscale variability with limited spatial coverage. Unlike groundbased instruments, the monitoring of AOPs on regional and global scales has been conducted using satellite measurements.
A previous study by Park et al. (2023) focused on AOD retrievals by considering spatial variability. Specifically, Park et al. (2023) used the IDW algorithm to regrid L2 products and estimated the mean field of L3 products by considering spatial variability. Compared to this previous work, our contributions are as follows. First, we integrate quality flag information into the IDW algorithm to rule out unreliable grid points. In this study, we employ the IDW algorithm to compute mean field estimates for AOD data. The IDW method is a widely used gapfilling algorithms that imputes missing observations through a linear combination of neighboring observations. Integrating quality flags (QFs) as weights into the IDW algorithm allows for a quantitative assessment of the influence of data from various quality levels on the final product. It can mitigate the impact of lowquality data and lead to a lower mean squared error (MSE) when implementing the IDW algorithm, as demonstrated in Sect. 4. By considering variability in L2 AOD products, we can obtain more reliable L3 AOD products in this step. Second, we use the spatiotemporal merging method (Kikuchi et al., 2018) to obtain L3 AOD mean field estimates. Numerous studies have demonstrated that aerosol optical depth (AOD) exhibits significant spatiotemporal variability due to natural and anthropogenic factors. This is particularly evident in regions like northwestern China, where understanding spatiotemporal dynamics is essential for effective atmosphericpollution management and control as aerosol concentrations are heavily influenced by seasonal and graphical variations (Zhang, 2023). Moreover, the characteristics of the data collection device, the satellite, play a crucial role. The AOD data analyzed in this study were obtained from the Geostationary Environment Monitoring Spectrometer (GEMS). When collecting data using such devices, it is imperative to consider spatiotemporal variability to ensure the high quality of the aerosol data. For instance, a study conducted in East Asia optimized the spatiotemporal ranges used for validating satellite products, such as total ozone and NO_{2}, by leveraging longterm data from both groundbased and satellite observations (Park et al., 2020). Incorporating quality flag information can reduce the influence of lowquality data on the IDW algorithm, resulting in a lower MSE, as demonstrated in our simulation study. Nevertheless, it remains crucial to identify and exclude pixels with unreliable AOD values as they can introduce substantial uncertainties into mean field estimates. In our analysis, we apply a filter with a cloud radiance fraction (CRF) exceeding 0.4 to remove pixels, ensuring that the included AOD values are not significantly impacted by cloud contamination. We observe that our method can provide smoother AOD surfaces than the simple averaging method that does not consider spatiotemporal variability.
The outline of the remainder of this paper is as follows. In Sect. 2, we describe the GEMS data used in our analysis. In Sect. 3, we describe our method for computing the mean field of L3 AOD products. In Sect. 4, we conduct simulation studies to validate our method. We apply the proposed method to GEMS data in Sect. 5. We conclude with a discussion in Sect. 6.
The GEMS is the first ultraviolet–visible (UV–Vis) hyperspectral satellite instrument aboard Geostationary Korea MultiPurpose Satellite2B (GEOKOMPSAT2B), launched on 19 February 2020. Its mission is to monitor air quality across Asia (5° S–45° N, 75–145° E) using high temporal resolutions (1 h) and high spatial resolutions (3.5 × 7.7 km^{2} for Seoul, South Korea) using hyperspectral measurements in the 300–500 nm range.
The GEMS aerosol retrieval algorithm (AERAOD) retrieves AOD, single scattering albedo (SSA), and aerosol layer height (ALH) using GEMS L1 data from six wavelengths (354, 388, 412, 443, 477, and 490 nm). This algorithm solves the issues regarding the limited degree of freedom for signal problems in the GEMS wavelength range by utilizing twochannel inversion to retrieve initial guesses for AOD and SSA and then input them into the optimalestimation method. This retrieval method was tested for sensitivity to aerosol absorption in the UV–Vis region and using Ozone Monitoring Instrument (OMI) Level1 data (Kim et al., 2018; Go et al., 2020a, b). Initially developed from synthetic OMI data by Kim et al. (2018) and Go et al. (2020b), the operational version was later improved by Cho et al. (2023) using real GEMS Level1 data. An update to the aerosol algorithm, version 2.0, was released in November 2022 and included the reprocessing of earlier data.
In this study, the following variables are used for calculating L3 AOD mean fields (Table 1).
In addition to the solar zenith angle (SZA) and viewing zenith angle (VZA), due to the unavailability of cloud fraction in the GEMS AOD product, we utilize the GEMS L2 cloud product as a masking criterion. The GEMS L2 cloud product is obtained from the GEMS cloud retrieval algorithm (Kim et al., 2024). With the same hyperspectralmeasurement range, temporal resolution, and spatial resolution as the GEMS aerosol retrieval algorithm (AERAOD), the GEMS cloud retrieval algorithm retrieves the effective cloud fraction (ECF) and provides the cloud radiance fraction (CRF) via the CRF conversion process (Choi et al., 2024).
To filter out pixels biased toward high cloud fraction, we leverage CRF at wavelengths corresponding to the AOD product.
We apply a threestage procedure to calculate the mean field of the L3 AOD products. First, we regrid L2 AOD products using the IDW method with neighboring spatial information to obtain the L3 AOD products. Then, we merge the L3 AOD products by considering spatiotemporal variability in the products, following Kikuchi et al. (2018). Specifically, we merge the L3 AOD products using the previous time T products from the target products of interest. Lastly, we produce the mean field of the L3 AOD products by taking a simple mean of the spatiotemporally merged L3 data. The outline of the method is illustrated in Fig. 1.
3.1 Inverse distance weighting
In this section, we describe the inversedistanceweighting (IDW) algorithm used to obtain L3 AOD products. Several methods, including the nearestneighbor method (Lotrecchiano et al., 2021), the linear interpolation method (Abdullah et al., 2019; Shepard, 1968), and the spline interpolation method (Kuhlmann et al., 2014), have been proposed to interpolate air quality mass. The IDW algorithm (Zimmerman et al., 1999) is one of the most popular linear interpolation methods due to its computational simplicity. Our goal is to obtain the L3 AOD products for each longitude–latitude location.
Let (x_{0},y_{0}) be the target longitude–latitude location for calculating the L3 AOD product. Suppose $({x}_{\mathrm{1}},{y}_{\mathrm{1}}),\mathrm{\dots},({x}_{n},{y}_{n})$ represent the neighboring longitude–latitude locations of (x_{0},y_{0}), endowed with their corresponding L2 AOD products, denoted by AOD$({x}_{\mathrm{1}},{y}_{\mathrm{1}}),\mathrm{\dots},$ AOD(x_{n},y_{n}). To calculate the L3 AOD product in our application, we use grid points within the fourthorder neighbor of (x_{0},y_{0}). This means that, for the given (x_{0},y_{0}), we use the locations that satisfy ${x}_{i}\in ({x}_{\mathrm{0}}r,{x}_{\mathrm{0}}+r)$ and ${y}_{i}\in ({y}_{\mathrm{0}}r,{y}_{\mathrm{0}}+r)$. Specifically, we set the resolution to 0.1° for the East Asia region and 0.05° for the Korean Peninsula region. Then, for the fixed observed time point t_{0}, the IDW estimate is
where the weight of each location is defined as
where d_{i} is the Euclidean distance from (x_{0},y_{0}) to (x_{i},y_{i}) and p is the power parameter. Therefore, Eq. (1) is based on the weighted average of the L2 AOD values from neighboring locations; a larger weight is assigned to grid points close to the location of interest x_{0}.
Depending on the choice of p, the IDW estimates yield different outcomes. As p approaches 0, Eq. (2) indicates equal weight, and the IDW estimate approaches a simple average from neighboring locations. On the other hand, as p approaches ∞, the larger weight is concentrated on the locations near x_{0}, and the IDW estimate converges to the estimate obtained through the nearest neighbor. Although the optimal choice of p can vary depending on the study (Liu et al., 2006), p = 2 is the most commonly used value since when p = 2, the weight for the distance between grid points decays faster to avoid discontinuity (Webster and Oliver, 2007). Following this convention (Isaaks and Srivastava, 1989), we also use p = 2 in our analysis.
Enhanced inverse distance weighting with quality flag information
A quality flag is an indicator that contains data quality information for each grid point. Such indicators are widely used for data cleaning and selection. In our study, we have quality flag information in the L2 AOD products (coded as a 16bit unsigned integer value). The quality flag used in our study is described in Table 2.
To incorporate quality flag information into the analysis, we convert the 16bit integer value into a binary value. For instance, the number 196 can be expressed as 000000011000100, which indicates that the features of bit 2, bit 6, and bit 7 are contained. According to Table 2, pixels with an algorithmic quality flag of 196 are likely to have the features of an AOD value smaller than −0.05 or larger than 3.6. In addition, they may have an SSA value smaller than 0.82 or larger than 1.0 in the presence of clouds, when the solar or viewing zenith angle is above the threshold. Table 3 provides further details about the quality flags.
Then, we define an uncertainty metric u_{i} corresponding to a weight λ_{i} used in the IDW method. The calculation of the uncertainty metric, denoted as u_{i}, is based on a quality flag that is represented by a 16bit unsigned integer. As mentioned, this integer is first converted into a binary format. We then add all the problematic bits with a value of 1 to compute u_{i}. With this quality flag information, the IDW weight used for our method is expressed as
where ${u}_{i}=\sum $(bit values of the quality flag) + 1. Since the high values of u_{i} imply low quality of the data, we take the inverse of u_{i} in the weights of the IDW algorithm. In Eq. (3), q is a power parameter that controls the amount of quality flag information. The larger the q value, the higher the penalty assigned to the grid with a large u_{i} value. In Sect. 4, we observe that incorporating quality flags into the weight can improve the accuracy of the IDW method. Furthermore, we validate the choice of quality flags from the simulation study.
3.2 Spatiotemporal merging algorithm
In this section, we describe a merging algorithm (Kikuchi et al., 2018) that can account for spatiotemporal variability in L3 AOD products. Using spatiotemporal information, we can adjust the weights to produce a more robust and accurate L3 AOD mean field output.
3.2.1 Spatiotemporal variability in AOD_{IDW}
It is crucial to consider the spatiotemporal variability in the IDW estimates when computing the mean field of L3 AOD products (Kikuchi et al., 2018). However, we only have a single IDW estimate, ${\mathbf{AOD}}_{\mathbf{IDW}}(x,y,t)$, at a specific location and time. Since we do not have repeated measures of ${\mathbf{AOD}}_{\mathbf{IDW}}(x,y,t)$, the spatiotemporal variability should be computed using neighboring information. Let $({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ be the location and time of interest and $({x}_{i},{y}_{i},{t}_{i})$ be the neighboring location i. Then, spatiotemporal variability is defined as a rootmeansquare difference (RMSD) of AOD_{IDW} estimates and is expressed as
where N is the number of neighboring pixels within r, which denotes the distance between the point of interest and its neighboring set and previous time T from $({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$. In this work, we consider the fourthorder neighboring grids and a window size of T=3, meaning that r = 0.1, 0.2, 0.3, and 0.4° and that t = 0, 1, 2, and 3.
3.2.2 Hourly merged AOD estimates
Using the spatiotemporal variability discussed in Sect. 3.2.1, we compute hourly combined AOD products. The procedure is summarized as follows. First, we obtain AOD_{pure} by filtering out unreliable grid points. Then, we compute AOD_{merged} at the location of interest by interpolating AOD_{pure}. From this, we can retrieve reliable AOD products.
3.2.3 Computing AOD_{est}
We first introduce ${\mathbf{AOD}}_{\mathbf{est}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$, which is a weighted average of the IDW estimates obtained in Sect. 3.1. The AOD estimate for a target grid $({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ is given as
where
Here, n denotes the number of effective pixels within r and past time T from $({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$, whose AOD_{IDW} values are greater than or equal to 0. In Eq. (5), AOD_{est} is the weighted average of AOD_{IDW}, and weights are defined by the inverse of the spatiotemporal variability discussed in Sect. 3.2.1. Note that the inverse of the variability quantifies the accuracy and reliability of the IDW estimate for each grid; w_{i} indicates the sum of the accuracies over the neighboring region of the target point.
3.2.4 Estimating the error variance
Our goal is to filter out grid points with high variability. Note that the spatiotemporal variability in Eq. (4) increases as the distance between grids increases. Utilizing this relationship, we estimate the spatial and temporal variabilities separately through a regression model. The combined variability, denoted as σ_{0}, for the currently considered grid point $({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ is then calculated as the mean of the spatial and temporal variabilities.
Before estimating the variability, we categorize the values of AOD_{IDW} into different classes. This is because the pattern of spatiotemporal variability varies depending on the magnitude of the AOD values (Kikuchi et al., 2018). Specifically, we categorize AOD_{IDW} values into six bins corresponding to 0.1, 0.25, 0.5, 0.75, 0.9, and 1.0. Note that previous work (Kikuchi et al., 2018) used 12 classes. However, we use six classes because certain classes are rarely observed; using 12 classes can lead to unreliable computation.
Let ${\mathit{\sigma}}_{\mathrm{dist}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ be the spatial variability and ${\mathit{\sigma}}_{\mathrm{time}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ be the temporal variability in the IDW estimate at $({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$. We first compute the average of the spatiotemporal variability ${\mathit{\sigma}}_{\mathrm{IDW}}({x}_{i},{y}_{i},{t}_{i})$ based on r and the class. Then, we regress the averaged values obtained for each component of the vector r = (0.1°, 0.2°, 0.3°, 0.4°) on a design matrix $[\mathbf{1},\mathit{r},{\mathit{r}}^{\mathbf{2}}]$, which is a secondorder design matrix for r pertaining to each class. Lastly, we obtain the spatial variability ${\mathit{\sigma}}_{\mathrm{dist}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ from the intercept estimate of the quadratic regression model fitting. We can obtain the temporal variability ${\mathit{\sigma}}_{\mathrm{time}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ in a similar manner. We first compute the average of the spatiotemporal variability ${\mathit{\sigma}}_{\mathrm{IDW}}({x}_{i},{y}_{i},{t}_{i})$ based on the time point and class. We regress them on a design matrix $[\mathbf{1},\mathit{t},{\mathit{t}}^{\mathbf{2}}]$ for each class and obtain ${\mathit{\sigma}}_{\mathrm{time}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ from the intercept estimate of the quadratic regression model fitting.
Finally, we compute the error variance by taking the average of σ_{dist} and σ_{time}, given as
The calculated error variance contains a measurement error caused by sensor noise that varies over time and space.
3.2.5 Computing AOD_{pure}
Here, we obtain AOD_{pure} by filtering out uncertain AOD_{IDW} values. For this, we introduce the estimated error in AOD_{pure}, which is defined as
Here, ${\mathit{\sigma}}_{\mathrm{est}}^{\mathrm{2}}$ is the error variance of AOD_{est} from Eq. (6) and ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ is the combined error variance from Eq. (7). To filter out uncertain AOD_{IDW} values, we consider an upper threshold of AOD_{IDW}, expressed as
Assuming a Gaussian distribution, if ${\mathbf{AOD}}_{\mathbf{IDW}}({x}_{\mathrm{0}},{y}_{\mathrm{0}},{t}_{\mathrm{0}})$ exceeds the upper threshold of the 99 % confidence interval, we consider the value to be unreliable and thus exclude it from the mean field calculation.
3.2.6 Computing AOD_{merged}
Based on Sect. 3.2.5, we calculate the total variability within the neighboring grids of the target grid. Then, we use the ratio of the inverse of this total variability to assign weights for calculating a weighted average of AOD_{pure}, resulting in AOD_{merged}, given as follows:
where
This merging process not only utilizes the reliable value (AOD_{pure}) but also incorporates the reliability (σ_{pure}) as a weight, resulting in a more trustworthy gapfilling outcome. In fact, in the study by Kikuchi et al. (2018), the root mean squared error (RMSE) of AOD_{merged} was notably lower (0.11) compared to the RMSE of AOD_{IDW} (0.20).
In this section, we conduct a simulation study to validate the choice of quality flags. The data generation procedure is summarized as follows.
4.1 Generating simulation data

For our simulation study, we constructed a 70 × 70 lattice over a 1 × 1 square domain with a grid spacing of 0.1°. Each grid within this lattice represents a location in our simulated dataset. In total, the lattice comprises 4900 grid cells, covering an area of 7° × 7°. Consequently, each unit grid on the lattice represents an area of 0.1° × 0.1°. We generate each element of $\mathbf{X}\in {\mathbb{R}}^{\mathrm{4900}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}\mathrm{2}}$ from Unif(0,1), i.e., uniform distribution with support from 0 to 1. This means that we have 4900 locations, each containing two pieces of coordinate information: longitude and latitude. We use ${\mathit{\beta}}_{\mathrm{1}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}\mathrm{2}}=(\mathrm{1},\mathrm{1})$ as a true coefficient vector.

For each location, we simulate a zeromean Gaussian process W from $N(\mathrm{0},(\mathit{\tau}{\mathbf{M}}^{\prime}\mathbf{\text{QM}}{)}^{\mathrm{1}})$, where M is obtained by taking the first k eigenvectors of the Moran operator (Hughes and Haran, 2013) using smoothness parameter τ. Here, Q = diag(A1) − A is a precision matrix calculated from the adjacency matrix A. Note that 1 is an allones vector.

We simulate AOD datasets from $\mathbf{Y}=\mathbf{X}\mathit{\beta}+\mathbf{W}\in {\mathbb{R}}^{\mathrm{4900}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}\mathrm{2}}$. In our simulation, Xβ represents the fixed effect, while W∈ℝ^{4900} accounts for spatial correlation in AOD products.

To generate missing data for the simulated Y that resemble the GEMS L2 product, we apply the observed missing pattern from the GEMS AOD data to the simulated Y from step 3. Specifically, we apply the missing data pattern observed in the GEMS AOD data for the 7° × 7° region in East Asia, with a grid spacing of 0.1°, collected on 1 April at 04:45 UTC. This selected dataset aligns with the 70 × 70 lattice in our simulation and has approximately 20 % of its values missing, effectively replicating the realistic missing data characteristics found in the actual GEMS AOD observations.

For a realistic simulation that incorporates the physical implications of quality flags, we adapt the observed quality flags from the GEMS AOD data to the simulated Y from step 3. Similar to step 4, we use the quality flag from the same data (i.e., data from the 7° × 7° region in East Asia (with a grid spacing of 0.1°), collected on 1 April at 04:45 UTC) for the simulated Y.
We repeat steps 1–3 100 times to generate different realizations of Y. Then, for each simulated Y, we apply the identical missing patterns (step 4) and quality flags (step 5) obtained from the real dataset. Figure 2 illustrates an example of a simulated dataset.
4.2 Sensitivity analysis of q
Here, we investigate the performance of the IDW method with regard to GEMS data by varying q in Eq. (3). Specifically, we consider q values of 0.5, 1, and 1.5 in our experiment. We first examine whether there is a significant difference in IDW estimates with difference choices of q. Figure 3 indicates that the IDW estimates are comparable with various q values. Table 4 also shows that the summary statistics of AOD_{IDW} are quantitatively similar across different q choices. Therefore, we conclude that the IDW algorithm is robust regardless of the q value used. To simplify the calculation, we set q = 1 in our analysis.
4.3 Quality flag simulation
As explained in Sect. 4.2 and Eq. (3), quality flag indicators weigh the uncertainty in the IDW algorithm. To improve the accuracy of the IDW algorithm, it is necessary to find the optimal bit combination of the quality flag by performing simulation studies. Therefore, we first identify bits that show substantially lower mean squared errors (MSEs) and then combine these specific bits into groups. Here, the MSE is calculated between the simulated AOD dataset and the IDW results, given by
where AOD_{i} is the simulated AOD value at location i and IDW_{i} is the IDW result at location i. Note that we refer to these bit combinations of quality flags as “cases”. For example, we may find that bits 0, 1, and 2 have a significantly low MSE. This means that we can create various combinations of these bits, such as 0 ⋅ 1, 0 ⋅ 2, or even 0 ⋅ 1 ⋅ 2, and each combination can be denoted as a particular “case”.
Before finding the optimal case of the quality flag through a simulation study, we examine whether the result differs depending on the value of r in the IDW algorithm. Note that this simulation study defines the unit of r as one unit grid and not as the distance based on coordinates, as described in Sect. 3.1. After careful consideration and analysis, we have chosen to set the neighboring order to 3 for our study because varying the order of the neighbors did not yield any significant differences in the mean squared error (MSE). However, to provide a comprehensive understanding of the impact of the order on the interpolation process, we have included additional IDW results with increasing order sizes in Appendix A. These supplementary results demonstrate that larger order sizes can potentially lead to oversmoothing. Figure 4 shows the generated AOD simulation data, with panel (b) indicating the data incorporating the missing value pattern and panel (b) visualizing the result of applying the IDW algorithm to the data.
We then follow the following procedure. After selecting specific bits with a lower MSE than bits fixed with an order of 3, we repeat the experiment to account for diverse uncertainty term calculations by considering various cases in the IDW algorithm. We then find the optimal case that obtains the highest accuracy by comparing the accuracies among the cases. Here, MSE values are evaluated between the simulation and the imputed data, treating the simulated data as real data.
To select the quality flag bit for making the combination case, we first calculate the MSE for each quality flag bit. Figure 5 displays the MSE values as boxplots for each bit. Although the boxplots indicate the results for bits 0 through 15, bit 2 and bit 6 exhibit significantly lower MSEs than other bits. The medians for the two bits are 0.122 and 0.123, respectively. Bit 2 indicates whether SSA or AOD is outside a specific value range, while bit 6 shows the presence of clouds. We then compile six different candidate cases, including bit 2 and bit 6. Table 5 summarizes these six candidate cases.
To find the optimal quality control case that yields the lowest MSE among the six candidate cases, we calculate the MSE for each case using two different τ values, τ = 1 and τ = 6, and present boxplots in Fig. 6. In Fig. 6a and b, each case shows a different MSE, while case 5, which contains bits 0, 2, and 6, shows the lowest median (0.387) in Fig. 6a. Therefore, for our real application, we use the quality flag for case 5 (bits 0, 2, and 6) when calculating the uncertainty term for the IDW algorithm.
We examine whether the optimal case for the IDW algorithm changes depending on the smoothing parameter τ of the simulation data. Figure 6a and b show the average MSE, represented as boxplots, with different smoothing parameters. We find that while the overall magnitude of the MSE value differs, the hierarchy of the MSEs across the cases remains unchanged. This suggests that the algorithm's performance is robust regardless of the smoothing parameter used.
5.1 Spatial resolution of output product
A geostationary satellite instrument, such as the GEMS, observes a fixed position. Therefore, if we use grids that are too small, the output will have many missing values; an appropriate grid size should be selected to ensure adequate coverage of the mean fields. Furthermore, we need to consider the effective range, accuracy, and computation time needed to obtain mean fields when choosing a grid size. If we use too fine a spatial resolution, computational costs will exponentially increase. On the other hand, grid sizes that are too large can lead to inaccurate output results. Accordingly, in our study, we set the mean field spatial resolution of the longitude–latitude grid to 0.1° × 0.1° and 0.05° × 0.05° for the East Asia and Korean Peninsula regions, respectively.
5.2 Level2 aerosol optical depth data
As evidenced by various studies (Kaufman et al., 2005; Loeb and ManaloSmith, 2005; Matheson et al., 2005), AOD exhibits a positive correlation with cloud fraction (CF), implying that proximity to clouds can result in a statistical increase in AOD measurements. Additionally, a dependency of AOD on solar and viewing zenith angles was observed (de Miguel et al., 2011), highlighting the complexities involved in accurate AOD estimation under varied atmospheric conditions. To address these complexities, this study involved masking data based on three variables, given below. Cloud radiance fraction is incorporated within the GEMS L2 cloud product, while the others are incorporated within the GEMS L2 AOD product. These variables are meticulously designed to account for variability and uncertainties in satellite data (Choi et al., 2024), addressing factors such as solar and viewing zenith angles and cloud radiance fraction. This methodology ensures a more precise and reliable estimation of AOD, which is crucial for understanding atmospheric dynamics and environmental monitoring. Unreliable values are treated as missing values when we apply the IDW algorithm. The thresholds are as follows:

cloud radiance fraction ≥ 0.4

solar zenith angle ≥ 70°

viewing zenith angle ≥ 70°.
In our study, we set two spatial domains for the mean field outputs – one with latitude (30–43° N) and longitude (123–131° E) corresponding to the vicinity of the Korean Peninsula and the other with latitude (32–43° N) and longitude (115–131° E) including the vicinity of the Shandong Peninsula in China.
5.3 Mean field estimates of the GEMS L3 AOD product
Figures 7 and 8 compare the spatiotemporally merged products with the simply averaged products. For spatiotemporally merged L3 AOD products, we apply the procedure described in Sect. 3. On the other hand, we take the mean of the IDW estimates for simply averaged L3 AOD products. As mentioned in Kikuchi et al. (2018), grid points with AOD values exceeding 1.0 are extremely rare. Therefore, to focus on the majority of values for detailed characterization, we set the threshold to 1.0. In Fig. 7, we observe that more missing values occur in the area of latitude (20–30° N) and longitude (130–140° E) for the spatiotemporally merged products compared to that for the simply averaged products. This can also be verified by the ratio of missing values for each product shown in each figure. This is due to the fact that a spatiotemporal merging procedure only considers reliable AOD estimates, while a simple averaging method does not. Therefore, the simply averaged products can be regarded as more unreliable, although they have fewer missing values. In Fig. 7, we also observe that the spatiotemporal merging method can provide smoother mean field outputs; for example, there is a significant difference in the area of latitude (35–45° N) and longitude (125–135° E). Similar trends are also observed in Fig. 8. Mean field estimates at all three wavelengths for the East Asia and the Korean Peninsula regions are provided in Appendix B and C, respectively.
As described in Sect. 3, we obtain AOD_{pure} from AOD_{est} by considering both the spatiotemporal variability σ_{0} and the estimation uncertainty σ_{est}. This procedure allows us to obtain more reliable AOD estimates. Furthermore, we can use more robust AOD estimates when computing the mean field, resulting in a smoother output.
5.4 Qualitative evaluation of the mean field products
Directly evaluating the accuracy of the L3 AOD mean field products is challenging because there are no true values for these products. Therefore, we compare our results with previous studies using qualitative aspects.
In Fig. 9a, we observe that our mean field products obtained at a 550 nm wavelength on 1 April are similar to the springtime global distribution of AOD at the same wavelength. Furthermore, high values due to dust are observed near the Taklamakan Desert, and high values due to biomass burning are observed in Southeast Asia during spring. This indicates that the computed mean fields can effectively capture realworld phenomena. In addition, we observe that the overall trends in AOD values are similar to those from Moderate Resolution Imaging Spectroradiometer (MODIS) data (Tian et al., 2023), although there are some discrepancies among values pertaining to the vicinity of the Taklamakan Desert and certain areas in Southeast Asia.
5.5 Quantitative evaluation of the mean field products
To measure the smoothness of the mean field products, we compute the secondorder gradients of the AOD estimates from the longitude and latitude. Specifically, we compute gradients at each pixel point based on its neighboring AOD values, following previous studies (Fornberg, 1988; Quarteroni et al., 2007; Durran, 1999). Then, we compare the absolute mean of the gradients from the simple and spatiotemporal averaging methods. Table 6 indicates that gradients in both directions are smaller in the spatiotemporal merging method. Considering that AOD exhibits significant spatiotemporal variability, our method provides more realistic mean field products than the simple averaging method.
In remote sensing, data are processed according to different levels. For example, the L2 dataset, which contains atmosphericstate retrievals, is converted into the L3 dataset. Specifically, we focus on AOD retrievals from the GEMS instrument through gap filling and noise filtering. We improve the quality of L3 AOD mean field products by considering quality flag information and spatiotemporal variability (Kikuchi et al., 2018). Specifically, the contribution of our work is summarized as follows.
First, we improve the performance of the IDW algorithm by including quality flag information in the weight calculation. We assign weights that are inversely proportional to the number of poorquality indicators. To validate the choice of different quality flags, we conduct simulation studies. We observe that including bits 0, 2, and 6 from the quality flags significantly improves the accuracy of the IDW algorithm. We apply this novel approach to GEMS AOD data covering various regions and wavelengths.
Second, we apply a spatiotemporal merging method (Kikuchi et al., 2018) to GEMS AOD data. Unlike previous work (Park et al., 2023) that only considered spatial variability, our method also accounts for temporal variability from previous time points. We observe that our mean field products exhibit similar trends to the previous studies, indicating that the products are reliable.
Although our current study has made notable progress in enhancing the accuracy of AOD mean field estimation, several avenues for future research remain. One potential direction involves integrating additional data, such as cloud information and distinctions between oceanic and terrestrial regions, to further refine our results by considering the impact of cloud cover on aerosol retrievals. Our method can account for variability due to cloud contamination by utilizing quality flag information in the IDW estimates. Note that we cannot use physical mechanisms (e.g., aerosols produced by wildfires) in the interpolation step due to limited data sources. Developing extensions of our approach by incorporating physical mechanisms may provide interesting future research avenues. Additionally, with the available information on quality flags, one direction for future work could involve developing methodologies for adaptively weighting or selecting quality flag bits by employing statistical variable selection methodologies. Validating our AOD mean field products against groundbased measurements or other satellite datasets could also offer valuable insights into their reliability and consistency, thereby helping to identify any potential biases or uncertainties. Lastly, conducting a sensitivity analysis to choose hyperparameters (e.g., the order of neighboring grids and time windows) would be useful for improving the performance of the method.
We compare the results of applying the IDW algorithm with varying orders of neighboring grids used for a weighted average. Observing the visualization, we find a remarkable difference in the imputed area as fewer missing values remain with an order of 9 compared to an order of 3. However, the MSE values are 0.255 and 0.25, respectively, showing no significant difference in numbers. Despite the negligible difference between the two window sizes, we determined that the window size should be 3 since an order of 9 could cause oversmoothing.
The GEMS Level2 products are available at https://nesc.nier.go.kr/en/html/datasvc/data.do?pageIndex=1&outputInnb=41&atrb=AERAOD443 (National Institute of Environmental Research Environmental Satellite Center, 2023a), which provides direct access to the latest 7 d observation data. For access to previous observations, readers may wish to use the application programming interface (API), available at https://nesc.nier.go.kr/ko/html/index.do (National Institute of Environmental Research Environmental Satellite Center, 2023b).
GEMS data were provided by YesC, JK, and WL. All authors participated in developing the overall algorithm design, which was based on a conceptual idea of a spatiotemporal merging algorithm proposed by SK. Specifically, the spatiotemporal merging algorithms were developed by SK and SP with the guidance of JP and IHJ. Simulation studies of quality flags incorporated into the spatiotemporal merging algorithms were conducted by SL and DO. Data analyses were performed by YeoC and HK. The paper was written, edited, and proofread by all the authors.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This article is part of the special issue “GEMS: first year in operation (AMT/ACP interjournal SI)”. It is not associated with a conference.
The authors are grateful to the reviewers for their careful reading and valuable comments.
This work was supported by the National Research Foundation of Korea (grant nos. 2020R1C1C1A0100386814, RS202300217705, and RS202300218377) and the ICAN (ICT Challenge and Advanced Network of HRD) support program (grant no. RS202300259934), supervised by the IITP (Institute of Information & Communications Technology Planning & Evaluation).
This paper was edited by Rokjin Park and reviewed by Won Chang and one anonymous referee.
Abdullah, S., Ismail, M., Ali Najah Ahmed, A.M., and Abdullah, A.: Forecasting Particulate Matter Concentration Using Linear and NonLinear Approaches for Air Quality Decision Support, Atmosphere, 10, 667, https://doi.org/10.3390/atmos10110667, 2019. a
Alexander, L., Allen, S., Bindoff, N., Breon, F.M., Church, J., Cubasch, U., Emori, S., Forster, P., Friedlingstein, P., Gillett, N., Gregory, J., Hartmann, D., Jansen, E., Kirtman, B., Knutti, R., Kanikicharla, K., Lemke, P., Marotzke, J., MassonDelmotte, V., and Xie, S.P.: IPCC, 2013: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp., ISBN 9781107057991, ISBN 9781107661820, 2013. a
Brauer, M., Freedman, G., Frostad, J., Donkelaar, A., Martin, R., Dentener, F., Van Dingenen, R., Estep, K., Amini, H., Apte, J., Balakrishnan, K., Barregard, L., Broday, D., Feigin, V., Ghosh, S., Hopke, P., Knibbs, L., Kokubo, Y., Liu, Y., and Cohen, A.: Ambient Air Pollution Exposure Estimation for the Global Burden of Disease 2013, Environ. Sci. Technol., 50, 79–88, https://doi.org/10.1021/acs.est.5b03709, 2015. a
Charlson, R.: Atmospheric visibility related to aerosol mass concentration – A review, Environ. Sci. Technol., 3, 913–918, https://doi.org/10.1021/es60033a002, 1969. a
Charlson, R. J., Schwartz, S. E., Hales, J. M., Cess, R. D., Coakley, J. A., Hansen, J. E., and Hofmann, D. J.: Climate Forcing by Anthropogenic Aerosols, Science, 255, 423–430, https://doi.org/10.1126/science.255.5043.423, 1992. a
Chen, C., Dubovik, O., Schuster, G., Chin, M., Henze, D., Lapyonok, T., Derimian, Y., and Ying, Z.: Multiangular polarimetric remote sensing to pinpoint global aerosol absorption and direct radiative forcing, Nat. Commun., 13, 7459, https://doi.org/10.1038/s4146702235147y, 2022. a
Cheng, M. and Tsai, Y. I.: Characterization of visibility and atmospheric aerosols in urban, suburban, and remote areas, Sci. Total Environ., 263, 101–114, https://doi.org/10.1016/S00489697(00)006707, 2001. a
Cho, Y., Kim, J., Go, S., Kim, M., Lee, S., Kim, M., Chong, H., Lee, W.J., Lee, D.W., Torres, O., and Park, S. S.: First atmospheric aerosolmonitoring results from the Geostationary Environment Monitoring Spectrometer (GEMS) over Asia, Atmos. Meas. Tech., 17, 4369–4390, https://doi.org/10.5194/amt1743692024, 2024. a
Choi, H., Liu, X., Jeong, U., Chong, H., Kim, J., Ahn, M. H., Ko, D. H., Lee, D.W., Moon, K.J., and Lee, K.M.: Geostationary Environment Monitoring Spectrometer (GEMS) polarization characteristics and correction algorithm, Atmos. Meas. Tech., 17, 145–164, https://doi.org/10.5194/amt171452024, 2024. a, b
Cressie, N.: Mission CO_{2}ntrol: A Statistical Scientist's Role in Remote Sensing of Atmospheric Carbon Dioxide, J. Am. Stat. Assoc., 113, 152–168, https://doi.org/10.1080/01621459.2017.1419136, 2018. a
de Miguel, A., Mateos, D., Bilbao, J., and Román, R.: Sensitivity analysis of ratio between ultraviolet and total shortwave solar radiation to cloudiness, ozone, aerosols and precipitable water, Atmos. Res., 102, 136–144, 2011. a
Durran, D. R. (Ed.): Numerical methods for wave equations in geophysical fluid dynamics, in: Texts in applied mathematics, Springer, ISBN 0387983767, 1999. a
Fornberg, B.: Generation of Finite Difference Formulas on Arbitrarily Spaced Grids, Math. Comput., 51, 699–706, https://doi.org/10.2307/2008770, 1988. a
Go, S., Kim, J., Mok, J., Irie, H., Yoon, J., Torres, O., Krotkov, N., Labow, G., Kim, M., Koo, J.H., Choi, M., and Lim, H.: Groundbased retrievals of aerosol column absorption in the UV spectral region and their implications for GEMS measurements, Remote Sens. Environ., 245, 111759, https://doi.org/10.1016/j.rse.2020.111759, 2020a. a
Go, S., Kim, J., Park, S. S., Kim, M., Lim, H., Kim, J.Y., Lee, D.W., and Im, J.: Synergistic Use of Hyperspectral UVVisible OMI and Broadband Meteorological Imager MODIS Data for a Merged Aerosol Product, Remote Sens., 12, 3987, https://doi.org/10.3390/rs12233987, 2020b. a, b
Hughes, J. and Haran, M.: Dimension reduction and alleviation of confounding for spatial generalized linear mixed models, J. Roy. Stat. Soc. B, 75, 139–159, 2013. a
Isaaks, E. H. and Srivastava, R. M.: An Introduction to Applied Geostatistics, Oxford University Press, New York, 561 pp., ISBN10: 0195050134, ISBN13: 9780195050134, 1989. a
Kaufman, Y., Tanré, D., and Boucher, O.: A satellite view of aerosols in the climate system, Nature, 419, 215–23, https://doi.org/10.1038/nature01091, 2002. a
Kaufman, Y., Koren, I., Remer, L., Tanré, D., Ginoux, P., and Fan, S.: Dust transport and deposition observed from the TerraModerate Resolution Imaging Spectroradiometer (MODIS) spacecraft over the Atlantic Ocean, J. Geophys. Res., 110, D10S12, https://doi.org/10.1029/2003JD004436, 2005. a
Kikuchi, M., Murakami, H., Suzuki, K., Nagao, T. M., and Higurashi, A.: Improved hourly estimates of aerosol optical thickness using spatiotemporal variability derived from Himawari8 geostationary satellite, IEEE T. Geosci. Remote, 56, 3442–3455, 2018. a, b, c, d, e, f, g, h, i, j
Kim, B.R., Kim, G., Cho, M., Choi, Y.S., and Kim, J.: First results of cloud retrieval from the Geostationary Environmental Monitoring Spectrometer, Atmos. Meas. Tech., 17, 453–470, https://doi.org/10.5194/amt174532024, 2024. a
Kim, M., Kim, J., Torres, O., Ahn, C., Kim, W., Jeong, U., Go, S., Liu, X., Moon, K., and Kim, D.R.: Optimal EstimationBased Algorithm to Retrieve Aerosol Optical Properties for GEMS Measurements over Asia, Remote Sens.Basel, 10, 162, https://doi.org/10.3390/rs10020162, 2018. a, b
Kuhlmann, G., Hartl, A., Cheung, H. M., Lam, Y. F., and Wenig, M. O.: A novel gridding algorithm to create regional trace gas maps from satellite observations, Atmos. Meas. Tech., 7, 451–467, https://doi.org/10.5194/amt74512014, 2014. a
Liu, X., Chance, K., Sioris, C., Kurosu, T., Spurr, R., Martin, R., Fu, T.M., Logan, J., Jacob, D., Palmer, P., Newchurch, M., Megretskaia, I., and Chatfield, R.: First Directly Retrieved Global Distribution of Tropospheric Column Ozone from GOME: Comparison with the GEOSCHEM Model, J. Geophys. Res., 111, D02308, https://doi.org/10.1029/2005JD006564, 2006. a
Loeb, N. G. and ManaloSmith, N.: Topofatmosphere direct radiative effect of aerosols over global oceans from merged CERES and MODIS observations, J. Climate, 18, 3506–3526, 2005. a
Lotrecchiano, N., Sofia, D., Giuliano, A., Barletta, D., and Poletto, M.: Spatial Interpolation Techniques For innovative Air Quality Monitoring Systems, Chem. Engineer. Trans., 86, 2021, https://doi.org/10.3303/CET2186066, 2021. a
Matheson, M. A., Coakley Jr., J. A., and Tahnk, W. R.: Aerosol and cloud property relationships for summertime stratiform clouds in the northeastern Atlantic from Advanced Very High Resolution Radiometer observations, J. Geophys. Res., 110, D24204, https://doi.org/10.1029/2005JD006165, 2005. a
National Institute of Environmental Research Environmental Satellite Center: Data Services, https://nesc.nier.go.kr/en/html/datasvc/data.do?pageIndex=1&outputInnb=41&atrb=AERAOD443 (last access: 29 August 2024), 2023a. a
National Institute of Environmental Research Environmental Satellite Center: GEMS OpenAPI, https://nesc.nier.go.kr/ko/html/index.do (last access: 29 August 2024), 2023b. a
Park, J., Kang, S., Jin, I. H., Kim, S., and Oh, D.: Geostationary Environment Monitoring Spectrometer (GEMS) User Manual Level S3 Mean Field Algorithm, National Institute of Environmental Research, Republic of Korea, https://nesc.nier.go.kr/ko/html/satellite/guide/guide.do (last access: 25 August 2024), 2023. a, b, c
Park, S. S., Kim, S.W., Song, C.K., Park, J.U., and Bae, K.H.: Spatiotemporal variability of aerosol optical depth, total ozone and NO_{2} over East Asia: Strategy for the validation to the GEMS Scientific Products, Remote SensBasel, 12, 2256, https://doi.org/10.3390/rs12142256, 2020. a
Pöschl, U.: Atmospheric Aerosols: Composition, Transformation, Climate and Health Effects, Angew. Chem. Int. Edit., 44, 75207540, https://doi.org/10.1002/anie.200501122, 2006. a
Quarteroni, A., Sacco, R., and Saleri, F.: Numerical Mathematics, Springer, 37, 442–577, ISBN 9781475773941, https://doi.org/10.1007/b98885, 2007. a
Shepard, D.: A TwoDimensional Interpolation Function for IrregularlySpaced Data, ACM National Conference, 23, 517–524, https://doi.org/10.1145/800186.810616, 1968. a
Stocker, T.: IPCC, 2013: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp., ISBN 9781107057991, ISBN 9781107661820, 2013. a
Tager, I. B.: Health effects of aerosols: Mechanisms and epidemiology, in: Aerosols Handbook, CRC Press, 635–712, ISBN: 9780429209222, 2004. a
Tian, X., Tang, C., Wu, X., Yang, J., Zhao, F., and Liu, D.: The global spatialtemporal distribution and EOF analysis of AOD based on MODIS data during 2003–2021, Atmos. Environ., 302, 119722, https://doi.org/10.1016/j.atmosenv.2023.119722, 2023. a
Webster, R. and Oliver, M. A.: Geostatistics for environmental scientists, John Wiley & Sons, https://doi.org/10.1002/9780470517277, 2007. a
Zhang, F.: Factors Influencing the Spatio–Temporal Variability of Aerosol Optical Depth over the Arid Region of Northwest China, Atmosphere, 15, 54, https://doi.org/10.3390/atmos15010054, 2023. a
Zimmerman, D., Pavlik, C., Ruggles, A., and Armstrong, M. P.: An experimental comparison of ordinary and universal kriging and inverse distance weighting, Math. Geol., 31, 375–390, 1999. a
 Abstract
 Introduction
 GEMS data
 Methodology
 Simulation
 GEMS data application
 Conclusions
 Appendix A: Sensitivity analysis for the radius of the IDW algorithm
 Appendix B: Mean field estimates of the GEMS L3 AOD product for East Asia
 Appendix C: Mean field estimates of the GEMS L3 AOD product for the Korean Peninsula region
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 GEMS data
 Methodology
 Simulation
 GEMS data application
 Conclusions
 Appendix A: Sensitivity analysis for the radius of the IDW algorithm
 Appendix B: Mean field estimates of the GEMS L3 AOD product for East Asia
 Appendix C: Mean field estimates of the GEMS L3 AOD product for the Korean Peninsula region
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References