An Overview of Theoretical and Practical Issues in Spatial Downscaling of Coarse Resolution Satellite-derived Products

Abstract

This paper presents a comprehensive overview of recent model developments and practical issues in spatial downscaling of coarse resolution satellite-derived products. First, theoretical aspects of spatial downscaling models that have been applied when auxiliary variables are available at a finer spatial resolution are outlined and discussed. Based on a thorough literature survey, the spatial downscaling models are classified into two categories, including regression-based and component decomposition-based approaches, and their characteristics and limitations are then discussed. Second, open issues that have not been fully taken into account and future research directions, including quantification of uncertainty, trend component estimation across spatial scales, and an extension to a spatiotemporal downscaling framework, are discussed. If methodological developments pertaining to these issues are done in the near future, spatial downscaling is expected to play an important role in providing rich thematic information at the target spatial resolution.

1. Introduction

In recent years, the rapid growth of computer and sensor technology has led to a rapid increase in the development and operation of various Earth observation satellites. Korea is also developing and operating its own Earth observation satellites, including the KOMPSAT (KOrea Multi-Purpose SATellite) and Geo-KOMPSAT (Geostationary KOMPSAT) series. Various satellites observe the Earth’s atmosphere, ocean, land, and polarregions using their own onboard sensors. Multiple satellite-derived outputs are also being produced frommodeling and integrated analyses, including land-cover, vegetation index,rainfall intensity, soil moisture, and land/sea surface temperatures (Cho and Suh, 2013; Kim et al., 2013; Park et al., 2015; Skofronick-Jackson et al., 2017; Xue and Su, 2017; Son and Kim, 2019). The thematic information from these satellite-derived products is often used as input for environmental modeling in various fields (Crow and Wood, 2002; Hong et al., 2007; Oh et al., 2012; Kim et al., 2019a).

The availability ofsatellite-derived products depends heavily on both spatial and temporal resolutions of input data. Satellite-derived products to be applied for global, regional, or local analyses differ according to theirspatial resolutions. The period of time-series data construction and analysis also varies, depending on temporalresolutions(Park and Kyriakidis, 2019). Most geostationary satellite productsfor global and regional monitoring are obtained with high temporal but low spatial resolutions. The low spatial resolution of satellite-derived products is not adequate for detailed local analysis. In contrast, high spatial resolution products are not very often available. Fig. 1 illustrates the difference of information content provided by satellite data having different spatial resolutions. As spatial resolution decreases, the amount of thematic information decreases accordingly. Therefore, an increase of spatial resolution is often required to fully use thematic information with high temporalresolution for environmental monitoring.

Fig. 1. Comparison of satellite data with different spatial resolutions. Numbers in parentheses indicate spatial resolution.

This kind ofscale conversion, also known as change of support within the geostatistics or spatial statistics community (Cressie, 1993), is required for combined analysis of multi-source spatial data sets with different spatial units, because spatial data are not directly obtained at the spatial resolution required for a particular analytic purpose (Zhang et al., 2014). One possible way to integrate data sets with differentspatial resolutions is to convert the fine resolution data to a coarse resolution, thereby equalizing the same spatial resolutions across data sets. However, this procedure inevitably results in loss of fine scale information. Therefore, it is also necessary to generate fine resolution spatial data for integrated analysis of multi-sensor/ source spatial data sets.

Spatial downscaling, also called disaggregation, is a scale conversion process that increases the spatial resolution (Atkinson, 2013), whereas a decrease of spatial resolution is referred to as spatial upscaling, or aggregation (Fig. 2). As illustrated in Fig. 2, spatial downscaling or upscaling can be applied to regularly spaced raster data (e.g.,satellite data) or areal data that have irregular shapes and different sizes.

Fig. 2. Illustration of spatial downscaling and upscaling.

Since the late 2000s, several spatial downscaling models have been developed and applied to obtain fine resolution thematic information from coarse resolution satellite-derived products. Although much effort has been made to develop advanced spatial downscaling models, there still remain important issues that greatly affect the performance and interpretation of spatial downscaling results; these have not yet been fully considered in previous studies.

The objectives of this paper are to provide a comprehensive overview of the state-of-the-art in spatial downscaling of coarse resolution satellitederived products and outline major issues for future research. First, we discuss the advantages and limitations of current spatial downscaling models. From a methodological viewpoint, particular attention is paid to a geostatistical downscaling framework (Kyriakidis, 2004) that can predict attribute values by reflecting the characteristics of coarse resolution input data. Second, open issues, which have been ignored in previous research, but should be considered for methodological developments and practical applications, are highlighted and discussed.

2. An Overview of Spatial Downscaling Models

In this paper, we consider spatial downscaling modelsthat predictfine resolution thematic information from coarse resolution satellite-derived products with fine resolution auxiliary variables that are associated with the coarse resolution data.

For the formulation ofspatial downscaling tasks, let z^C(v) and z^F(u) denote the attribute value of the coarse resolution satellite-derived product and its value at a fine resolution, respectively, and v and u denote a coarse resolution pixel and a point within the coarse pixel in a study area of interest, respectively. The superscriptsCand F denote coarse and fine resolutions, respectively, throughout this paper. Other inputs for spatial downscaling are N auxiliary variables that are available at a fine resolution {y^F_n(u), n=1,…,N}. Under this setting, the aim of spatial downscaling isto predict the target attribute at a fine resolution (zˆF (u)) using all available information sources. Note that the downscaling output is obtained at a resolution that is the same as that of the auxiliary variables. The coarse resolution data can be linked to the attribute values at a fine resolution via a continuous kernel or sampling function, g (Kyriakidis and Yoo, 2005): where g can take several forms such as equal or arbitrary weighting functions.

\(z^C(v) = \int_v g(u) z^F (u) du\)        (1)

In this paper, spatial downscaling models are classified into two categories, based on a literature survey: one category uses a regression approach and the other is based on component decomposition. Table 1 lists some research papers on spatial downscaling of coarse resolution satellite-derived products according to the target variables and the applied methodology. In the following section,salient theoretical aspects ofthe two approaches are given and discussed.

Table 1. Summary and classification of research papers on spatial downscaling of coarse resolution satellite-derived products

Table 1. Continued

1) Regression-based Downscaling

(1) Basic Theory

In this approach, attribute values at a fine resolution are predicted via regression modeling with fine resolution auxiliary variablesthat are closely associated with the target variable as:

\(\widehat{z} ^F (u) = f(y_n^F (u))n =1, ..., N\)         (2)

 where f(·) denotes an applied regression function.

The whole workflow of the regression-based downscaling approach is presented in Fig. 3. The coarse resolution target variable and fine resolution auxiliary variables are regarded as dependent and independent variables,respectively. Since the dependent variable for regression modeling is available only at a coarse resolution, the fine resolution auxiliary variables are first upscaled to the coarse resolution and regression modeling is then conducted at a coarse resolution. f(·) in equation (2) is constructed at a coarse resolution under the practical assumption that the statistical relationship between dependent and independent variablesremain unchanged acrossspatial resolutions. Linear/non-linearmodels or advancedmachine learning models, including geographically weighted regression, random forest,support vector regression, and artificial neural networks, have been applied as regression models(Table 1). The regression model constructed at a coarse resolution is then applied directly to the fine resolution auxiliary variables to obtain downscaling results at a fine resolution, as shown in equation (2).

Fig. 3. Workflow of a regression-based spatial downscaling model.

(2) Limitations

The regression-based approach is intuitively understandable and easy to apply for spatial downscaling. However, this approach has limitations. Even though advanced regression models are applied and auxiliary variables are very closely associated with the target variable, regression modeling seldom accountsforthe variation ofthe target variable perfectly (i.e., explanatory power cannot be 100%).Consequently, a residual component that cannot be fully explained by the auxiliary variables remains after regression modeling has been applied. However, the regressionbased approach ignores the residual component. If the auxiliary variablesfail to convey sufficient information on the target attribute, the residuals cannot be discarded, because most of the variability of the target variable is included in the residual component. Furthermore, the spatial patterns ofthe downscaling results based purely on regression modeling would be very different from those of the coarse resolution data, when the explanatory power ofthe regression model is not large.

2) Component Decomposition-based Downscaling

(1) Basic Theory

In this approach, the target variable is regarded as a realization of a random function (Z)that is decomposed into a deterministic trend component(M) and a stochastic residual component (R) across the scale (Goovaerts, 1997; Kyriakidis, 2004; Park, 2013) as:

\(Z^c(v) = M^c(v) +R^c(v)\)        (3)

\(Z^F(u) = M^F(u) + R^F(u)\)        (4)

Using the relationshipsin equations(3) and (4), fine resolution attribute values are obtained by summing the trend and residual components estimated at a fine resolution (m^F(u) and r^F(u), respectively) as:

\(\widehat{z}^F (u) = m^F (u) +r^F(u)\)        (5)

Regarding the decomposition of the target attribute into the two components, it should be noted that since there are no trend data, the true trend component cannot be available, which indicatesthat different trend models(i.e., differentregressionmodels) can be applied (Kyriakidis et al., 2004). As the residual component is computed by subtracting the trend component from the input attribute value, different trend models yield differentresidual components.Therefore, the selection of a regression model for estimation of the trend component affects both the estimation of the residual component and the final downscaling results.

As shown in Fig. 4, the fine resolution trend component is estimated using the fine resolution auxiliary variables like a regression-based spatial downscaling approach (see equation (2)). Since the regression model is constructed at a coarse resolution, the residual component is also obtained at a coarse resolution. The fine resolution residual component is then predicted from the coarse resolution residual component by applying interpolation methods. This component decomposition approach can be regarded as a hybrid approach in thatregression modeling isfirst applied for estimation of the trend component and residual correction (i.e., adding the fine resolution residual component) is then applied. Unlike the regression-based approach, there is a possibility that because of residual correction, spatial patterns of the final downscaling results would be the same as or similar to the coarse resolution data, even when the explanatory power of the applied regression model is not large.

Fig. 4. Workflow of a component decomposition-based spatial downscaling model.

(2) Prediction of Fine Resolution Residuals

Once the fine resolution trend component has been estimated, the next task isto estimate the fine resolution residual componentfrom the coarse resolution residual component (i.e., downscaling of coarse resolution residuals). Conventional point interpolation methods, including splines, point kriging, and inverse distance weighting (IDW), have been routinely applied for downscaling of residuals (Table 1). Since splines and IDW are deterministic interpolators, the residual component in equations(3) and (4)should be regarded asthe deterministic component, not the stochastic one. When applying the point interpolators, it is assumed that attribute values within each coarse pixel can be collapsed into its centroid point. As discussed in Goovaerts (2006), this assumption is reasonable if the zoom ratio (i.e., the ratio between the spatialresolution of the input data and the target fine resolution for downscaling) is not too large. As the zoom ratio increases, however, this simplification of coarse resolution pixels to points fails to account for the variations within the coarse resolution pixel.

To predict fine resolution attribute values from areal data, Kyriakidis (2004) proposed a novel geostatistical framework, called area-to-point kriging (ATPK), which can explicitly and consistently account for scale differences between input data and point predictions. Unlike conventional point kriging, ATPK does not treat coarse resolution pixels or blocks as points. Instead, it predicts attribute values at a target fine resolution by a linear combination of neighboring coarse resolution pixels or blocks. Thanksto its ability to explicitly account for spatial correlation structures across spatial scales, its application to downscale residuals has increased (see Table 1).

The practical issue ofthe application ofATPK isthat a point-support variogram model is required to compute area-to-area and area-to-point covariances in the block kriging system (Kyriakidis, 2004). Since spatial correlation structures differ acrossspatialscales, the parameters(e.g.,range and sill) ofthe areal or block variogram model are different from those of the unknown point-support variogram model.Thisimplies that the block variogram model should not be used for ATPK, particularly when heterogeneous local variations are dominant in the coarse resolution data. Variogram deconvolution is usually applied to estimate the point-support variogram model from the areal or block variogram model (Pardo-Igúzquiza and Atkinson, 2007; Goovaerts, 2008). If all area-to-area and area-to-point covariances are properly modelled, ATPK yields coherent or mass-preserving predictions, which means that upscaled ATPK predictions are the same asthe coarse resolution data values, asin equation (1). Other deterministic interpolators and point kriging cannot satisfy the coherence property perfectly.

(3) Area-To-Point Regression Kriging

Park (2013) proposed a geostatistical spatial downscaling framework that combined multiple linear regression-based trend component estimation with ATPK-based residual component estimation. This model waslater named area-to-pointregression kriging (ATPRK) by Wang et al. (2015). Kwak et al. (2018) also developed the prototype R-based tool, named R4ATPRK, for the implementation of ATPRK.

ATPRKpredictions can perfectly satisfy the coherence property because linearregression andATPKare adopted for estimation of trend and residual components, respectively (Park, 2013). Coherent predictions of ATPRK indicate thatATPRK can generate downscaling results that are consistent with the coarse resolution data.Recently, non-linearregression models, including geographically weighted regression and random forest, have been widely employed within an ATPRK framework (Table 1). However, if a non-linear model is applied, coherence cannot be perfectly satisfied. When a non-linear model is applied but coherence should be satisfied, the following normalization of the fine resolution trend component can be combined with downscaling of coarse resolution residuals byATPRK. The trend component estimated at a fine resolution is upscaled to the coarse resolution and the coarse resolution residual component is then computed by subtracting the new upscaled trend componentfromthe coarse resolution data (Hutengs and Vohland, 2016). If ATPK is employed to estimate the fine resolution residual component, the final downscaling result can satisfy the coherence property.

3. Open Issues and Challenges

Advanced modelssuch asATPRK can be promising for spatial downscaling of coarse resolution satellitederived product. However, practical issuesto be further developed or considered stillremain.These are pointed out and discussed, with examples, in the following section.

1) Comparison of Different Regression and Interpolation Models

In both regression-based and component decomposition-based approaches, different regression models can be applied to estimate the fine resolution trend component (Table 1). The performance of a regression model is generally dependent on the quality of the coarse resolution input data and associated fine resolution auxiliary variables, as well as the physical and geographical conditions of the study area. Thus, there is no single best regression model for estimation ofthe fine resolution trend component and a quantitative comparison of different regression models is required for spatial downscaling.

The performance offine resolution residual estimates also depends on an applied interpolator. To compare the characteristics of different interpolators for spatial downscaling of residuals, an experiment on spatial downscaling of normalized difference vegetation index (NDVI) data was undertaken. First, MODIS NDVI data at 250 m were upscaled to four different spatial resolutions(1 km, 2 km, 5 km, and 10 km). The 10 km NDVI data (Fig. 5(a)) were then downscaled to the 5 km, 2 km, and 1 km resolutions(i.e., zoom ratios are 2, 5, and 10, respectively). The NDVI data upscaled to 1 km, 2 km, and 5 km were regarded as unknown true data at each resolution and used to compute error statistics(e.g., root mean square error (RMSE)) for the downscaling results. Fourinterpolators, including IDW, spline with tension, point ordinary kriging, andATPK, were compared for this experiment because they have been commonly applied to residual correction.To apply the first three interpolators, 10 km pixels were treated as points and then used for spatial interpolation. The point-support variogram model required for the application of ATPK was estimated using variogram deconvolution and the points at target resolutions within each 10 km pixel were used to compute area-toarea and area-to-point covariance values.

Fig. 5. Spatial downscaling results at 1 km resolution using 10 km NDVI data and different interpolators: (a) NDVI data at 10 km resolution, (b) IDW predictions, (c) spline predictions, (d) point kriging predictions, and (e) ATPK predictions. The black rectangles denote the 10 km pixels.​​​​​​​
 

The downscaling results at 1 km resolution are presented in Figs. 5(b)-(e) and the summary statistics of the downscaling results are listed in Table 2. Since only distances between the centroids of 10 kmresolution pixels andthe1kmprediction-gridpointswere considered in IDW, an unrealistic bull’s-eye effect around the centroid is shown in the IDW predictions (Fig. 5(b)). The point kriging-based predictions show strong smoothing effects with less spatial variability (Fig. 5(d)). Similar spatial patterns were observed in spline and ATPK predictions.

 Table 2. Summary statistics of 10 km NDVI data, true NDVI data, and downscaling results from different interpolators with respect to different zoom ratio values. Correlation refers to the linear correlation coefficient between 10 km NDVI values and upscaled predictions

The summary statistics in Table 2 clearly reveal the different characteristics of the four interpolators. All predictionsshowed mean valuesthat were the same as, or similar to, those of coarse resolution data and true fine resolution data, indicating unbiased predictions. However, the standard deviations that quantify the variability of predicted values were quite different. As expected from Fig. 5(d), the standard deviation of point kriging was the smallest and IDW was the next smallest.ATPK exhibited the largeststandard deviation, yielding more variability than spline predictions, as shown in Figs. 5(c) and (e). The different standard deviationsresulted in different prediction ranges of the four interpolators. IDW predictions exhibited the minimum and maximum values that were similar to those of the 10 km input data, whereas ATPK and splines returned attribute values that were either larger than the maximum orsmallerthan the minimum values of the input data. Because the variability of attributes at a finer resolution is usually greater than at a coarse resolution, this characteristic of bothATPK and splines is more appropriate for spatial downscaling. It is noteworthy that the range of point kriging predictions was the smallest, although kriging is a non-convex interpolator. The smallest range of point kriging may be explained by the variogram model of NDVI values at the centroids of which the total sill value was the same as the variance of 10 km NDVI data. When looking into the correlation between the upscaled predictions and the 10 km NDVI data,ATPK perfectly satisfied the coherence property and splinesshowed the second largest correlation. When comparing the prediction accuracy using RMSE, the most accurate interpolator was ATPK and the second most accurate was a spline interpolator. Point kriging was the worst interpolator, indicating that the application of point kriging after simplification of areal data or pixels to their centroids was not appropriate. This experiment confirms that ATPK can be a useful interpolator for downscaling coarse resolution residual components, but it is still necessary to compare the effects ofATPK and a spline interpolator on the final downscaling results.

2) Scale Invariant Assumption in the Trend Component Estimation

An important issue that has not been fully considered in traditional regression-based spatial downscaling approaches is the assumption of preservation of statistical relationships across spatial scales, (i.e., the scale invariant assumption). In other words, the statisticalrelationships quantified at a coarse resolution are assumed to remain unchanged at a finer resolution, thus, these relationships have been directly applied to fine resolution auxiliary variables. This practical but strong assumption is valid only when the linear regression model is applied, but non-linear models cannot satisfy this assumption. If specific conditions are not satisfied, regression coefficients may have serious bias (Woojoo Lee, personal communication). The resulting biasis often referred to as ecological bias (Waller and Gotway, 2004). Possible solutionsto solve this problem include a continuous stochastic Gaussian process model (Tanaka et al., 2018) or a Bayesian hierarchical model (Keil et al., 2013), and the performance of the models should be investigated.

3) Quantification of Uncertainty

Downscaled thematic maps are usually used as inputs for environmental modeling. If significant uncertainties or errors are included in the downscaling results, the model outputs are greatly affected by those uncertainties or errors. Therefore, it is critical to quantify uncertainty attached to spatial downscaling.

Similar to conventional point kriging,ATPK returns error variance as well as prediction values. Kriging error variance reflects the sample data configuration, but isindependent ofsample values(Goovaerts, 1997). Consequently, ATPK error variance cannot be used as ameasure of uncertainty attached to spatial downscaling. Any kriging algorithm provides prediction results that are optimal in a least-square sense. As discussed in Boucher and Kyriakidis (2006), spatial downscaling should be considered as an under-determined inversion problem, because there are many possible attribute values at a fine resolution that satisfies the coherence property.

In this context, stochastic simulation, which can generate multiple alternative realizations of unknown true values at a fine resolution (Goovaerts, 1997; Kyriakidis and Yoo, 2005), may be a possible solution to quantification of uncertainty. When applying stochastic simulation to spatial downscaling, the target for simulation is the residual component at a coarse resolution and multiple realizations of the residual component are obtained at a fine resolution.By adding them to the deterministic trend component at a fine resolution, multiple realizations of the unknown target attribute value, which are conditional to all available information and reasonably reproduce the spatial correlation structure, are obtained at a fine resolution.

Fig. 6 illustrates a simulation-based downscaling approach (Park et al., 2018). The AMSR-2 rainfall intensity at 25 km resolution was downscaled using COMS brightness temperatures from infrared and water vapor channels at 5 km resolution. A total of 50 realizations were experimentally generated using simulation of kriging errors(Chilès and Delfiner, 1999) and the 10th and 30th realizations out of 50 realizations are shown in Figs. 6(b) and (c).The differences between multiple realizations (e.g., Figs.6(b) and (c)) can be used to quantify the uncertainty attached to spatial downscaling. If the variability of downscaled rainfall intensity values among multiple realizations is large at a certain pixel, that pixel may be regarded as one with large uncertainty. Since multiple realization values are already obtained at each pixel, a probability of exceeding a certain thresholding value can be easily computed for a probabilistic interpretation. Figs. 6(d) and (e) represent a probability of exceeding a rainfall intensity of 10 mm/h and 20 mm/h, respectively.

Fig. 6. Example of simulation-based downscaling: (a) AMSR-2 rainfall intensity at 25 km, (b) 10th simulation at 5 km, (c) 30th simulation at 5 km, (d) probability of exceeding 10 mm/h, and (e) probability of exceeding 20 mm/h. The black rectangles and polylines denote the 25 km AMSR-2 pixels and administrative boundaries, respectively​​​​​​​.

The simulation-based approach, discussed above, quantifies the uncertainty only from the residual component estimation, because the trend component is assumed to be deterministic.When the trend component is regarded as stochastic, the uncertainty attached to spatial downscalingmay also come fromthe uncertainty or errors from the trend component estimation. It is, therefore, necessary to develop a stochastic simulation framework that can simultaneously accountforthe two uncertainty sources.

4) Errors of Coarse Resolution Data

Even though satellite-derived products have been produced with high quality, any products derived from satellite remote sensing data inevitably include systematic, random, and modeling errors. In turn, the errors ofthe satellite-derived products greatly affect the downscaling results, because the performance ofspatial downscaling is subject to the accuracy or reliability of coarse resolution input data. If the coarse resolution data include significant errors, the over- and/or mis-fitted trend component seriously affects the performance of downscaling. As discussed in Park et al. (2016) and Kim and Park (2017a), the regression models with higher explanatory power do not always lead to an improvement of prediction performance, and fitting with higher-order polynomials and machine learning algorithms may degrade the accuracy of downscaling. Therefore, the explanatory power of regression models alone cannot be used for the selection of an optimal trend estimationmodel in spatial downscaling (Kim and Park, 2017a) and comparisons of different regression models should be made with caution. Furthermore, when the explanatory power of the regression models is not large, the residual component is proportionally large and most of the variability of the input data isincluded in the residuals. In this case, the errors of input data also affect the prediction of the residuals. Therefore, it is of great importance to quantify and considerthe errors contained in the satellite-derived products when comparing differentregression models and interpolators(Kim and Park, 2017b).

In relation to error correction, if information concerning error or uncertainty is available at the coarse resolution pixel (e.g., random errors in TRMM and GPM IMERG products), this information can be first used for error correction of the satellite-derived products and the error-filtered data are then used as inputs for spatial downscaling. If information concerning error of satellite-derived products is not available, it can be estimated using geostatistical modelingwith ground observations(Park andKyriakidis, 2019). Christensen (2011) proposed filtered kriging that accounts for error variance information to predict error-filtered values in climate model outputs. However, filtered kriging treats satellite-products as their centroid points; thus, it cannot account for variations within coarse resolution pixels. Recently, Kim et al. (2019b) proposed area-to-area filtered kriging (ATAFK) to filter out an error or uncertainty component in coarse scale satellite-based products without the simplification of coarse resolution pixelsto points. This can be regarded as a hybrid kriging algorithm that combines the basic concept of filtered kriging with area-to-area block kriging. Fig. 7 illustrates the error-corrected GPM product by ATAFK, in which the attribute values of any pixel with large error values have been filtered out. Although ATAFK can be a promising tool as a preprocessor for spatial downscaling, it should be noted that the error information at the coarse resolution pixel may not be related directly to point-support ground observations, thusit is necessary to use theATAFK result as an input forspatial downscaling and then investigate itsimpact on the performance of downscaling.

Fig. 7. Illustration of ATAFK to filter out the errors in the coarse resolution satellite-derived products: (a) TRMM rainfall intensity data, (b) corresponding error information, and (c) filtered result by ATAFK. The black rectangles and polylines denote the 25 km TRMM pixels and administrative boundaries, respectively​​​​​​​.

5) Evaluation of Downscaling Results

The quality of spatial downscaling results and the performance of different downscaling models have been evaluated via direct comparisons against ground observations or estimated true maps. A small number of the ground observations are usually available for comparison, thus, not all pixels in the downscaling results can be assessed. Forspatial downscaling ofland surface temperature (LST) products such as MODIS LST products, LST values at a fine spatial resolution are firstretrieved using Landsat thermal band data.The retrieved LSTvalues are then regarded astrue data and compared with the downscaling results (Hutengs and Vohland, 2016;Yang et al., 2017). However, errors are inevitably included during the retrieval of LST values at a fine spatial resolution. As a result, quantitative evaluation of the downscaling resultsis affected by the errors contained within the estimated LST values.

Anotherlimitation oftraditional evaluation approaches is that they fail to account for the spatial resolution or scale difference between the point-support observations and areal data or pixels (Park and Kyriakidis, 2019). The output ofspatial downscaling is obtained at a finer spatialresolution, but it isstillregarded asthe areal data. For example, when satellite-derived products at 10 km resolution are downscaled to 1 km resolution, the downscaling results are obtained at 1 km pixels or blocks, not at 1 km-interval points, implying a direct comparison of these areal data against point-support ground observation data cannot lead to a proper quantitative evaluation. For example,soil moisture and land surface temperature have much larger variability depending on observation times and locations. Thus, location-specific observation values may be significantly differentfrom the aggregated pixel values and this fact should be considered when interpreting quantitative evaluation results.

In this context, a supplementary quantitativemeasure should be further considered to evaluate downscaling results in addition to traditional accuracy measures. When reminding that spatial downscaling aims at predicting thematic information at a finer resolution, an increase of information content through spatial downscaling should be considered to evaluate downscaling results. For example, information entropy that indicates the richness of information (Zou et al., 2015) can be a quantitative measure. Large information entropy values indicate that more detailed spatial variations are depicted in the downscaling results.Thus, when comparing differentspatial downscaling models, the information entropy can be used as a supplementary measure. Park et al.(2016) demonstrated the usefulness of information entropy as an evaluation measure in spatial downscaling of AMSR-2 precipitation with COMSobservations.Therefore,the depiction of detailed spatial patternsshould be evaluated along with the error statistics based on comparison with point-level ground observations.

6) Consideration of Data Characteristics

In addition to the data that are unbounded in the real-numberspace,some data have their own constraints to their values. For example,soil moisture valuesrange between 0 (0 %) and 1 (100 %), and precipitation intensity takes non-negative values.Land-coverfractions have both non-negative values and a constant sum of 100 % or 1. If these data are used as inputs for spatial downscaling, any spatial downscaling model should take into account these variable-specific characteristics. However, traditionalspatial downscalingmodels do not fully account for the constraints of the input data. For example, negative values that are physically meaningless are usually generated around and within the coarse resolution pixels having relatively small values.

For the proper consideration of the characteristics of input data, a procedure for optimization of the initial downscaling results should be included into the downscaling process to satisfy the constraints, as quadratic programming conjunction with ATPK in Yoo and Kyriakidis(2006).In addition to the inclusion of a post-processing procedure, a series of data transformationsspecific to the input data type can also be applied to spatial downscaling. For example, a log-ratio approach that has been frequently applied to compositional data (Pawlowsky-Glahn and Buccianti, 2011) can be employed forspatial downscaling oflandcover fractions (Park et al., 2018).

7) Spatio-temporal Downscaling

As mentioned in the Introduction, products derived from geostationary satellites are temporally dense but spatially coarse, whereas productsfrom low Earth orbit satellites are spatially fine but temporally sparse. For the synergistic use ofthese complementary productsin terms of both spatial and temporal resolutions, spatial downscaling can be extended to a spatio-temporal downscaling framework that aims at generating dense time-series fine resolution datasets (Fig. 8). The key idea ofspatio-temporal downscaling isto first quantify temporal correlation information from dense timeseries coarse resolution data sets and then to use this information to predict fine resolution thematic information at a time when fine resolution data are not acquired.

Fig. 8. Illustration of a spatio-temporal downscaling framework.​​​​​​​

Several fusion models, including a spatial and temporal adaptive reflectance fusionmodel(STARFM: Gao et al., 2006), an enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM: Zhu et al., 2010), and flexible spatio-temporal data fusion (FSDAF: Zhu et al., 2016), have been proposed and applied to spatio-temporal downscaling (e.g., blending reflectance valuesfrom MODIS andLandsat). Despite their great potential for spatio-temporal downscaling, there remain some aspects requiring further improvement. Traditional models only quantified the temporal correlation information when the dense time-series coarse resolution data and the sparse time-series fine resolution data were acquired simultaneously. That is, they could not take into account temporal correlation information from the dense time-series coarse resolution datasets. As a result, the prediction performance may degrade as the difference between the acquisition time and the prediction time increases. Therefore, advanced spatiotemporal downscalingmodelsshould be developed that can fully take into account the temporal correlation information.

4. Conclusions

With the increasing availability of multi-sensor satellite data, satellite-derived products have been widely applied across various fields related to Earth observation. However, it is not always easy to acquire satellite-derived products at an appropriate spatial resolution. Thus, there is a strong need for scale conversion to extract thematic information at a target resolution, as well as to integrate multiple spatial datasets with different supports. Spatial downscaling, which is a process to predict attributes at a finer resolution from coarse resolution data, should be regarded not as a conventional spatial estimation or prediction task, but as a challenging research topic, since many aspects need to be taken into account.

In this context, this paper presented a comprehensive framing of theory and issues in spatial downscaling of coarse resolution satellite-derived products. One particular focus was on a geostatistical framework that properly treats coarse resolution data and generates coherent predictions.Afterthe state-of-the-art ofspatial downscaling models had been outlined, several practical issues and challenges for future research, including quantification of uncertainty attached to spatial downscaling, rigorous estimation of trend and residual components, impacts of erroneous input data on downscaling performance, model developments considering data constraints, and spatio-temporal downscaling, were discussed and demonstrated, with examples provided based on synthetic and real data sets. These future research issues were individually detailed and discussed, but each issue can be unified with others. For example, a stochastic downscaling framework presented for quantification of uncertainty can be extended to the spatio-temporal downscaling framework and combined to correct the intrinsic errors contained in coarse resolution data.As anotherresearch topic not discussed in this paper, determining the optimalfinestspatialresolution forspatial downscaling should also be investigated.

One misconception about spatial downscaling is to expect an improvement in quality or accuracy after spatial downscaling.Itshould be noted that the ultimate goal ofspatial downscaling isto obtain detailed spatial patterns at a finerresolution, not to generate results with better accuracy than the input data, because the quality or accuracy of spatial downscaling results is subject to the quality of the input data. If this basic concept of spatial downscaling is properly adapted to the application at hand and unresolved issues outlined in this paper are further investigated,spatial downscaling will contribute to the expansion of the potential of satellite-derived products across various fields related to Earth observation.​​​​​​​​​​​​​​