Next Article in Journal
Evaluation of ISS-RapidScat Wind Vectors Using Buoys and ASCAT Data
Next Article in Special Issue
Reconstruction of Single Tree with Leaves Based on Terrestrial LiDAR Point Cloud Data
Previous Article in Journal
Evaluation of Heavy Precipitation Simulated by the WRF Model Using 4D-Var Data Assimilation with TRMM 3B42 and GPM IMERG over the Huaihe River Basin, China
Previous Article in Special Issue
Spatio-Temporal Analysis and Uncertainty of Fractional Vegetation Cover Change over Northern China during 2001–2012 Based on Multiple Vegetation Data Sets
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Spatial Downscaling of Gross Primary Productivity Using Topographic and Vegetation Heterogeneity Information: A Case Study in the Gongga Mountain Region of China

1
Research Center for Digital Mountain and Remote Sensing Application, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China
2
University of Chinese Academy of Sciences, Beijing 100049, China
*
Author to whom correspondence should be addressed.
Submission received: 14 March 2018 / Revised: 16 April 2018 / Accepted: 20 April 2018 / Published: 22 April 2018

Abstract

:
Due to the spatial heterogeneity of land surfaces, downscaling is an important issue in the development of carbon cycle models when evaluating the role of ecosystems in the global carbon cycle. In this study, a downscaling algorithm was developed to model gross primary productivity (GPP) at 500 m in a time series over rugged terrain, which considered the effects of spatial heterogeneity on carbon flux simulations. This work was carried out for a mountainous area with an altitude ranging from 2606 to 4744 m over the Gongga Mountain (Sichuan Province, China). In addition, the Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product at 1 km served as the primary dataset for the downscaling algorithm, and the 500 m MODIS GPP product was used as the reference dataset to evaluate the downscaled GPP results. Moreover, in order to illustrate the advantages and benefits of the proposed downscaling method, the downscaled results in this work, along with ordinary kriging downscaled results, spline downscaled results and inverse distance weighted (IDW) downscaled results, were compared to the MODIS GPP at 500 m. The results showed that (1) the GPP difference between the 500 m MODIS GPP and the proposed downscaled GPP results was primarily in the range of [−1, 1], showing that both vegetation heterogeneity factors (i.e., LAI) and topographic factors (i.e., altitude, slope and aspect) were useful for GPP downscaling; (2) the proposed downscaled results (R2 = 0.89, RMSE = 1.03) had a stronger consistency with the 500 m MODIS GPP than those of the ordinary kriging downscaled results (R2 = 0.43, RMSE = 1.36), the spline downscaled results (R2 = 0.40, RMSE = 1.50) and the IDW downscaled results (R2 = 0.42, RMSE = 1.10) for all Julian days; and (3) the inconsistency between MODIS GPP at 500 m and 1 km increased with the increase in altitude and slope. The proposed downscaling algorithm could provide a reference when considering the effects of spatial heterogeneity on carbon flux simulations and retrieving other fine resolution ecological-physiology parameters (e.g., net primary productivity and evaporation) over topographically complex terrains.

Graphical Abstract

1. Introduction

The gross primary productivity (GPP), which corresponds to carbon fixation by vegetation at the ecosystem level, is an important indicator to assess the photosynthetic capacity of vegetation and the function of the ecosystem. The accurate estimation of regional, continental and global GPP plays an important role when monitoring vegetation growth conditions [1], the terrestrial carbon budget [2] and the interactions in the soil-vegetation-atmosphere continuum [3].
Carbon cycle models have been developed over many years to predict carbon-climate feedbacks at multiple scale levels [4,5,6], which can be used to estimate spatially and temporally continuous GPP. However, most carbon cycle models are generated from the site scale, without considering the spatial heterogeneity within each modeling grid [7]. Moreover, due to the paucity of the data and computing complexities, the carbon cycle model is usually executed at coarse resolutions, which are based on the simplification of landscape complexities. Such a simplification results in large uncertainties due to spatial heterogeneity, especially in topographically complex terrains. Some researchers have shown that the oversimplification of landscape complexities may inevitably cause model simulations to be considerably biased [8,9,10]. Landscape complexities typically present high spatial heterogeneity, and spatial heterogeneity is scaled by multiple factors, including endogenous and exogenous factors [11,12]. Endogenous heterogeneity refers to the spatial variability of vegetation types and density (i.e., land cover and leaf area index), while exogenous heterogeneity is associated with surface topography (i.e., altitude, slope and aspect). Therefore, it is crucial to consider the effects of spatial heterogeneity on carbon flux simulation by models.
The carbon flux simulation through carbon cycle models usually requires massive forcing data. Atmospheric data are one of the most important inputs for carbon cycle model simulation. Theoretically, atmospheric forcing data, such as air temperature, air humidity and solar radiation, describe the external environment for the carbon cycle process and exhibit strong spatial variation characteristics over complex terrain [13,14,15]. Dense field observations are required to determine the spatial patterns of atmospheric forcing data, especially for rugged regions. However, due to insufficient meteorological stations located in mountainous area, the atmospheric forcing data could not be provided accurately at fine resolution [12]. To overcome the limitation of field meteorological observations for running the carbon cycle model, downscaling is a potential method to obtain high resolution carbon flux simulation in topographically complex terrains. Using the relationships between carbon flux estimation and the surrounding landscape, the disaggregation of low spatial resolution carbon flux estimation and the information from ancillary data could be introduced in the downscaling scheme.
Spatial scaling is a process of using information available at one scale to derive processes at another scale, including upscaling and downscaling. Downscaling refers to an increase in spatial scale following disaggregation of the coarse-resolution dataset and information from ancillary data at a finer resolution, which restores the variations at a finer scale by assuming the values of the coarser are the mean values at the finer scale. In the last decade, the downscaling of the carbon cycle process has been one of the most challenging issues in environmental science, especially over mountainous areas [16,17,18]. In the literature, some studies have been carried out to examine the effects of spatial heterogeneity on the upscaling process of carbon fluxes [8,19], while little attention has been shown towards the issue of downscaling carbon fluxes over complex terrains. Therefore, it is necessary to downscale carbon fluxes using the subpixel information of topography and vegetation heterogeneity, which is an important issue when evaluating the role of ecosystems in the global carbon cycle.
Vegetation heterogeneity and surface topography are important factors that introduce biases into carbon modeling and water fluxes. Vegetation heterogeneity directly influences the photosynthetic capacity per unit of surface area. Topography affects the redistribution of water, radiation and heat, which have significant influences on the process of carbon assimilation. In this paper, the main purposes are (1) to develop a GPP downscaling algorithm using the subpixel information of topography and vegetation heterogeneity and (2) to assess the GPP inconsistency at 1 km and 500 m, to analyze the relationship between spatial heterogeneity and the GPP difference at two scales.

2. Materials and Methods

2.1. Study Area

As shown in Figure 1, the experimental region is northwest of Gongga Mountain on the quaternary sections of the Tibetan plateau. The summit of Gongga Mountain is 7556 m above sea level, which is the highest mountain in the Sichuan Province of China. Mt. Gongga has a mountainous climate that ranges from a cool plateau climate to a subtropical lowland climate [20]. The annual precipitation is approximately 3000 mm at the summit and declines to 1300 mm at lower altitudes [21]. The average annual air temperature decrease with altitude, from approximately 11 °C at lower altitudes to −1 °C at higher altitudes [21]. This experimental region, with an area of 2584 km2, has an altitude ranging from 2606 to 4744 m above sea level, and the dominant types of land cover are forestlands (i.e., evergreen forests and mixed forests), grasslands and shrublands.

2.2. Data and Processing

2.2.1. GPP Images at 500 m and 1 km

The Moderate Resolution Imaging Spectroradiometer (MODIS) has produced a regular estimation of GPP since 2000, including the MOD17A2 (8-day, 1 km) and MOD17A2H (8-day, 500 m) products. The algorithm for the MODIS GPP (gC m−2 day−1) is based on the concept of radiation conversion efficiency and can be described as follows:
GPP = ε max × T M I N _ s c a l a r × V P D _ s c a l a r × S W R a d × 0.45 × F P A R
where εmax represents the maximum light use efficiency (gC MJ−1) related to the plant functional type (PFT), which can be obtained from the Biome Properties Look-Up Table (BPLUT); TMIN_scalar and VPD_scalar are the scalars for minimum air temperature and vapor pressure deficits (VPD), respectively, ranging from 0 to 1; SWRad represents the incoming short wave radiation (MJ m−2 day−1), and FPAR represents the fraction of absorbed photosynthetic active radiation, which can be obtained from MODIS product. The MODIS GPP were scaled by multiple environmental factors, and each environmental factor has a comparable level of influence on the final simulated GPP. More details regarding the MODIS GPP algorithm can be found in the literature [22,23]. In addition, some research has been done to evaluate MODIS GPP products across multiple biomes, suggesting that the MODIS GPP works effectively for the majority of PFTs [23,24,25].
The MODIS GPP products MOD17A2 (version 5.0) and MOD17A2H (version 6.0) used in our study were obtained from Julian day 169, 2010 to Julian day 209, 2010, which can be freely acquired at National Aeronautics and Space Administration Distributed Active Archive Center (http://www.modis.ornl.gov/modis/index.cfm). Using spatial heterogeneity information within each 1 km modeling grid, our study was designed to estimate downscaled GPP from a 1 km resolution to a 500 m resolution. Because the 1 km GPP product is executed at a coarse resolution without considering the spatial heterogeneity within each 1 km modeling grid, the GPP product at 1 km served as the primary dataset for the downscaling algorithm in this work. The 500 m GPP product was used as a reference dataset to evaluate the downscaled GPP results at 500 m.

2.2.2. LAI Data

The MODIS LAI products MCD15A2 (version 5.0) at 1 km and MCD15A2H (version 6.0) at 500 m, were selected to describe the vegetation density in this study. The 8-day MODIS LAI product is available at http://wist.echo.nasa.gov in the sinusoidal projection. According to the quality control layer (FparLai-QC) involved in this product, the MODIS LAI consists of a main look-up-table (LUT) algorithm (QC < 64) and a back-up algorithm (64 ≤ QC < 128). The LUT algorithm exploits surface reflectance information in the red and near-infrared bands based on a 3-D radiative transfer equation. When the LUT algorithm fails, the back-up algorithm adopts NDVI-LAI empirical relationships to estimate the LAI [26]. The LUT algorithm usually has a higher reliability than that of the back-up algorithm [26]; the clumping effects at the canopy scale are considered in the LUT algorithm [27].

2.2.3. DEM Data

Some researches have illustrated that the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM) data have a good elevation accuracy over mountainous area [28,29,30], ASTER GDEM were adopted to describe the complex terrain in this work. GDEM data available at a 30 m were averaged to obtain the values at 500 m and 1 km resolutions [7]. The surface slope and aspect were derived from the GDEM data, according to a fast algorithm based on a 3 × 3 pixels window [31].

2.3. Algorithm for Downscaling

As shown in Figure 2, the proposed downscaling algorithm contains two phases: the regression process and area-to-point kriging (ATPK) process. It firstly conducts regression analysis between GPP and heterogeneity factors (i.e., altitude, slope, aspect and LAI), and then ATPK is performed on the downscaling residuals from the regression process.

2.3.1. Problem Formulation

Let GPPR(xj) represent the random variable of pixel R centered at xj (j = 1, 2, …, N, where N is the number of pixels) in 1 km resolution image of GPP, whereas GPPr(xi) represents the random variable of pixel r centered at xi (i = 1, 2, …, M, where M is the number of pixels) in the 500 m resolution image of GPP. The notations R and r represent the 1 km and 500 m pixels, respectively. The purpose of the downscaling method based on topographic factors and the LAI is to predict GPPr(x) at a 500 m resolution. GPPr(x) can be described as
G P P r ( x ) = G P P r 1 ( x ) + G P P r 2 ( x )
where GPPr1(x) and GPPr2(x) represent the predictions of the regression and ATPK processes, respectively. More descriptions about their calculation can be found in the following sections.

2.3.2. Regression between GPP and Topographic Factors

The regression takes advantage of valuable spatial heterogeneity information from the ancillary data at 500 m spatial resolution. The relationship between topographic factors and MODIS GPP at 1 km is modeled by multiple linear regression:
G P P R _ o b s ( x ) = b 0 + b 1 H R + b 2 C O S ( S R ) + b 3 C O S ( A R ) + b 4 L A I R
where GPPR_obs(x), HR, SR, AR and LAIR represent the MODIS GPP(gC d−1), altitude (m), slope (deg), aspect (deg) and LAI (m2 m−2) at a 1 km resolution, respectively. Since GPPR_obs(x), HR, SR, AR and LAIR in Equation (3) are known, b0, b1, b2, b3 and b4 can be estimated. Currently, a wide range of fitting methods has been applied in the theoretical framework of regression, such as generalized least squares (GLS) and ordinary least squares (OLS). In this paper, the OLS algorithm was adopted to estimate the regression coefficients.
The relationship in Equation (3) is assumed to be universal at different resolutions [32,33], and thus the relationship modeled at a 1 km resolution can be implemented at a 500 m resolution. Based on this assumption, the regression prediction GPPr1(x) is calculated as
G P P r 1 ( x ) = b 0 + b 1 H r + b 2 C O S ( S r ) + b 3 C O S ( A r ) + b 4 L A I r
where Hr, Sr, Ar and LAIr represent the altitude (m), slope (deg), aspect (deg) and LAI (m2 m−2) at a 500 m resolution, respectively. Variables b0, b1, b2, b3 and b4 are regression coefficients.

2.3.3. ATPK for Downscaling Residuals

If the result of the first phase is perfect, there should be no bias between it and the original GPP at a 1 km resolution. However, it is unrealistic to obtain such an ideal result, and there are inevitable residuals from the first process. The residuals at a 1 km resolution are calculated as:
G P P R _ r e s i d u a l ( x ) = G P P R _ o b s ( x ) ( b 0 + b 1 H R + b 2 C O S ( S R ) + b 3 C O S ( A R ) + b 4 L A I R )
In this work, the residuals at a 1 km resolution were considered for downscaling, and the ATPK process served as the second phase when downscaling the residuals, GPPR_residual(x), at 1 km to those at 500 m, GPPr_residual(x). In the ATPK approach, GPPr2(x) represents a linear combination of residuals from N 1 km pixel, which can be described as:
G P P r 2 ( x ) = i = 1 N λ i G P P R _ r e s i d u a l ( x i ) , s . t . i N λ i = 1
where λi represents the weight for the residual of the ith 1 km pixel, which is centered at xi. As shown in Equation (6), the ATPK phase takes the spatial correlation among pixels into account, which is not considered in the above equations.
The main issue of the ATPK phase is obtaining the weights, λ, which are estimated by minimizing prediction error spatial variations at a 1 km resolution. The kriging system can be expressed as:
[ C RR ( x 1 , x 1 ) C RR ( x 1 , x N ) 1 C RR ( x N , x 1 ) C RR ( x N , x N ) 1 1 1 0 ] [ λ 1 λ N μ ] = [ C r R ( x , x 1 ) C r R ( x , x N ) 1 ]
where CRR(xi, xj) is the residual covariance between a 1 km pixel centered at xi and a 1 km pixel centered at xj, and CrR(x, xj) is the residual covariance between a 500 m pixel centered at x and a 1 km pixel centered at xj. Variable μ represents the range multiplier. To estimate weight λ according to Equation(7), these two types of residual covariance need to be calculated in advance.
Assuming each pixel as a point, s is intended to represent the Euclidean distance between two pixels. The 1 km to 1 km residual covariance, CRR(s), and the 500 m to 1 km residual covariance, CrR(s), can be estimated as:
C R R ( s ) = C r r ( s ) h R ( s ) h R ( s )
C r R ( s ) = C r r ( s ) h R ( s )
where Crr(s) represents the 500 m to 500 m residual covariance, and hR(∗) represents the point spread function (PSF). Variable −s indicates that the distance from a 500 m pixel (A) to another 500 m pixel (B) within a 1 km pixel is opposite to the distance from a 500 m pixel (B) to a 500 m pixel (A) (i.e., s).
Supposing the value of the coarse pixel is the mean of the 500 m pixels within it, and F is the ratio between the resolutions of 1 km pixel and 500 m pixel, then, PSF can be described as:
h R ( x ) = { 1 F 2 x R ( x ) 0 o t h e r w i s e
Given Equations (8)–(10), the calculation of CRR(xi, xj) and CrR(x, xj) can be simplified as:
C R R ( x i , x j ) = 1 F 4 m = 1 F 2 m = 1 F 2 C r r ( s m m )
C r R ( x , x j ) = 1 F 2 m = 1 F 2 C r r ( s m )
where Smm’ is the distance between each pair of 500 m pixels from within two 1 km pixels centered at xi and at xj, respectively, and sm is the distance between 500 m pixels centered at x and any 500 m pixels within the 1 km pixel centered at xj. Therefore, the critical issue when estimating the weight is the estimation of the 500 m to 500 m residual covariance, Crr(s), which is essentially semivariogram modeling [32,34].

2.4. Result Validation and Method Evaluation

In this study, the performance of the proposed downscaling method was evaluated by MODIS GPP at a 500 m resolution. The coefficient of determination (R2) and the root mean square error (RMSE) between the downscaled results and the 500 m MODIS GPP were adopted as indicators of validation. In addition to the proposed downscaling approach, other three downscaling methods, including ordinary kriging (OK) [35], spline [36] and inverse distance weighting (IDW) [37], were used to provide a systematic comparison and illustrate the advantages of the proposed downscaling method.

3. Results and Analyses

3.1. Topographic and Vegetation Heterogeneities at Two Scales

Figure 3 shows the spatial distributions of topographic factors and the LAI on Julian day 193 at 500 m and 1 km resolutions. Although the values of these variables at two resolutions presented similarity on the whole distributions, the spatial variations at a 500 m resolution are more pronounced than those at a 1 km resolution.
The statistics of the topographic factors and the LAI on Julian day 193 at 500 m and 1 km are listed in Table 1. The altitude at a 500 m resolution varied between a minimum of 2606 m and a maximum of 4744 m; at a 1 km resolution, the altitude varied between a minimum of 2647 m and a maximum of 4671 m. The mean altitudes at the two resolutions are similar, while the range and deviation in altitude at 500 m is larger than that at a 1 km resolution. The slope had a maximum of 38° at a 500 m resolution, which varied from the 28° maximum at a 1 km resolution. The mean slope value at 500 m (18°) was greater than that at a 1 km resolution (12°), showing that there was a considerable loss of topographic information at a coarse resolution. The statistical results also show that the LAI on Julian day 193 at 1 km had a smaller deviation than that at 500 m, and the mean value of the LAI at 500 m (7.0 m2 m−2) was greater than that at a 1 km resolution (5.7 m2 m−2). In addition, the aspect at 500 m had a similar range, mean value and deviation to those at 1 km.

3.2. GPP Difference of the Two MODIS Products

To assess the GPP inconsistency at 1 km and 500 m, a pixel-by-pixel comparison between the MODIS GPP at two resolutions was implemented with a time interval of 8 days during Julian days 169–209. Figure 4a presents the relationships between MODIS products at 500 m and 1 km during the study period. From Julian day 169 to Julian day 209, the 1 km MODIS GPP explained between 24% and 47% of MODIS GPP at 500 m, with RMSEs varying between 1.34 and 2.14 gC m−2d−1. In general, the MODIS GPP at 1 km had a significant inconsistency with MODIS GPP at 500 m (R2 = 0.36, RMSE = 1.62) during the whole study period, showing that the simplification of complex terrain within each modeling pixel might cause the model results to be considerably biased.
To investigate the effects of topography on GPP, a mismatched MODIS GPP between a 1 km resolution and a 500 m resolution was analyzed. Figure 5 illustrates the relationships between topographic factors (i.e., altitude, slope) and the inconsistency in the MODIS GPP at 500 m and 1 km during Julian days 169–209. Altitude and slope were divided by 500 m and 10° intervals, respectively, and the relative error (RE) was adopted to evaluate the inconsistency in the MODIS GPP at 500 m and 1 km. It was found that the inconsistency in MODIS GPP at two resolutions increased with the increase of slope and altitude. In addition, there was not an obvious relationship between the aspect and the inconsistency in MODIS GPP at two resolutions.

3.3. Downscaled Results

To evaluate the quality of the downscaled GPP, a pixel-by-pixel comparison between the MODIS GPP at 500 m and the downscaled GPP results was implemented during Julian days 169–209 (Figure 4b). From Julian day 169 to 209, the downscaled GPP for all pixels could be explained between 85% and 91% of 500 m MODIS GPP, with RMSEs varying between 0.71 and 1.30 gC m−2d−1. The points in the scatterplots clustered around the 1:1 lines (R2 = 0.89, RMSE = 1.03), and the downscaled GPP were slightly higher than MODIS GPP at high GPP values. In general, the downscaled GPP showed a strong consistency with the MODIS GPP at 500 m.
Figure 6 illustrates the absolute differences of GPP values between the 500 m MODIS GPP and the downscaled GPP results. During Julian days 169–209, the percentages of absolute differences in ranges of [0, 1], [1, 2], [2, 3], [3, 4], [4, 5] and [5, 6] were 70.78%, 23.02%, 3.47%, 1.61%, 0.89% and 0.23%, respectively. The pixels whose absolute GPP differences between the 500 m MODIS GPP and the downscaled GPP results were in the range of [0, 2], accounting for 93.80% of the total pixels in the study area, showing that the downscaling method proposed in this study could be effectively used to obtain the GPP at 500 m resolution. Figure 7 shows the density distributions concerning the GPP differences between the 500 m MODIS GPP and the downscaled GPP results. The mean values of differences between downscaled GPP and 500 m MODIS GPP varied between 0.60 and 1.03 gC m−2d−1, with standard deviations changing between 0.47 and 1.10 gC m−2d−1, which were smaller than those between MODIS GPP at 500 m and 1 km (mean = 1.36, std. = 1.34). It was worthwhile to note that, the proposed downscaling method could effectively solve the GPP underestimation of 1 km MODIS GPP.
To illustrate the advantages and benefits of the proposed downscaling method, the downscaled results in this work, along with the kriging downscaled results, spline downscaled results and IDW downscaled results, were compared against the MODIS GPP at 500 m (Figure 4). In general, the proposed downscaled results (R2 = 0.89, RMSE = 1.03) had a stronger consistency with the 500 m MODIS GPP than the kriging downscaled results (R2 = 0.43, RMSE = 1.36), the spline downscaled results (R2 = 0.40, RMSE = 1.50) and the IDW downscaled results (R2 = 0.42, RMSE = 1.10) for all Julian days. Hence, it can be concluded that the proposed downscaling method has a potential application for GPP downscaling to obtain 500 m GPP distributions, especially in mountainous areas. The results also illustrated that vegetation heterogeneity factors and topographic factors were effective for downscaling land surface fluxes, showing good matches with other related studies [7,8,38].

4. Discussions

4.1. The Topographic Effects on GPP Inconsistency at Two Resolutions

The time series GPP estimation is very important for monitoring vegetation growth conditions. Mountainous regions typically present high spatial heterogeneity, which is characteristic of steep slopes, altitude variations and notably biological diversity [11,39,40]. The results in this study showed the topography had a significant relationship with MODIS GPP inconsistency at two resolutions, as the topography influences the surrounding environment of the carbon assimilation procedure, such as air humidity, air temperature, water redistribution and solar radiation. Generally, air temperature increases linearly with the decrease of altitude. However, in most carbon cycle models, air temperature is assumed to be homogeneous within each modeling grid. This assumption neglects the impact of inner altitude variations in each grid, bringing uncertainties to the simulation over complex terrain. Besides air temperature, solar radiation, air humidity and precipitation also face the same problems. Related studies have illustrated that the spatial variations of altitude and slope greatly influenced hydrologic procedure [13,14,41], and hence affect the carbon assimilation procedure [42,43,44]. In general, the patchy ecosystem structures in mountain areas usually follow topographic and climatic gradients, and the remarkably spatial heterogeneity has considerable effects on the spatial distributions of GPP estimation.

4.2. Evaluations of the Proposed Downscaling Method

Serving as an effective tool in the estimation of GPP, the carbon cycle model is typically executed at coarse resolutions due to the paucity of the data and computing complexities. Landscape complexities are assumed to be simplified in a carbon cycle model by ignoring spatial heterogeneities within each modeling grid. It has been illustrated that the oversimplification of landscape complexities may inevitably cause model simulations to become considerably biased [8,9,10]. In this approach, massive forcing data were required to run the carbon cycle model, including the spatial distributions of temperature, radiation, vapor pressure deficits and so on. In mountainous areas, the regional estimation of the meteorological element was a challenge, especially at fine resolution [12,13,45].
In this study, a downscaling algorithm of GPP was developed to obtain the GPP distribution at 500 m resolution. The proposed downscaling method contains two phases: a regression process and an ATPK process. The first phase was used to reduce uncertainties caused by topographic effects and vegetation heterogeneities. A regression analysis was first conducted between the GPP estimation and spatial heterogeneity factors; then, the ATPK was performed to downscale residuals from the regression process.
This study tested the thesis that the 500 m GPP could be obtained from 1 km GPP and spatial heterogeneity information, which avoided complex calculation and the limit of forcing data through the carbon cycle models. It was found that the proposed method could be effectively used to obtain the GPP at 500 m resolution, showing significant advantages over the kriging method, spline method and IDW method. Our results also illustrated that both endogenous factors (i.e., LAI) and exogenous factors (i.e., altitude, slope and aspect) were useful for GPP downscaling. Exogenous factors have significant effects on the various instantaneous land surface fluxes due to the redistribution of water and their interactions with atmospheric elements [43,46,47]. The endogenous factors can be regarded as the accumulated outcomes of the surrounding environmental conditions [48], such as topography, climate and soil property. Therefore, both exogenous and endogenous factors should be considered in the downscaling scheme to remove biases within the coarse modeling grid.

4.3. Limitations of the Current Work and Prospects for Future Studies

The proposed downscaling method in this work fully utilized valuable topographic and vegetation information from ancillary data at a 500 m spatial resolution to reduce uncertainties in spatial heterogeneities over rugged terrain. However, due to practical constraints, there were also several limitations in this study. First, the proposed downscaling algorithm was based on the MODIS GPP, DEM and LAI data, and any errors in the source data would propagate into the final downscaled results [49,50]. Second, the resampling method for 30 m GDEM data influenced the quality of DEM data at 500 m and 1 km resolutions, and hence affected the accuracy of the topographic information [51,52]. In addition, the choice of algorithm for slope and aspect also resulted in uncertainty [53,54,55], besides slightly affecting the downscaling procedure. Third, although our algorithm considered the main factors (i.e., vegetation density, altitude, slope and aspect) for GPP downscaling, several other factors (e.g., land cover, soil texture, and soil depth) were not included in the algorithm, which might result in slight uncertainties. Finally, due to the lack of reference data for validation at smaller resolution, MODIS GPP at 1 km resolution were scaled only down to a 500 m resolution. This study tested the thesis that the 500 m GPP could be obtained from 1 km GPP and spatial heterogeneity information. Based on this thesis, the relevant work of obtaining the GPP at a finer resolution (e.g., 30 m) will be presented in an additional study.
To improve the performance of the proposed downscaling method, future studies should consider more factors in the downscaling algorithm, and attention should be given to the effects of accurate ancillary data (i.e., DEM and LAI data) on the final GPP downscaled results. Moreover, the proposed downscaling method should be applied in more heterogeneous areas to test the feasibility and stability of the algorithm.

5. Conclusions

In this work, a downscaling algorithm was developed to retrieve the GPP at 500 m resolution in a time series over rugged terrain, which considered the effects of spatial heterogeneity on carbon flux simulations. The GPP product at 1 km served as the primary dataset of the downscaling algorithm, and the 500 m GPP product was used as the reference dataset to evaluate the downscaled GPP results. It was found that the proposed downscaling method had advantages over the kriging method, spline method and IDW method. Our results also showed that vegetation heterogeneity factors (i.e., LAI) and topographic factors (i.e., altitude, slope and aspect) were useful for GPP downscaling.
Our study developed a reliable downscaling method to derive a 500 m resolution GPP for time series over mountainous areas, which is an important issue when evaluating the role of ecosystems in the global carbon cycle. Furthermore, our proposed downscaling algorithm could provide a reference when considering the effects of spatial heterogeneity on carbon flux simulations and retrieving other fine resolution parameters (e.g., net primary productivity and evaporation) in topographically complex terrains.

Acknowledgments

The research was supported from the National Key Research and Development Program of China (Grant No. 2016YFA0600103), the National Natural Science Foundation of China (Grant No. 41631180, 41571373, 41671376). We are grateful to all the data providers and the anonymous reviewers for their valuable suggestions.

Author Contributions

Xinyao Xie, Ainong Li and Huaan Jin conceived and designed experiments; Xinyao Xie performed the experiments, analyzed the data and wrote the paper; Gaofei Yin provided guidance and helped write the codes of algorithm; Jinhu Bian downloaded the data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Li, A.; Bian, J.; Lei, G.; Huang, C. Estimating the maximal light use efficiency for different vegetation through the casa model combined with time-series remote sensing data and ground measurements. Remote Sens. 2012, 4, 3857–3876. [Google Scholar] [CrossRef]
  2. Zhou, Y.; Hilker, T.; Ju, W.; Coops, N.C.; Black, T.A.; Chen, J.M.; Wu, X. Modeling gross primary production for sunlit and shaded canopies across an evergreen and a deciduous site in canada. IEEE Trans. Geosci. Remote Sens. 2017, 55, 1859–1873. [Google Scholar] [CrossRef]
  3. Liu, L.; Guan, L.; Liu, X. Directly estimating diurnal changes in gpp for c3 and c4 crops using far-red sun-induced chlorophyll fluorescence. Agric. For. Meteorol. 2017, 232, 1–9. [Google Scholar] [CrossRef]
  4. Ju, W.; Chen, J.M.; Black, T.A.; Barr, A.G.; Liu, J.; Chen, B. Modelling multi-year coupled carbon and water fluxes in a boreal aspen forest. Agric. For. Meteorol. 2006, 140, 136–151. [Google Scholar] [CrossRef]
  5. Alemu, G.; Chen, J.M.; Price, D.T.; Kurz, W.A.; Liu, J.; Céline, B.; Hember, R.A.; Wu, C.; Chang, K.H. Improved assessment of gross and net primary productivity of canada’s landmass. J. Geophys. Res. Atmos. 2013, 99, 1089–1104. [Google Scholar]
  6. Xie, X.; Li, A.; Jin, H. The simulation models of the forest carbon cycle on a large scale : A review. Acta Ecol. Sin. 2018, 38, 1–14. [Google Scholar]
  7. Chen, J.M.; Chen, X.; Ju, W. Effects of vegetation heterogeneity and surface topography on spatial scaling of net primary productivity. Biogeosciences 2013, 10, 4879–4896. [Google Scholar] [CrossRef]
  8. El Maayar, M.; Chen, J.M. Spatial scaling of evapotranspiration as affected by heterogeneities in vegetation, topography, and soil texture. Remote Sens. Environ. 2006, 102, 33–51. [Google Scholar] [CrossRef]
  9. Mitchell, S.W.; Csillag, F.; Tague, C. Impacts of spatial partitioning in hydroecological models: Predicting grassland productivity with rhessys. Trans. GIS 2005, 9, 421–442. [Google Scholar] [CrossRef]
  10. MacKay, D.S.; Ahl, D.E.; Ewers, B.E.; Gower, S.T.; Burrows, S.N.; Samanta, S.; Davis, K.J. Effects of aggregated classifications of forest composition on estimates of evapotranspiration in a northern wisconsin forest. Glob. Chang. Biol. 2002, 8, 1253–1265. [Google Scholar] [CrossRef]
  11. Jin, H.; Li, A.; Bian, J.; Nan, X.; Zhao, W.; Zhang, Z.; Yin, G. Intercomparison and validation of modis and glass leaf area index (lai) products over mountain areas: A case study in southwestern china. Int. J. Appl. Earth Observ. Geoinf. 2017, 55, 52–67. [Google Scholar] [CrossRef]
  12. Zhao, W.; Li, A.N. A review on land surface processes modelling over complex terrain. Adv. Meteorol. 2015, 2015, 607181. [Google Scholar] [CrossRef]
  13. Immerzeel, W.W.; Petersen, L.; Ragettli, S.; Pellicciotti, F. The importance of observed gradients of air temperature and precipitation for modeling runoff from a glacierized watershed in the nepalese himalayas. Water Resour. Res. 2014, 50, 2212–2226. [Google Scholar] [CrossRef]
  14. Mizukami, N.; Clark, M.P.; Slater, A.G.; Brekke, L.D.; Elsner, M.M.; Arnold, J.R.; Gangopadhyay, S. Hydrologic implications of different large-scale meteorological model forcing datasets in mountainous regions. J. Hydrometeorol. 2014, 15, 474–488. [Google Scholar] [CrossRef]
  15. Helbig, N.; Loewe, H. Shortwave radiation parameterization scheme for subgrid topography. J. Geophys. Res. Atmos. 2012, 117. [Google Scholar] [CrossRef]
  16. Li, A.; Bian, J.; Jin, H.; Zhao, W.; Zhang, Z.; Nan, X.; Yin, G.; Lei, G. Mountain Remote Sensing; Science Press: Beijing, China, 2016. [Google Scholar]
  17. Chen, J.M. Spatial scaling of a remotely sensed surface parameter by contexture. Remote Sens. Environ. 1999, 69, 30–42. [Google Scholar] [CrossRef]
  18. Liang, S.; Li, X.; Wang, J. Advanced Remote Sensing; Academic Press: Cambridge, MA, USA, 2012. [Google Scholar]
  19. Hong, S.-H.; Hendrickx, J.M.H.; Borchers, B. Up-scaling of sebal derived evapotranspiration maps from landsat (30 m) to modis (250 m) scale. J. Hydrol. 2009, 370, 122–138. [Google Scholar] [CrossRef]
  20. Thomas, A. Overview of the geoecology of the gongga shan range, sichuan province, china. Mt. Res. Dev. 1999, 19, 17–30. [Google Scholar] [CrossRef]
  21. Zhu, W.; Fu, X.; Feng, X.; Lu, J.Y. Annual time-series analyses of total gaseous mercury measurement and its impact factors on the gongga mountains in the southeastern fringe of the qinghai-tibetan plateau. J. Mt. Sci. Engl. 2008, 5, 17–31. [Google Scholar] [CrossRef]
  22. Running, S.W.; Nemani, R.R.; Heinsch, F.A.; Zhao, M.S.; Reeves, M.; Hashimoto, H. A continuous satellite-derived measure of global terrestrial primary production. Bioscience 2004, 54, 547–560. [Google Scholar] [CrossRef]
  23. Heinsch, F.A.; Zhao, M.; Running, S.W.; Kimball, J.S.; Nemani, R.R.; Davis, K.J.; Bolstad, P.V.; Cook, B.D.; Desai, A.R.; Ricciuto, D.M.; et al. Evaluation of remote sensing based terrestrial productivity from modis using regional tower eddy flux network observations. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1908–1925. [Google Scholar] [CrossRef]
  24. Jung, M.; Reichstein, M.; Margolis, H.A.; Cescatti, A.; Richardson, A.D.; Arain, M.A.; Arneth, A.; Bernhofer, C.; Bonal, D.; Chen, J.; et al. Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations. J. Geophys. Res. Biogeo 2011, 116. [Google Scholar] [CrossRef]
  25. Turner, D.P.; Ritts, W.D.; Cohen, W.B.; Gower, S.T.; Running, S.W.; Zhao, M.S.; Costa, M.H.; Kirschbaum, A.A.; Ham, J.M.; Saleska, S.R.; et al. Evaluation of modis npp and gpp products across multiple biomes. Remote Sens. Environ. 2006, 102, 282–292. [Google Scholar] [CrossRef]
  26. Myneni, R.B.; Hoffman, S.; Knyazikhin, Y.; Privette, J.L.; Glassy, J.; Tian, Y.; Wang, Y.; Song, X.; Zhang, Y.; Smith, G.R.; et al. Global products of vegetation leaf area and fraction absorbed par from year one of modis data. Remote Sens. Environ. 2002, 83, 214–231. [Google Scholar] [CrossRef]
  27. Knyazikhin, Y.; Martonchik, J.V.; Myneni, R.B.; Diner, D.J.; Running, S.W. Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from modis and misr data. J. Geophys. Res. Atmos. 1998, 103, 32257–32275. [Google Scholar] [CrossRef]
  28. Zhang, C.; Liu, Q.; Liu, G.; Ding, S.; Li, H.; Dong, J. Elevation quality evaluation of aster gdem data in china’s eastern region. J. Geomat. Sci. Technol. 2014, 31, 67–72. [Google Scholar]
  29. Zhang, C.; Liu, Q.; Liu, G.; Ding, S.; Wang, P.; Dong, J. Data processing and application progress of srtm3 and aster gdem. Geogr. Geo-inf. Sci. 2012, 28, 29–34. [Google Scholar]
  30. Frey, H.; Paul, F. On the suitability of the srtm dem and aster gdem for the compilation of topographic parameters in glacier inventories. Int. J. Appl. Earth Obs. Geoinf. 2012, 18, 480–490. [Google Scholar] [CrossRef]
  31. Dozier, J.; Frew, J. Rapid calculation of terrain parameters for radiation modeling from digital elevation data. IEEE Trans. Geosci. Remote Sens. 1990, 28, 963–969. [Google Scholar] [CrossRef]
  32. Wang, Q.; Shi, W.; Atkinson, P.M.; Zhao, Y. Downscaling modis images with area-to-point regression kriging. Remote Sens. Environ. 2015, 166, 191–204. [Google Scholar] [CrossRef]
  33. Gao, F.; Kustas, W.P.; Anderson, M.C. A data mining approach for sharpening thermal satellite imagery over land. Remote Sens. 2012, 4, 3287–3319. [Google Scholar] [CrossRef]
  34. Wang, Q.; Shi, W.; Atkinson, P.M. Area-to-point regression kriging for pan-sharpening. ISPRS J. Photogramm. Remote Sens. 2016, 114, 151–165. [Google Scholar] [CrossRef]
  35. Knotters, M.; Brus, D.J.; Voshaar, J.H.O. A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations. Geoderma 1995, 67, 227–246. [Google Scholar] [CrossRef]
  36. Tait, A.; Henderson, R.; Turner, R.; Zheng, X. Thin plate smoothing spline interpolation of daily rainfall for new zealand using a climatological rainfall surface. Int. J. Climatol. 2006, 26, 2097–2115. [Google Scholar] [CrossRef]
  37. Watson, D.F.; Philip, G.M. A refinement of inverse distance weighted interpolation. Geo-Processing 1985, 2, 315–327. [Google Scholar]
  38. Simic, A.; Chen, J.M.; Liu, J.; Csillag, F. Spatial scaling of net primary productivity using subpixel information. Remote Sens. Environ. 2004, 93, 246–258. [Google Scholar] [CrossRef]
  39. Jin, H.; Li, A.; Bian, J.; Zhang, Z.; Huang, C.; Li, M. Validation of global land surface satellite (glass) downward shortwave radiation product in the rugged surface. J. Mt. Sci. Engl. 2013, 10, 812–823. [Google Scholar] [CrossRef]
  40. Jin, H.A.; Li, A.N.; Bian, J.H. Comparative analysis of hj-1, spot, and tm data for leaf area index estimation in a mountainous area. In Proceedings of the 2013 IEEE International Geoscience and Remote Sensing Symposium, Melbourne, Australia, 21–26 July 2013; pp. 2782–2785. [Google Scholar]
  41. Bonan, G.B. The land surface climatology of the ncar land surface model coupled to the ncar community climate model. J. Clim. 1998, 11, 1307–1326. [Google Scholar] [CrossRef]
  42. Qiu, Y.; Fu, B.J.; Wang, J.; Chen, L.D. Soil moisture variation in relation to topography and land use in a hillslope catchment of the loess plateau, china. J. Hydrol. 2001, 240, 243–263. [Google Scholar] [CrossRef]
  43. Govind, A.; Chen, J.M.; Margolis, H.; Ju, W.; Sonnentag, O.; Giasson, M.-A. A spatially explicit hydro-ecological modeling framework (beps-terrainlab v2.0): Model description and test in a boreal ecosystem in eastern north america. J. Hydrol. 2009, 367, 200–216. [Google Scholar] [CrossRef]
  44. Sabetraftar, K.; Mackey, B.; Croke, B. Sensitivity of modelled gross primary productivity to topographic effects on surface radiation: A case study in the cotter river catchment, australia. Ecol. Model. 2011, 222, 795–803. [Google Scholar] [CrossRef]
  45. Konzelmann, T.; Calanca, P.; Muller, G.; Menzel, L.; Lang, H. Energy balance and evapotranspiration in a high mountain area during summer. J. Appl. Meteorol. 1997, 36, 966–973. [Google Scholar] [CrossRef]
  46. Huang, W.; Zhang, L.; Furumi, S.; Muramatsu, K.; Daigo, M.; Li, P. Topographic effects on estimating net primary productivity of green coniferous forest in complex terrain using landsat data: A case study of yoshino mountain, japan. Int. J. Remote Sens. 2010, 31, 2941–2957. [Google Scholar] [CrossRef]
  47. Gelybo, G.; Barcza, Z.; Kern, A.; Kljun, N. Effect of spatial heterogeneity on the validation of remote sensing based gpp estimations. Agric. For. Meteorol. 2013, 174, 43–53. [Google Scholar] [CrossRef]
  48. Chen, J.M.; Mo, G.; Pisek, J.; Liu, J.; Deng, F.; Ishizawa, M.; Chan, D. Effects of foliage clumping on the estimation of global terrestrial gross primary productivity. Glob. Biogeochem. Cycle 2012, 26. [Google Scholar] [CrossRef]
  49. Lin, S.; Jing, C.; Coles, N.A.; Chaplot, V.; Moore, N.J.; Wu, J. Evaluating dem source and resolution uncertainties in the soil and water assessment tool. Stoch. Environ. Res. Risk Assess. 2013, 27, 209–221. [Google Scholar] [CrossRef]
  50. Xu, F.; Dong, G.; Wang, Q.; Liu, L.; Yu, W.; Men, C.; Liu, R. Impacts of dem uncertainties on critical source areas identification for non-point source pollution control based on swat model. J. Hydrol. 2016, 540, 355–367. [Google Scholar] [CrossRef]
  51. Grohmann, C.H. Effects of spatial resolution on slope and aspect derivation for regional-scale analysis. Comput. Geosci. 2015, 77, 111–117. [Google Scholar] [CrossRef]
  52. Wechsler, S.P.; Kroll, C.N. Quantifying dem uncertainty and its effect on topographic parameters. Photogramm. Eng. Remote Sens. 2006, 72, 1081–1090. [Google Scholar] [CrossRef]
  53. Oksanen, J.; Sarjakoski, T. Error propagation of dem-based surface derivatives. Comput. Geosci. 2005, 31, 1015–1027. [Google Scholar] [CrossRef]
  54. Zhou, Q.M.; Liu, X.J. Analysis of errors of derived slope and aspect related to dem data properties. Comput. Geosci. 2004, 30, 369–378. [Google Scholar] [CrossRef]
  55. Liu, H.H.; Kiesel, J.; Hormann, G.; Fohrer, N. Effects of dem horizontal resolution and methods on calculating the slope length factor in gently rolling landscapes. Catena 2011, 87, 368–375. [Google Scholar] [CrossRef]
Figure 1. Study area.
Figure 1. Study area.
Remotesensing 10 00647 g001
Figure 2. Flow chart of the proposed downscaling algorithm.
Figure 2. Flow chart of the proposed downscaling algorithm.
Remotesensing 10 00647 g002
Figure 3. Spatial distributions of altitude (a); slope (b); aspect (c) and LAI (d) at 500 m and 1 km.
Figure 3. Spatial distributions of altitude (a); slope (b); aspect (c) and LAI (d) at 500 m and 1 km.
Remotesensing 10 00647 g003
Figure 4. Density scatterplot between MODIS GPP and downscaled GPP at 500 m during Julian days 169–209. (a) represents the comparison between MODIS GPP products; and (be) represent the comparison between downscaled GPP and MODIS GPP at 500 m. The solid lines are the regression lines, while the dashed lines are the 1:1 lines. The green and red represent the low-density and high-density areas, respectively.
Figure 4. Density scatterplot between MODIS GPP and downscaled GPP at 500 m during Julian days 169–209. (a) represents the comparison between MODIS GPP products; and (be) represent the comparison between downscaled GPP and MODIS GPP at 500 m. The solid lines are the regression lines, while the dashed lines are the 1:1 lines. The green and red represent the low-density and high-density areas, respectively.
Remotesensing 10 00647 g004
Figure 5. Relationships between topographic factors and inconsistency in the MODIS GPP at 500 m and 1 km during Julian days 169–209.
Figure 5. Relationships between topographic factors and inconsistency in the MODIS GPP at 500 m and 1 km during Julian days 169–209.
Remotesensing 10 00647 g005
Figure 6. The spatial distribution of absolute GPP differences between the 500 m MODIS GPP and the downscaled GPP results for Julian days 169–209. Figures (af) indicate Julian days 169–209, respectively.
Figure 6. The spatial distribution of absolute GPP differences between the 500 m MODIS GPP and the downscaled GPP results for Julian days 169–209. Figures (af) indicate Julian days 169–209, respectively.
Remotesensing 10 00647 g006
Figure 7. The density distributions concerning the GPP differences between the 500 m MODIS GPP and the downscaled GPP results for Julian days 169–209. Figures (af) indicate Julian days 169–209, respectively.
Figure 7. The density distributions concerning the GPP differences between the 500 m MODIS GPP and the downscaled GPP results for Julian days 169–209. Figures (af) indicate Julian days 169–209, respectively.
Remotesensing 10 00647 g007
Table 1. Comparisons of topographic factors and at 500 m and 1 km.
Table 1. Comparisons of topographic factors and at 500 m and 1 km.
Pixel Resolution (m)MinMaxMeanSTD
5001000500100050010005001000
Altitude (m)260626474744467137843787501489
Slope (deg)003828181286
Aspect (deg)00360360178177101100
LAI (m2 m−2)007.05.71.881.581.931.61

Share and Cite

MDPI and ACS Style

Xie, X.; Li, A.; Jin, H.; Yin, G.; Bian, J. Spatial Downscaling of Gross Primary Productivity Using Topographic and Vegetation Heterogeneity Information: A Case Study in the Gongga Mountain Region of China. Remote Sens. 2018, 10, 647. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10040647

AMA Style

Xie X, Li A, Jin H, Yin G, Bian J. Spatial Downscaling of Gross Primary Productivity Using Topographic and Vegetation Heterogeneity Information: A Case Study in the Gongga Mountain Region of China. Remote Sensing. 2018; 10(4):647. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10040647

Chicago/Turabian Style

Xie, Xinyao, Ainong Li, Huaan Jin, Gaofei Yin, and Jinhu Bian. 2018. "Spatial Downscaling of Gross Primary Productivity Using Topographic and Vegetation Heterogeneity Information: A Case Study in the Gongga Mountain Region of China" Remote Sensing 10, no. 4: 647. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10040647

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop