Prediction of High-Quality MODIS-NPP Product Data

Abstract

Net primary productivity (NPP) is a key vegetation parameter and ecological indicator for tracking natural environmental change. High-quality Moderate Resolution Imaging Spectroradiometer Net primary productivity (MODIS-NPP) products are critical for assuring the scientific rigor of NPP analyses. However, obtaining high-quality MODIS-NPP products consistently is challenged by factors such as cloud contamination, heavy aerosol pollution, and atmospheric variability. This paper proposes a method combining the discrete wavelet transform (DWT) with an extended Kalman filter (EKF) for generating high-quality MODIS-NPP data. In this method, the DWT is used to remove noise in the original MODIS-NPP data, and the EKF is applied to the de-noised images. The de-noised images are modeled as a triply modulated cosine function that predicts the NPP data values when excessive cloudiness is present. This study was conducted in South China. By comparing measured NPP data to original MODIS-NPP and NPP estimates derived from combining the DWT and EKF, we found that the accuracy of the NPP estimates was significantly improved. The MODIS-NPP estimates had a mean relative error (RE) of 13.96% and relative root mean square error (rRMSE) of 15.67%, while the original MODIS-NPP had a mean RE of 23.58% and an rRMSE of 24.98%. The method combining DWT and EKF provides a feasible approach for generating new, high-quality NPP data in the absence of high-quality original MODIS-NPP data.

1. Introduction

Net primary productivity (NPP) is a key variable in the carbon exchange between the biosphere and the atmosphere [1]. NPP is defined as the amount of above- and below-ground organic matter accumulated in vegetation per unit of area during a certain time period (usually a year), representing the productivity of an ecosystem [2]. NPP with 1 km spatial resolution and an 8 day interval has been produced operationally with the Moderate Resolution Imaging Spectroradiometer Net Primary Productivity product (MOD17) algorithm based on observations from the MODIS sensor (National Aeronautics and Space Administration, America) [3]. These MODIS-NPP data provide consistent spatial and temporal measures of crop yield, range forage, forest production, and other socio-economically significant products related to vegetation growth [4,5]. Thus, accurate NPP estimation is critical. The MOD17 algorithm was developed as a result of past achievements and lessons learned from a general ecosystem model, Biome-BGC [5]. A detailed description of the algorithm can be found elsewhere [6,7]. Although the MOD17 algorithm has undergone several improvements, MODIS-NPP satellite data still include various sources of noise, such as clouds, atmospheric variability, and bidirectional reflectance distribution factors [8,9]. These distribution factors can affect the monitoring of crop yield, range forage, and forests. Therefore, it is crucial to reduce noise in images and obtain high-quality NPP data.

Several methods for reducing noise and constructing high-quality MODIS datasets have been proposed, applied, and evaluated in recent years [10,11,12,13,14,15,16]. Fast Fourier transform (FFT) and wavelet transform (WT) approaches are often used when evaluating MODIS time-series data [17,18,19,20]. Several mathematical filters have been applied to reduce noise and obtain high-quality time-series of remote sensing datasets [21,22]. Recently, data assimilation methods have also been used to reconstruct remote-sensing data dynamically. Gu et al. [23] adopted an optimal interpolation method to generate appropriate Moderate Resolution Imaging Spectroradiometer leaf area index (MODIS-LAI) data. Gu et al. [24] developed a simplified data assimilation method based on Normalized Difference Vegetation Index (NDVI) quality assurance (QA) datasets to reconstruct high-quality, spatially, and temporally continuous MODIS NDVI products. However, these methods also suffer from drawbacks. For example, they can only produce a high-quality output when the data contain little noise, so they are challenged when used with very cloudy images.

The objective of this study was to generate high-quality MODIS-NPP data by combining discrete wavelet transform (DWT) with an extended Kalman filter (EKF). Taking an area in South China as an example, we used the DWT to reduce cloud and atmospheric variability contamination in MODIS-NPP data. Then, the EKF was applied to the de-noised MODIS-NPP data to predict NPP data when the original data was not valid.

2. Materials

The study area was located in South China, and covered approximately 65,536 km² (256 km × 256 km) within N22°36’30″–N24°44’0″ and E112°53’21.25"–E113°22’13.28″ (Figure 1). The area has a humid, subtropical, monsoon climate characterized by warm winters, hot summers, little frost and snow, and sufficient rain and sunshine, with an annual mean temperature of 22 °C. The annual average precipitation is 1400–2000 mm. The main vegetation types in the area are subtropical evergreen broadleaf forests, evergreen needleleaf forests, savannas, and common crops such as rice, peanut, and rape (Figure 1).

NPP measurements were taken from a total of 20 sample plots located in evergreen broadleaf and evergreen needle leaf forests within the study area in October 2013 (Figure 1). Each of the sample plots had an area of 1 km × 1 km, matching the spatial resolution of the MODIS products. Within each sample plot, the tree species, height, and diameter at breast height (DBH), with breast height being 130 cm, and age class (young, middle age, and mature) were measured. Estimates of above-ground tree biomass (AGB) were predicted with empirical regression models using given the following general formula:

$${\mathrm{AGB}=a\times \mathrm{DBH}}^b,$$

(1)

where $a$ and $b$ are empirical coefficients that vary by tree species [25,26,27,28]. Estimates of tree biomass were calculated using the ratio of AGB to the below-ground tree biomass [29]. Estimates of tree biomass in the sample plots were obtained by adding the tree biomass values within each plot. The value of NPP for each sample plot was estimated using the plot biomass estimate divided by the average age, in years, of the trees within each plot [30]. The NPP within 1 km was an average of five sample NPPs. The average age of trees for each plot was estimated by calculating a weighted mean of the average age classes and the number of trees within each age class. According to previous studies [31], the six representative sample plots, including Schima, Pinus massoniana, Acacia confusa, Camphora, Ficus, and Eucalyptus, were selected to collect the measured NPP data, from 18 October 2013 to 27 December 2015, as shown in Figure 1.

To evaluate the performance of the method proposed here, we used 8 day composites of MOD17 data with a spatial resolution of 1 km from 9 January 2013 to 18 December 2015. The original MODIS-NPP data composites with heavy cloud cover do not provide useful data for predicting high-quality MODIS-NPP data during periods of active phenological change. Thus, the time-distance weight function was used to obtain these NPP data using NPP composites with no heavy clouds at adjoining time periods. The time-distance weight function can be expressed as

$${\mathrm{NPP}}_i{=c\times \mathrm{NPP}}_k{+d\times \mathrm{NPP}}_j,$$

(2)

where ${\mathrm{NPP}}_i$ is the ith day NPP, ${\mathrm{NPP}}_k$ represents the original MODIS-NPP value with no heavy clouds on the closest kth day preceding the ith day, and ${\mathrm{NPP}}_j$ represents the original MODIS-NPP value with no heavy clouds on the closest jth day after the ith day. The time distance weight, equal to $(i-k)/(j-k)\mathrm{is}c$, and $d$ is the time distance weight, equal to $(j-i)/(j-k).$

3. Methods

3.1. The Discrete Wavelet Transform for De-Noising MODIS-NPP Images

The purpose of de-noising original MODIS-NPP data is to reduce the noise in the DWT domain while preserving the NPP characteristics. The $N$ NPP values of a noisy image are represented as a two-dimensional matrix:

$$\left[y_{ij}\right],i=1\cdots N,j=1\cdots M,$$

(3)

where $N$ and $M$ denote the number of rows and columns in the MODIS-NPP image. The noisy version of the image is modelled as [32]

$$y_{ij}{=x}_{ij}{+n}_{ij}i=1\cdots N,j=1\cdots M,$$

(4)

where $\left[n_{ij}\right]$ are $iid$ as $N\left({0,\sigma }^2\right)$ and $x_{ij}$ denotes the de-noised data value.

Wavelet de-noising involved three main steps: calculation of the DWT of the image, completion of a nonlinear thresholding step, and completion of a linear inverse DWT to get the de-noised estimate. One of the most important steps in the de-noising development process was finding the optimal threshold. Encouraged by results in recent studies [33,34,35,36,37], we tested six thresholding functions (Figure 2) to de-noise the original MODIS-NPP images in this study.

Taking a MODIS-NPP image from 9 January 2013 as an example, the signal-to-noise ratios (SNRs) and root mean square difference (RMSDs) from de-noising the image using the six thresholds varied from 9.02 to 17.20 Mg ha⁻¹ a⁻¹ and from 0.38 to 0.99, respectively (

Table 1. Signal-to-noise ratios (SNRs) and root mean square difference (RMSDs) for the different threshold functions.

Threshold function	SNR	RMSD
ST	12.11	0.69
HT	9.02	0.99
HST	10.86	0.82
ISHST	9.95	0.89
ISHT	17.20	0.38
ISIT	11.38	0.75

). The highest SNR and lowest RMSE were obtained using the improved soft and hard threshold (ISHT) threshold function, suggesting that this was the most suitable thresholding function for de-noising the images. Thus, the ISHT threshold function was selected as the optimal threshold in this study.

3.2. The Extended Kalman Filter for Predicting High-Quality MODIS-NPP Data

The formula to predict MODIS-NPP time series in a given pixel is written as [38]

$$y_k{=\mu }_k{+\alpha }_k\cos \left({\omega }_k{+\Phi }_k\right){+v}_k$$

(5)

where the subindex k indicates the time step, $y_k$ represents the observed value of the predicted MODIS-NPP time series, $v_k$ is the noise sample, ${\omega }_k$ denotes the angular frequency, ${\alpha }_k$ is the amplitude, ${\mu }_k$ is the mean, and ${\Phi }_k$ is the phase. The angular frequency can be calculated as ${\omega }_k=2\pi f$, where $f$ is determined by the annual vegetation growth cycle. Because MODIS-NPP images were calculated based on the MODIS combined Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation (FPAR) products (MOD15A2 products), $f$ was set at 8/365. The values of ${\mu }_k$, ${\alpha }_k$ and ${\Phi }_k$ are functions of time, and must be estimated given $y_k$ for k∈1,...,N. The estimation of these parameters is nontrivial and requires a non-linear estimator. Thus, the EKF was driven to estimate the three parameters for every increment of k [38].

The EKF is considered the standard framework to solve estimation problems in navigation [39]. The EKF prediction formula can be obtained by linearizing the process model. The observation model at time $k+1$ is given by [39]

$$y_{k+1}=H·x_{k+1}{+v}_k,$$

(6)

where $y_{k+1}$ is the estimated NPP and $x_{k+1}$ denotes the a priori estimate, defined as $x_{k+1}={\left[{\mu }_{k+1} {\alpha }_{k+1} {\Phi }_{k+1}\right]}^T$. $H$ is an observation operator used to relate model state variables to observations (it does not need be linear), and $v_k$ is the zero-mean Gaussian observation noise $v_k~N\left(O,R\right)$. The discrete-time prediction EKF algorithm has two main stages, namely a prediction and an update (filtering) step as follows [39]:

Prediction equation

The predicted state estimate is defined as:

$$x_{k+1}=F·{\widehat{x}}_k+w_k.$$

(7)

The predicted covariance estimate is defined as

$$P_{k+1}=\nabla F_x·{\widehat{P}}_k\nabla F_x^T+Q_k,$$

(8)

where $F$ denotes the state transition matrix, $P_{k+1}$ is the a priori estimate error covariance matrix, ${\widehat{x}}_k$ represents the a posteriori estimate, ${\widehat{P}}_k$ is the a posteriori error covariance matrix, $w_k$ is the zero-mean Gaussian process noise $w_k~N\left({O,Q}_k\right)$, $\nabla F_x$ is the Jacobian of function $F$ with respect to x evaluated at ${\widehat{x}}_k$, and $Q_k$ is the covariance matrix of the process noises.

Update equation

The updated state estimate is defined as

$${\widehat{x}}_{k+1}=x_k+K_{g\left(k\right)}\left(y_k-H_k{x}_k\right).$$

(9)

The near-optimal Kalman gain is defined as

$$K_g\left(k\right){=P}_k\nabla H_k^T{\left[\nabla H_k{P}_k\nabla H_k^T{+R}_k\right]}^{-1}.$$

(10)

The updated covariance estimate is defined as

$${\widehat{P}}_{k+1}=\left[{I-K}_g\left(k\right)\nabla H_k\right]P_k,$$

(11)

where $K_g\left(k\right)$ is the Kalman gain matrix at instant $k$, $\nabla H_k$ is the Jacobian of function $H$ evaluated at $x_k$, $R_k$ is the covariance matrix of the measurement noises, and $I$ is the identity matrix.

The values of ${\mu }_k,$ ${\alpha }_k,$ and ${\Phi }_k$ are functions of time and must be estimated given $y_k$ for k∈1. This study used 137 time-series images so that M = 137. The a priori initial state parameters (${\mu }_1{\alpha }_1{,\Phi }_1$) in state vector ($x_1$) were set as

$${\mu }_1={{\displaystyle \sum }}_{i=1}^N\frac{{\mathrm{NPP}}_i}{N},$$

(12)

$$a_1=\frac{{\mathrm{NPP}}_{\max }{-\mathrm{NPP}}_{\min }}{2},$$

(13)

where $N$ is the number of remote-sensing pixels. ${\mathrm{NPP}}_{\max }$ and ${\mathrm{NPP}}_{\min }$ are the maximum and minimum NPP values in the studied area, respectively. An a priori initial estimate ${\Phi }_1{=0}^{°}$ was chosen assuming that the start date of 18 January 2013 represented a stage after vegetation was exchanged. Regarding parameters R, P₁, and Q₁, their values were first taken from the literature [39] and then tested and adjusted using the EKF model and the layout of NPP samples. Thus, they were finally set as $R=0.01$, $P_1=\mathrm{diag}\left[0.015,0.015,0.02963\right]$, $Q_1=\mathrm{diag}\left[0.01,0.001,0.03441\right]{,\mathrm{and}x}_1=\left[4.5612,2.98425,0\right]$.

3.3. Evaluation of the NPP Accuracy

The accuracy of the original MODIS-NPP data, de-noised MODIS-NPP data, and predicted MODIS-NPP data was evaluated by comparing the NPP in the 20 sample plots (Section 2) with the NPP in the pixels where the plots were located. The relative error (RE) was calculated as follows for each of the 20 plots in the dates when the ground measurements were performed:

$$\mathrm{RE}=100\times \frac{|\widehat{y}-y|}{y},$$

(14)

where $\widehat{y}$ is the MODIS-NPP and $y$ is the plot NPP. For the six representative plots, the RE was also analyzed across time. Then, the relative root mean square error (rRMSE) for the 20 plots as a whole was also calculated:

$$\mathrm{rRMSE}=100\times \frac{\sqrt{\frac{{\sum }_{i=1}^n{\left({\widehat{y}}_i-y_i\right)}^2}{n}}}{\overline{y}},$$

(15)

where ${\widehat{y}}_i$ is the MODIS-NPP, $y_i$ is the plot NPP, $\overline{y}$ is mean value of the plot NPP and i indicates the plot.

The NPP products elaborated in this study (de-noised MODIS-NPP data and predicted MODIS-NPP data) were also compared to the original MODIS-NPP data by means of the root mean square difference (RMSD) and, for the specific case of the de-noised MODIS-NPP data, through the signal-to-noise ratio (SNR):

$$\mathrm{RMSD}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{\left({\widehat{y}}_i-y_i\right)}^2}}{n}},$$

(16)

where ${\widehat{y}}_i$ is the NPP product elaborated in this study, and $y_i$ is the original MODIS-NPP.

$${\mathrm{SNR}=10*\log }_{10}\left[\frac{{{\displaystyle \sum }}_{i=1}^M{{\displaystyle \sum }}_{j=1}^N{g(i,j)}^2}{{{\displaystyle \sum }}_{i=1}^M{{\displaystyle \sum }}_{j=1}^N{\left[g\left(i,j\right)-f(i,j)\right]}^2}\right],$$

(17)

where, $g(i,j)$ represents original MODIS-NPP,$f(i,j)$ is de-noised MODIS-NPP.

4. Results

4.1. Image De-Noising

In this study, we used the DWT algorithm to reduce the noise in the MODIS-NPP images and improve the data quality of MODIS imagery. The wavelet decomposition of an image by DWT is shown in Figure 3.

The wavelet decomposition of a MODIS image using DWT (Figure 3) was completed as follows [31]: At the first level of decomposition, the image was divided into four sub-bands, $H{H}_k,H{L}_k,L{H}_k,$ and $L{L}_{k,}k=1\cdots J$, where $k$ is the scale, with $J$ being the largest (or coarsest) scale in the decomposition. The $H{H}_k$ sub-band represented the diagonal structures of the image, the $H{L}_k$ sub-band represented the horizontal features, and the $L{H}_k$ sub-band gave the vertical details. The $L{L}_k$ sub-band was the coarser resolution residual consisting of lower frequency components and was further split into higher levels of decomposition. According to the ISHT threshold function, the DWT was driven to acquire de-noised MODIS-NPP images. For visual comparison, we selected two MODIS-NPP images from 1 January 2014 and 27 December 2015 (Figure 4) with different amounts of noise. The 27 December 2015 image was heavily contaminated by clouds, while the 1 January 2014 image only contained a little noise.

The white speckles in the original image from MOD17 on 1 January 2014 were notably suppressed in the de-noised results, suggesting that the DWT could reduce noise effectively. While comparing the original image with heavy cloud noise from 27 December 2015, the results of removing noise from the original image had little change, which indicated the de-noised (post DWT) results could not remove heavy clouds. To quantify the de-noising performance, RMSD and SNR were calculated between original the MODIS-NPP images and the de-noised ones. For 1 January 2014 and 27 December 2015, the respective measures were as follows: RMSD 0.36 Mg ha⁻¹ a⁻¹ and SNR 17.42, and RMSD 0 Mg ha⁻¹ a⁻¹ and SNR ∞. Thus, the ISHT thresholding method achieved good speckle suppression along with improving the 1 January 2014 NPP image quality, while the values of SNR and RMSD were ineffective for MODIS-NPP images from 27 December 2015 with heavy cloud contamination.

4.2. Extended Kalman Filter Prediction

De-noised MODIS-NPP data were provided as input to the EKF model, and the resulting NPP data for the representative sample plots in the study area were plotted along with the original MODIS-NPP data, de-noised MODIS-NPP data, and measured NPP data (Figure 5).

Figure 5a shows that the predicted data (in blue) were close to the measured data (in green), while the original data (in black) were generally much lower than the others. Figure 5b shows that REs between predicted MODIS-NPP data and measured NPP data were smaller than REs between the original MODIS-NPP data and measured NPP data and between the de-noised MODIS-NPP data and measured NPP data. These lower REs between the predicted MODIS-NPP data and measured NPP data suggest that the proposed method has a greater potential to estimate NPP data with high accuracy compared to the results of DWT alone.

Comparing the predicted MODIS-NPP data on the 361st day in 2015 with the de-noised MODIS-NPP data obtained from the DWT algorithm with the original MODIS-NPP data for the same date, the predicted NPP estimates showed better quality than the de-noised MODIS-NPP (Figure 6), indicating that the data quality of the de-noised data alone was too low in the presence of heavy clouds and heavy aerosols.

The results of comparing the different NPP products with the NPP in the 20 sample plots are shown in

Table 2. Comparison of NPP data between the predicted MODIS-NPP, de-noised MODIS-NPP, original MODIS-NPP, and measured NPP based on field observations of sample plots (NPP units: Mg ha⁻¹ a⁻¹).

Plot#	Species	Measured NPP	Original NPP	RE (%)	De-Noised NPP	RE (%)	Predicted NPP	RE (%)
1	Pinus massoniana	3.72	2.78	25.27	2.89	22.31	3.49	6.18
2	Bamboo	3.65	3.34	8.60	3.64	0.39	3.39	7.23
3	Eucalyptus	3.92	3.10	20.92	3.24	17.35	3.33	15.05
4	Acacia confusa	3.90	2.80	28.21	3.30	15.39	3.15	19.23
5	Ficus	3.88	2.90	25.31	3.31	14.75	3.27	15.78
6	Pinus massoniana	4.17	3.20	23.26	3.36	19.42	3.39	18.71
7	Mucuna Birdwoodiana	4.62	2.90	32.90	3.2	26.41	4.26	7.79
8	Mulberry	4.58	3.40	25.69	3.40	25.69	3.66	20.01
9	Pinus massoniana	4.64	3.37	27.37	3.45	25.65	3.34	28.02
10	Schima	4.07	3.12	28.10	3.57	12.39	3.18	21.96
11	Cheery blossom	3.70	2.92	21.13	3.68	0.60	3.40	8.16
12	Ficus	4.01	3.30	17.80	3.76	6.35	3.37	16.06
13	Camphor	4.70	3.25	30.85	3.18	32.34	3.76	20.00
14	Pinus massoniana	3.96	3.72	6.06	3.24	18.18	3.57	9.85
15	Camphor	4.04	3.30	18.29	3.30	18.29	3.35	17.05
16	Lychee	4.05	3.05	24.71	3.34	17.55	3.87	4.46
17	Camphor	4.50	3.17	29.56	3.50	22.22	3.85	14.44
18	Delonix regia	4.23	3.20	24.35	3.72	12.06	4.11	2.84
19	Tea tree	4.44	3.30	25.60	3.74	15.68	3.95	10.95
20	Alsophila spinulosa	4.96	3.80	23.38	4.02	18.94	4.19	15.52
Mean		4.19	3.19	23.58	3.45	17.10	3.59	13.96
SD		0.38	0.28		0.26		0.34
rRMSE (%)			24.98		19.50		15.67

. The average NPP of the sample plots was 4.19 Mg ha⁻¹ a⁻¹. The average NPP was 3.59 Mg ha⁻¹ a⁻¹ for the predicted NPP, 3.45 Mg ha⁻¹ a⁻¹ for the de-noised NPP, and 3.19 Mg ha⁻¹ a⁻¹ for the original MODIS-NPP data. The RE relative to plot NPP estimates varied from 2.84% to 28.02% for the predicted MODIS-NPP data and from 6.06% to 37.23% for the original MODIS-NPP data. The rRMSE was 15.67% for the predicted MODIS-NPP data, 19.5% for the de-noised MODIS-NPP data, and 24.98% for the original MODIS-NPP data, indicating that the method combining the DWT and the EKF notably increases the estimation accuracy of NPP in the study area. The increase in accuracy could be attributed to de-noised data from the DWT and effective prediction of the EKF. The results indicate that the predicted MODIS-NPP using the proposed method is more precise than the original MODIS-NPP and the de-noised MODIS-NPP when monitoring above- and below-ground organic accumulated vegetation matter.

5. Discussion

5.1. Comparison with Other Similar Studies

Previous research approaches for removing noise from images to construct high-quality MODIS NDVI and LAI time-series datasets [10,13,14,17,23,38] while providing valuable information have some limitations.

Firstly, these methods do not perform well under conditions of heavy cloud contamination, heavy aerosols, and atmospheric variability, among other factors. Figure 6 shows that the vegetation in the west of study area cannot be found due to heavy cloudiness, while the NPP estimates using the novel method presented in this study can be displayed clearly. To quantify the de-noising performance, we computed the REs between the original and predicted data, and the REs between the de-noised and measured data (Figure 5b). The results showed that big REs occurred during the rainy season (from April to September), which further indicates that de-noised results cannot improve the NPP quality under conditions of heavy cloud contamination, heavy aerosols, and atmospheric variability.

Secondly, these methods aim at constructing high-quality MODIS time-series datasets for existing MODIS data and do not predict the data for the next time period. The method proposed in this paper leverages an effective wavelet de-noising function and a robust EKF predicting technique, where the wavelet transform is used to remove noise in MODIS-NPP images, and the EKF is applied to deal with de-noised images modeled as a triply modulated cosine function to predict MODIS-NPP data for the next time period. Figure 7 shows that the high-quality MODIS-NPP data (the last blue point) was predicted using the novel method on the 361th day in 2015, which implies that the novel method can be used to predict the next eight-day NPP data in the future.

Thus, this novel method can also produce high-quality MODIS data in cases of heavy cloud contamination, heavy aerosol, and atmospheric variability, and can also predict MODIS-NPP values in a time series. Compared to original and de-noised NPP, the predicted values in the study area increased the estimation accuracy of forest NPP, reducing the rRMSE by 9.31% and 3.83%, respectively, which supports the novel method.

5.2. Prospects for Future Studies

In our research, a large amount of NPP data was needed to meet the high accuracy prediction requirements of the EKF algorithm. Some original MODIS-NPP data were invalid due to rain, heavy clouds, and heavy aerosol. Thus, we used the time-distance weighting approach to generate baseline NPP data before predicting the NPP data. However, this approach is not suitable for generating MODIS-NPP data during periods of rapid vegetation change, which may limit the generation of high-precision NPP products. Thus, the approach presented in this study should be used with caution. In addition, future work should include additional test sites to confirm the robustness of our proposed method.

6. Conclusions

In this study, we proposed a novel method combining the DWT with the EKF to accurately predict MODIS-NPP data when the original product is heavily contaminated by factors such as clouds or aerosols. For this purpose, the DWT algorithm of the wavelet transform successfully suppressed white speckles in the original MODIS-NPP data. The mean amplitude phase derived from EKF appeared to predict highly accurate NPP estimates based on the de-noised results. The results showed that compared to the original MODIS-NPP data, the method significantly decreased the mean RE and rRMSE values of the predicted MODIS-NPP data of the study area by about 9.62% and 9.31%, respectively, which indicates that this method has great potential for generating new high-quality MODIS-NPP data under conditions of invalid original MODIS-NPP values.

Even though the proposed method improved the accuracy effectively, some limitations still exist. One deficiency is that this method is not suitable for periods of rapid vegetation change, affecting the prediction accuracy of this period, and therefore the approach should be applied with caution.