A Geostatistical Approach to Spatial Quality Assessment of Coarse Spatial Resolution Remote Sensing Products

A geostatistical framework for spatial quality assessment framework of coarse resolution remote sensing products is presented that can account for either the scale difference or the uncertainty of reference value prediction at coarse resolutions. A set of multiple reference field realizations is first generated at a fine spatial resolution using geostatistical simulation to explore the uncertainty in the true unknown reference field. The upscaling of multiple reference field realizations to coarse resolution is then followed to match the spatial resolution of the target remote sensing product and create coarse resolution reference fields. The simulated reference values at each coarse pixel are compared to the corresponding reported value from the coarse resolution remote sensing product, yielding alternative error values, from which several location-dependent statistics such as mean error, mean absolute error, and probability of overestimation can be computed. An experiment involving monthly Tropical Rainfall Measuring Mission (TRMM) precipitation products and point-level rain gauge data over South Korea illustrates the applicability of the proposed approach. The spatially distributed error statistics are useful to identify areas with larger errors and the degree of overestimation in the study area, leading to the identification of areas with unreliable estimates within the TRMM precipitation products. Therefore, it is expected that the geostatistical assessment framework presented in this paper can be effectively used to evaluate the spatial quality of coarse resolution remote sensing products.


Introduction
Many environmental parameters such as precipitation, temperature, soil moisture, vegetation index, and land-cover are critical for monitoring the Earth's environment and are routinely obtained from remote sensing images and products [1][2][3][4][5][6]. Advantages of remote sensing products over field surveys are their ability of providing periodic and exhaustive thematic information on the Earth's environment. The spatial scale of the monitoring effort can be global, regional, or local, depending on the spatial resolution of the remote sensing products. Different periodic monitoring or time-series analysis are also possible depending on temporal resolutions [7,8].
Any products derived from remote sensing data, however, inevitably include errors. Particularly, global and regional satellite-based products are generated by tuning with sparse local data; thus the reliability may vary across regions. In turn, the error of the satellite-based products affects the environmental model predictions and leads to error propagation problems [19]. Thus, the quality of such products should be assessed prior to their utilization in environmental modeling.
Traditional error assessment of remote sensing products is typically based on a direct comparison against a small number of point-level data. Such a direct comparison, however, does not account for the spatial resolution difference between the point-based ground observations and the coarse resolution pixels. Most satellite-based products for global and regional monitoring have low spatial but high temporal resolution. Thus, the spatial resolution difference should be properly accounted for during the comparison. Furthermore, error statistics computed from this comparison, such as mean absolute errors and overall accuracy, provide only locationinvariant errors but no information on their spatial distribution [20]. Such information on spatial error distribution can be used for further field surveys and detailed analysis of local errors.
An alternative error assessment route involves the construction of a fine spatial resolution reference product from validation data and the subsequent comparison of that reference product with the coarse resolution remote sensing products [21][22][23][24]. The simplest way to construct the reference product is to interpolate or average the values from validation data within each coarse pixel. However, the quality of this approach depends heavily on the number and location of the validation data. Furthermore, it is not possible to generate an error-free reference product, even in cases where numerous validation data are available and advanced interpolation methods are applied. Thus, it is necessary to consider the uncertainty of the reference product, due to lack of knowledge on the spatial distribution of the parameter at the fine spatial resolution.
To overcome these limitations, this paper presents a geostatistical framework for spatial quality assessment of coarse resolution satellite-based products using point-level validation data and geostatistical simulation. More precisely, multiple fine spatial resolution reference maps are first generated using geostatistical simulation, instead of a deterministic single reference product with uncertainty. The simulated maps represent possible versions of the unknown reference product, hence accounting for its inherent uncertainty [25,26]. These alternative synthetic realizations are then upscaled at the coarse spatial resolution, and multiple error values are computed at each coarse pixel through a pixel-by-pixel product comparison. From the set of alternative error values at each coarse pixel, several useful location-dependent error statistics can be computed, such as the probability that attribute values in coarse resolution satellite-based products overestimate the reference values. The applicability of the proposed approach is demonstrated through an experiment with Tropical Rainfall Measuring Mission (TRMM) precipitation products and validation data acquired over South Korea.

Study Area and
Data. An experiment is conducted using coarse resolution TRMM monthly precipitation products and rain gauge data acquired over South Korea. Three TRMM 3B43 monthly precipitation rates (mm/h) at a spatial resolution of 0.25 ∘ , acquired in three representative seasons of South Korea in 2013 (i.e., April, July, and October for spring, summer, and autumn, respectively), were used as the coarse resolution remote sensing products. The winter season data were not considered in this study because heavy snow is common in Korean winters, and satellite-based products mainly provide rainfall information, not precipitation including rainfall and snowfall. It should be noted that the main purpose of this experiment is to demonstrate the procedures and potential of the geostatistical framework proposed in this paper, not to reveal detailed local characteristics of monthly precipitation in the study area.
The original monthly TRMM data at a spatial resolution of 0.25 ∘ were first reprojected using a Transverse Mercator projection with a spatial resolution of 25 km. The monthly precipitation rates were then converted to monthly accumulated precipitation at a mm scale for quantitative comparison with rain gauge data ( Figure 1). The rainy season in South Korea starts in late June and continues until late July [8]. Therefore, the precipitation amount of July is relatively greater than in other seasons.
Monthly accumulated rainfall measurements obtained from 606 automatic weather stations (AWS) over South Korea were used to evaluate the TRMM precipitation products. Since the AWS data are available mainly over land, not all blocks from the TRMM products over South Korea were considered for spatial quality assessment and comparison with AWS data, only the 25 km overland blocks were assessed (gray rectangles with blue boundary in Figure 2). Figure 3 shows all the process steps applied for the geostatistical spatial quality assessment framework proposed in this paper: (1) generation of multiple fine resolution reference maps using geostatistical simulation, (2) upscaling of the fine resolution reference maps to the resolution of the coarse scale remote sensing product, (3) generation of multiple error maps through a pixel-by-pixel comparison of the upscaled reference maps and the coarse resolution remote sensing product, and (4) derivation of the summary statistics from multiple error maps. The key idea of the proposed approach is to generate multiple reference maps from a small number of point-level reference data and then to compare them with the coarse resolution remote sensing products to account for the uncertainty attached to the reference map generation.

Geostatistical Spatial Quality Assessment Framework.
First, a fine spatial resolution reference set from validation data is generated using geostatistical simulation. Unlike kriging that provides one single optimal map in a leastsquares sense [25], geostatistical simulation aims at generating alternative equi-probable realizations of an unknown truth, conditioning all available data and preserving statistical properties of the input data (e.g., histogram and variogram) [26,27]. The AWS data used in this study did not show any strong skewed distributions, thus sequential Gaussian simulation (SGSIM) was applied to generate multiple realizations of the reference precipitation. The simulation was conducted at a spatial resolution of 1 km, considering the density and location of the AWS data.
Suppose there is a set of n precipitation measurements { (u ), = 1, ⋅ ⋅ ⋅ , } from AWS in the study area. Within the multi-Gaussian framework, SGSIM proceeds with the following steps [25]: (1) If necessary, the AWS data are first transformed to follow a standard Gaussian distribution using normal score transform that is a rank-preserving and invertible quantile transform [26]. The normal score transformed dataset { G (u ), = 1, ⋅ ⋅ ⋅ , } is used as an input for the SGSIM.  (2) A random path is defined to visit each simulation pixel only once.  Figure 2: Location of AWS data in the study area (black dots) and zones for error analysis. Black rectangles and polylines represent TRMM 25 km pixels and administrative boundaries, respectively. Gray rectangles with blue boundary denote TRMM pixels used for quality assessment. Step 2 Step 3 Step 4 Step 1 for the Gaussian ccdf, respectively. At the first pixel of the random path, the conditioning data are only the normal score transformed AWS dataset. Otherwise, the previously simulated values at pixels visited before are included into the conditioning information for the simulation of the other pixel.
(4) A simulation value is drawn from the Gaussian ccdf using a Monte-Carlo simulation.
(5) Until all pixels in the study area are visited, the above two steps are repeated at each pixel.
(6) The normal score back-transform is applied to compute simulated values in the original attribute space.
(7) A new realization is generated by repeating all the above steps with a different random path. One hundred realizations are generated in this study.
Once the set of 100 realizations of reference precipitation at 1 km resolution has been simulated, change of support is performed to generate the coarse resolution reference precipitation fields by averaging the precipitation values simulated at 1 km grids falling within each 25 km pixel. Through this upscaling procedure, the resolution of the simulated reference precipitation becomes the same as that of the coarse resolution remote sensing product. Note that different upscaling schemes can be applied for change of support.
The third step is to compute error values at each coarse resolution pixel. Multiple realizations of the unknown reference precipitation fields are already available through geostatistical simulation and upscaling. Thus, multiple error values can be readily computed from a pixel-by-pixel comparison. Suppose there is a set of 100 simulated coarse resolution reference precipitation fields for a given number N of 25 km pixels in the study area { ( ) (V ); = 1, ⋅ ⋅ ⋅ , 100, = 1, ⋅ ⋅ ⋅ , }, and { (V ), = 1, ⋅ ⋅ ⋅ , } denote the coarse spatial resolution TRMM precipitation product. An error value at each coarse pixel is obtained per each simulation, and 100 error values are finally obtained at each 25 km pixel, as in { ( ) (V ); = 1, ⋅ ⋅ ⋅ , 100, = 1, ⋅ ⋅ ⋅ , }. The derivation of multiple error values at each pixel is the main difference between the presented geostatistical framework and the conventional approach in which only a single error value is computed for each coarse pixel. Through computation of multiple error values for each coarse pixel, the uncertainty attached to the generation of reference precipitation fields can be reflected in the computation of error values.
The final step is to summarize the multiple error values and visualize the spatial distribution of the error statistics. Since multiple errors are already computed at each coarse pixel, several summary statistics useful for interpretation can be easily derived. For example, pixel-by-pixel mean errors (ME) and mean absolute errors (MAE) can be computed to quantify the degree of bias and the magnitude of errors, respectively, as Journal of Sensors 5 Note that conventional ME and MAE values are computed over the whole study area or the whole reference dataset; thus they are location-invariant. In this study, however, the ME and MAE values are computed at each pixel by comparing one coarse resolution remote sensing product with multiple reference values. Thus, spatial error distribution, which is useful for local error analysis, can be easily visualized.
Another useful measure is the probability that the attribute values in coarse resolution products either overestimate or underestimate the reference coarse values. The probability can be derived by computing the proportion of overestimation or underestimation cases among 100 cases through a comparison between 100 reference values and the coarse resolution remote sensing product. The positive ME value means that coarse resolution product tends to overestimate the reference precipitation field. However, its meaning is slightly different from the overestimation probability. ME is computed by averaging 100 error values at each coarse pixel. Thus, the magnitude of error values greatly affects the ME value computation, while the sign of error values is only considered for the computation of over-or underestimation probability.
In addition, the relative MAE (rMAE) is also computed to evaluate the reliability of TRMM precipitation products for different locations and seasons [28,29] as The average of 100 simulated reference precipitation fields is used as the reference for rMAE. Since the precipitation amount varies per season, rMAE can be used to compare the reliability of TRMM precipitation product for different seasons. The rMAE value is also a location-dependent measure because the reference varies across the study area. If the rMAE at any pixel is less than 0.5, the TRMM product at that pixel is considered to be reliable. Conversely, if any pixel has the rMAE value equal to or greater than 0.5 of the reference precipitation amount, that pixel is considered unreliable [28,29].
GSLIB [26] was used for all geostatistical analyses including normal score transform and back-transform, variogram analysis, and SGSIM. Spatial upscaling and computation of error statistics were implemented using Fortran programming. Additionally, ArcGIS (ArcMap version 10.3, ESRI Inc., Redlands, CA, USA) was used for basic spatial data processing and visualization.

Generation of Reference Precipitation Fields at Fine and
Coarse Spatial Resolutions. One hundred realizations of reference precipitation for each month were generated at 1 km pixels using 606 AWS data and SGSIM. The October's normal score transformed data showed an anisotropic pattern along the direction of NE 70 ∘ that was caused by Typhoon Danas, which passed through south and southeastern parts of South Korea on October 8 [30]. A geometric anisotropy variogram model was fitted for October data. As other season data did not show any distinctive anisotropic patterns, isotropic variogram models were considered for April and July. Figure 4(a) presents two alternative realizations out of the 100 reference precipitation simulation results for July. Each realization honors AWS data values and reproduces the sample histogram and input variogram model. Thus, simulated precipitation realizations provide realistic representations of the unknown reference precipitation field. The upscaled realizations at a 25 km spatial resolution of the two realizations are also shown in Figure 4(b). By applying an average operator for upscaling, local details are greatly diminished, and regional variations are only reflected in the upscaled reference precipitation. In addition, the upscaled realizations do not reproduce the AWS data nor their statistical properties.
The collection of 100 realization values at each 25 km pixel constitutes the histogram and can be regarded as an approximation of the local ccdf at that pixel. Figure 5 presents the average and standard deviation of 100 realizations of reference precipitation fields for each month. The ensemble of 100 realizations, such as the average and standard deviation, depict the overall smoothed precipitation pattern and spread of reference precipitation fields, respectively. The differences among 100 realizations reflect the uncertainty prevailing at that coarse pixel and affect the error computation results.

Generation of Spatial Error Distributions.
To quantitatively assess the quality of TRMM precipitation products, error statistics including ME, MAE, and a probability that TRMM precipitation overestimates true precipitation were computed and visualized. These error statistics values represent an ensemble or a summary of 100 error values at each coarse pixel. Figures 6-8 show spatial quality distributions for three TRMM precipitation products. Note the different error values with the same color for each month. The ensemble average and standard deviation of 100 reference precipitation fields in Figure 5 were used for visual comparison and interpretation. Most pixels in April's TRMM product had positive ME values (Figure 6(a)), which means that April's TRMM product overestimated the reference precipitation fields. This observation is also corroborated by the overestimation probability shown in Figure 6(c). The overestimation probability of most pixels was very high (i.e., more than 0.75). In addition, the magnitude of errors of underestimated pixels was relatively smaller than that of overestimated pixels (Figure 6(b)). The correlation coefficient (r) between MAE (Figure 6(b)) and the average of reference precipitation fields ( Figure 5(a)) is -0.47, which implies that the magnitude of errors in April's TRMM product was moderately greater for a relatively lower amount of reference precipitation.
As shown in Figure 7, the ME and MAE values for July's TRMM product were much greater than those for April's product, implying that a heavy rainfall was not well captured in the TRMM product. The large precipitation amount in northern parts of the study area in Figure 6(  was underestimated by the July's TRMM product, yielding a greater ME value. Overestimation of reference precipitation fields was also observed in some pixels of northwestern and northeastern coastal areas. In particular, some pixels in southwestern areas underestimated the reference precipitation fields (Figure 7(c)) and exhibited large MAE values shown in Figure 7(b). However, the standard deviation values of reference precipitation fields at these pixels are relatively larger, as shown in Figure 5(b). Thus, the larger MAE values at these underestimated pixels are caused by the larger uncertainty of synthetic reference precipitation fields, not merely by the TRMM product itself, which indicates that interpretations should be made with caution. Unlike the April's TRMM product, the July's product had a positive correlation between MAE and the average of reference precipitation fields (r =0.42), implying that the magnitude of     errors in July's TRMM product was moderately greater for a relatively larger amount of reference precipitation.
October's TRMM product exhibited lower MAE values in western areas and larger MAE values in eastern areas (Figure 8(b)), which was mainly caused by anisotropic characteristics of the reference precipitation fields. Overestimation was slightly more dominant than underestimation, as depicted in Figures 8(a) and 8(c). When considering the standard deviation of simulated reference precipitation fields in Figure 5(c), large errors were moderately correlated with the uncertainty of synthetic reference precipitation fields. MAE and the average of reference precipitation fields were positively correlated, but the strength was very weak (r=0.28).
The spatial distribution of rMAE values are presented in Figure 9. Most pixels in the April's product can be considered reliable, except for some pixels located at the western coastal areas and at the southern over-land. The unreliable pixels overestimated the reference precipitation fields and showed  larger ME and MAE values. For July's product, unreliable pixels are mainly located at western and eastern coastal areas. The magnitude of errors in the July's product was greater than that in April's product, but the reliability of over-land pixels in the July's product is likely to be similar to that for April, considering the reference precipitation amount. In contrast, October's product had more unreliable over-land pixels than April's and July's products. For the northern part of the study area, the reference precipitation amount at those unreliable pixels is less than 100 mm and the uncertainty is also low ( Figure 5(c)). Thus, the TRMM precipitation at those locations is regarded consistently unreliable. When comparing the three products simultaneously, it is concluded that the reliability of October's product is the worst in land. A possible interpretation of October's product unreliability is that it was caused mainly by discrete temporal sampling errors. Errors of monthly TRMM estimates are known to be due to both discrete temporal sampling and retrieval algorithm errors   Figure 8: Spatial quality distributions for October's TRMM product: (a) ME, (b) MAE, and (c) the probability that the attribute value of the TRMM product overestimates the reference precipitation value. [29]. Particularly, temporal sampling errors range from ±8 to ±12% per month relative to the mean precipitation [29,31]. For example, the effect of a typhoon could not be properly accounted for due to the temporal discrepancy between a typhoon passing and satellite observations.
In summary, all the interpretation results for spatial distributions of error statistics can be derived from multiple reference precipitation fields, which corroborates the benefits of the simulation-based spatial quality assessment framework presented in this study.

Discussion
Error statistics estimated from the proposed approach include both errors from the coarse resolution TRMM precipitation product and the uncertainty due to the construction of synthetic reference precipitation fields. To investigate the relationships between the magnitude of errors and the uncertainty of the reference precipitations, we computed correlation coefficients between them over the study area. More specifically, each pixel in the study area has both the absolute error value per each simulation and the standard deviation value as a summary statistic of the uncertainty. Thus, the correlation coefficient value can be computed from these two values at all pixels over the study area. Since we already have 100 simulated reference precipitation fields, the correlation coefficient was computed per each realization of reference fields, yielding the distribution of 100 correlation coefficients ( Figure 10).
As shown in Figure 10, no distinct linear relationship between the magnitude of errors and uncertainty of reference precipitation fields was found. As aforementioned, the errors of April's product were negatively correlated with reference precipitation. A weak negative correlation between  absolute errors and the uncertainty of reference fields was observed in April. April's product for a relatively lower amount of reference precipitation might be interpreted as relatively reliable, but this interpretation is unclear because of its very weak correlation strength. In contrast, a weak positive correlation was observed in July and October. A certain realization for July's product showed a maximum correlation coefficient value of 0.451. As discussed before, a moderate positive linear correlation was observed between the magnitude of errors and the average of 100 realizations of reference precipitation fields for July and October. In addition, the average of the reference precipitation fields was also positively correlated with the standard deviation for both months (r=0.67 and 0.88 for July and October, respectively). These relationships were intermingled and generated positive but weak correlations for July's and October's products. If there is a strong positive correlation between the magnitude of errors and the uncertainty of the reference fields, the computed errors do not result from the coarse resolution remote sensing product, but from the uncertainty attached to the reference field generation. Therefore, caution should be exercised when interpreting the error information.
The applicability of the proposed methodology has been illustrated via an experiment with coarse resolution precipitation products, but the proposed methodology provides a general and flexible framework for spatial quality assessment of general coarse resolution remote sensing products. Other coarse resolution remote sensing products such as satellitebased soil moisture and land surface temperature can be spatially and quantitatively evaluated if validation data are available. Furthermore, the geostatistical simulation-based approach can be applied to not only regularly shaped lattice data such as remote sensing data, but also irregularly shaped blocks with different sizes. When the shape of the blocks are irregular, simulated fine resolution reference fields can be easily upscaled using points discretizing the blocks.
Despite the above advantage, however, it is difficult or even impossible to quantitatively evaluate the quality or performance of the proposed approach, because exhaustive reference field values are usually unavailable. To enhance the applicability of the spatial quality assessment framework proposed in this study, we will consider the following aspects in future work.
The uncertainty attached to the reference field prediction usually depends on both the number and the spatial distribution of available validation data. Different geostatistical simulation algorithms besides SGSIM can be applied to generate synthetic reference realizations, depending on statistical characteristics of validation data and availability of covariates. For example, if the validation data have skewed or censored distributions, nonparametric sequential indicator simulation can be applied. In addition, if high resolution covariates that are closely associated with the target attribute are available exhaustively in the area of interest, multivariate kriging algorithms such as cokriging and simple kriging with local means [25] can be used for local ccdf modeling, and the uncertainty of the reference field prediction may be reduced. For example, land-cover types, elevation, and vegetation index can be integrated with sparse ground observation data during local ccdf modeling. Thus, more research is needed to tease out the effects of characteristics of validation data on error statistics and the potential of using covariates for the reference field prediction.
In this study, upscaling was done through the equalweighted average of the simulated fine resolution reference fields, under the assumption that coarse resolution data are linear averages of fine resolution data. However, this assumption may be unrealistic for certain attributes with varying degrees of reliability in space. Nonlinear and unequal weighting functions can be applied to upscaling and their applications may yield different error estimates at coarse resolution pixels. Thus, it is necessary to evaluate the impact of different upscaling schemes on spatial quality assessment results. In addition, we will also investigate the impact of the correlation among multiple realizations on the quality assessment procedure, as well as the possibility of quantitatively separating the two error terms (i.e., coarse resolution remote sensing products and uncertainty of reference fields).

Conclusion
A geostatistical spatial quality assessment framework based on geostatistical simulation and change of support has been presented in this paper. The framework can account for (a) the uncertainty of the reference map generated from validation data and (b) the spatial resolution difference between coarse resolution remote sensing products and validation data. The main novelty of the proposed framework lies in its ability to furnish the spatial distribution of error statistics that is useful for subsequent sampling design efforts, unlike conventional evaluation procedures that fail to consider the uncertainty of reference maps and provide location-invariant global error statistics information only.
From an experiment with monthly TRMM precipitation and AWS data over South Korea, we could pinpoint the areas or pixels in which overestimation was dominant and/or large errors were observed. By using the relative magnitude of errors, we also extracted pixels with large errors with respect to synthetic reference precipitation fields, implying less reliability of those pixels.
To extend the applicability of the geostatistical spatial quality assessment framework presented in this study, further issues such as the impact of uncertainty attached to the reference field prediction on error statistics, as well as application to other coarse resolution satellite-based products, will be investigated in future work.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
Preliminary results of the research presented in this paper were presented at the Fourth International Conference on Remote Sensing and Geoinformation of Environment (RSCy 2016) held in Cyprus from 4 to 8 April, 2016.

Conflicts of Interest
The authors declare that there are no conflicts of interest.