Multidecadal evaluation of changes in coffee-growing areas using Landsat data in Central Highlands, Vietnam

Abstract Vietnam is the world’s second-largest coffee exporter. Information on coffee-growing areas is thus important for production estimation. This research developed an approach for coffee mapping and change detection using Landsat and DEM data in Central Highlands, Vietnam. The data were processed for 1995, 2005, and 2020 using random forests (RF). The mapping results, verified with reference data, showed overall accuracy and Kappa coefficient higher than 86.9% and 0.74, respectively. These findings were reaffirmed by a close relationship between the satellite-based coffee area and official statistics, with the root mean squared percentage error (RMSPE) and mean absolute percentage error (MAPE) greater than 9.9% and 7.6%, respectively. From 1995 to 2020, the newly planted coffee area increased by 135,520 ha, due to the conversion of forests to coffee plantations. Such distributions of coffee areas and decadal changes in farming activities are critically useful for policymakers to formulate successful crop planning strategies.


Introduction
Coffee is one of the world's most traded agricultural commodities and popular beverages, consumed by at least one-third of the world's population.It is among the most valuable tropical export crops, cultivated on an estimated area of roughly 10.9 million hectares across more than 80 countries (Statista 2019), with approximately 125 million people relying on it for their economic survival (Krishnan 2017).Vietnam is globally the second largest coffee exporter after Brazil (Adriana 2019; ICO 2020; Santos et al. 2021), annually producing roughly 1763,500 tones (GSO 2020).Coffee is the country's second most valuable agricultural export after rice, playing considerable socioeconomic importance in the country because it provides income sources and rural livelihoods for at least two million people living in the Central Highlands.Dak Lak is the home of robusta coffee, and the recurrent neural networks (RNN) due to its less prone to overfitting problems (Caruana and Niculescu-Mizil 2006;Caruana et al. 2008;Miller et al. 2019;Martinez-Castillo et al. 2020;Ali 2021;Tan et al. 2021).It also produces competitive results with time-consuming deep-learning neural networks in different mapping applications using remotely sensed data (Mboga et al. 2017;Wright and Ziegler 2017;Shiferaw et al. 2019;Ge et al. 2020;Sheykhmousa et al. 2020;Vel asquez et al. 2020;Mohajane et al. 2021;Vel asquez et al. 2021).Other advantages of RF include easily implementing and computing acceleration in a parallel structure for a large geospatial dataset, handling numerous and complex input variables, being vigorous to noise and outliers, mitigating data variance without increasing prediction bias, and creating more generalized solutions (Belgiu and Dragut 2016;Miller et al. 2019;Jain et al. 2020;Lawson et al. 2021;Lu 2021).
In this work, rather than comparing classification accuracies among methods, our ultimate goal was to map coffee-growing areas and accordingly delineate changes in coffee farming activities in the Central Highlands of Vietnam during periods 1995À2005 and 2005À2020 using the RF algorithm.The findings obtained from this work in the form of quantitative information on spatiotemporal distributions and decadal changes in coffeegrowing areas could be beneficial to land-use planners and policymakers to formulate successful crop management and agricultural planning strategies in the province.

Study area and data collection
The study region (Dak Lak province), which was dissolved from Dak Lak and Dak Nong provinces in 2004, is situated in the Central Highlands of Vietnam, encompassing roughly 13,030.5 km 2 and lying between 12 40'N 108 3'E (Figure 1).This province was home to approximately 2.1 million people in 2019 with a density of 160 people/km 2 , in which more than 70% of the population was living in rural areas and relying on agriculture for their livelihoods.The region is blessed with vast natural resources, mainly well-drained basaltic soils that are suitable for the cultivation of industrial trees, especially coffee, cashew, rubber, and pepper.The climate is tropical with two distinct seasons.The wet season lasts from May to November, and the dry season is from December to April.The annual precipitation ranges between 1,500 to 2,400 mm, with more than 80% of the total precipitation occurring during the rainy season.The air temperature is between 20-25 C, with the highest temperatures of 27-31 C in April and May.Dak Lak is the largest coffee-producing province in Vietnam, playing a crucial important role in the country's economy of the Central Highlands.Coffee production in the province is typically practiced from October of the previous year to September, and harvested by the end of October to January (Table 1).
A suite of Landsat-5 and 8 imageries (collection 2 level-2 science products), whereas a single tile (path 124, and row 051) covers the entire study region, were collected for 1995, 2005, and 2020 and used in this work.The Landsat-5 data include seven spectral bands, with a spatial resolution of 30 m for bands 1À5 and 7, except band 6 (thermal infrared) acquired at 120 m resolution but resampled to 30 m.The Landsat-8 data comprises nine spectral bands 1À7 and 9 (spatial resolution of 30 m), one panchromatic band 8 (15 m resolution), and two thermal infrared bands 10 and 11 (100 m resolution, but resampled to 30 m).In this research, six spectral bands 1-5 and 7 (Landsat 5), and bands 2-7 (Landsat 8) from the surface reflectance climate data record product for each dataset were used to derive vegetation indices, while the thermal band 6 (Landsat 5) and band 10 (Landsat 8) were used for derivation of the land surface temperature (LST).Landsat imageries have been radiometrically and geometrically corrected to account for atmospheric effects and geometric errors (Feng et al. 2013;Roy et al. 2014;Claverie et al. 2015).The 30-m resolution Aster DEM data were also collected for the derivation of slope and aspect in the study region.

Satellite data preprocessing and derivation of predictive variables
The data of Landsat images (surface reflectance level-2 product) stored as digital numbers were firstly multiplied with a scaling factor of 0.0000275 and subtracted by 0.2 to convert it to the surface reflectance (scale from 0 to 1).The geometric corrections were accordingly performed for the 1995 and 2005 Landsat images using the 2020 Landsat image as a reference image using ArcGIS 10.3 package.The correction process was carried out using at least 30 control points uniformly chosen from district features throughout the study region.The results showed that the root mean square error (RMSE) value was smaller than half of the size of a Landsat pixel.The LST data were retrieved from Landsat's thermal bands 6 (Landsat 5) and 10 (Landsat 8) using the single-channel algorithm (Jim enez-Muñoz and Sobrino 2003; Sobrino et al. 2004).Because the Landsat data in the region were commonly polluted by clouds, we employed a two-step procedure to create yearly composite images for delineating coffee-cultivated areas in the region.First, the pixels contaminated by dense and cirrus clouds were masked out using the cloud mask product.Second, the temporal pixel-based image composite method using the median value was applied to generate yearly cloud-free composite images to reduce the effects caused by cloud contamination, atmospheric attenuation, and surface directional reflectance (Holben 1986;Huete et al. 2002;Guerschman et al. 2009;Roy et al. 2010;Flood 2013;Mountford et al. 2017).
For each year, a total of 43 variables, including 39 Landsat-based vegetation indices, one Landsat-derived LST, and three topographic indices (i.e.DEM, slope, and aspect), were used in this research for mapping coffee-planted areas (Table 2).The rationale for choosing these parameters was that the DEM, slope, and aspect are topographical parameters that characterize changes in topography and water collection, thus highly related to the land suitability for coffee production.The LST and vegetation indices (e.g.NDVI, EVI, SAVI, and NDWI) have been widely exploited for crop monitoring and yield estimation because they are strongly associated with vegetation health and water stress conditions (Huete 1997;Turner et al. 1999;Gao et al. 2000;Hancock and Dougherty 2007;Kanke et al. 2016;Nolasco et al. 2021).Some of these vegetation and topographic indices have been documented to improve the mapping accuracy in various applications, including crop, forest, and soil classification (Fahsi et al. 2000;Xue and Su 2017;Maleki et al. 2020;Chu et al. 2021;Al-Doski et al. 2022;Kordi and Yousefi 2022;Zeng et al. 2022).All these indices (i.e. 30 m resolution) were stacked into one composite scene (43 bands) and used as independent variables for RF regression analysis to estimate coffee-growing areas.

Ensemble machine learning algorithm
The RF regression method was selected in this research to investigate coffee-growing areas, based on its primary advantages over other machine-learning algorithms (e.g.linear regression, logistic regression, ANN, SVM, and RNN), including non-parametric nature, high degree of prediction accuracy, capability to determine the variable importance (i.e. the contribution of each feature variable to the model prediction), and ability to handle a large volume of datasets within a relatively short computing time (Ver Hoef et al. 2001;Liaw and Wiener 2002;Prasad et al. 2006;Ishwaran 2007;Roy and Larocque 2012;Belgiu and Dragut 2016;Mi et al. 2017).This algorithm, using the ensemble learning method for regression (Breiman 2001;Tibshirani et al. 2002), combines predictions from multiple machine-learning algorithms to make a more accurate prediction than a single model.It is known as the most commonly used method for solving complex prediction tasks.A set of decision trees is constructed from training samples using the bootstrapping technique to perform the regression, leaving the remaining samples (from the training dataset) for model testing.By intensively learning training samples, the algorithm can discover specific trends and patterns in the dataset.Once the optimal trends and patterns are determined, it is possible to predict accurate future instances (Breiman et al. 1984;Breiman 2001).To reduce the computation time and evaluate a variable's relevance for prediction, the random-forest predictor importance analysis based on a Gini impurity measure (Loh 2002, Loh andShih 1999;Louppe et al. 2013) was first applied using the training dataset to extract a minimal set of discriminative predictor variables or determine predictor variables important for the model.The greater importance estimates indicate more important predictors.The relationship between the dependent variable (i.e.coffee-planted areas) and independent variables (i.e.Landsat satellite and DEM-derived indices) can be expressed using the following equation: where, Y is the dependent variable, determined by independent variables X, b are unknowns, and e is the error term.
In this research, a total of 20,000 training pixels (i.e.10,000 samples for each class of coffee-planted areas and non-coffee areas) was randomly chosen from the government's inventory coffee-growing map, synchronized with those from the composite scene of independent variables.Of these training samples, 70% and 30% of the data were respectively used for model training and testing.Because the accuracy of RF relies on the optimal number of trees, to ensure the number of trees is sufficiently large to stabilize the mean squared error (MSE) rate, 1,000 trees were used in the RF model based on trial and error results.In addition, because the computation time increases linearly with the number of trees and variables, to efficiently reduce the computational cost of data processing, we performed the estimation of predictor importance using the permutation algorithm to determine the influence of each dependent variable in the model for predicting the response (Breiman et al. 1984;Loh 2002b;Ishwaran 2007).Only variables that greatly influence predictions were used for the modeling.The 10-fold cross-validation was also processed to evaluate the model performance by dividing the data into 10 folds of training and testing sets to train and test the model, respectively.The iteration is repeated until each fold of such 10 folds has been used as the testing dataset.The mean squared error (MSE) is used as an indicator for evaluating the model performance.

Threshold determination and accuracy evaluation
The prediction results of coffee-planted areas (obtained from RF regression) were an image scene, showing the pixel values ranging between 0 and 1, indicating the proportion of coffee-growing area in each pixel.A hardening process was thus applied to convert a mixed pixel to a pure pixel with respect to the desired coffee class or non-coffee class.The threshold was determined based on the intersection between the producer's and user's accuracies calculated from the error matrix.To avoid the inclination of selecting the same pixels for model training and mapping accuracy assessment, we excluded training pixels from the selection of pixels used for the accuracy assessment.1,000 pixels (500 pixels for each class) were randomly extracted from the ground reference map to calculate the producer's and user's accuracies with an incremental interval of 0.05 between 0.1 and 0.9.The intersection point between the two curves of the producer's and user's accuracies was utilized as a threshold value to classify coffee-planted areas.For assessing the mapping results, we utilized the overall accuracy and Kappa coefficient, which were widely used in the field of remote sensing and map comparison for accuracy assessment (Congalton 1991;Pontius and Millones 2011;Congalton and Green 2019).

Importance of predictive variables
The results of predictor importance analysis for the 1995, 2005, and 2022 datasets showed that 10 out of 43 variables both had variable importance values greater than 0.7 (Figure 3).These variables, including the atmospherically resistant vegetation index (ARVI), green atmospherically resistant vegetation index (GARI), green leaf index (GLI), global vegetation moisture index (GVMI), modified normalized difference water index (MNDWI), normalized difference infrared index (NDII), specific leaf area vegetation index (SLAVI), DEM, SLOPE, and LST, were thus chosen for RF regression analysis.Specifically, the ARVI is developed to correct atmospheric influences in regions with a high content of atmospheric aerosol and air pollution (Tanre et al. 1992).The GARI is more sensitive to the chlorophyll concentration than NDVI and ARVI and is useful to measure the rate of photosynthesis and monitor plant stress (Gitelson et al. 1996).The GLI originally designed to measure wheat cover characterizes soil or non-living features and green leaves or stems (Louhaichi et al. 2001).The GVMI is suggested for the retrieval of vegetation water content (Ceccato et al. 2002).The MNDWI is designed for the detection of water bodies while diminishing built-up areas (Xu 2006).The NDII is sensitive to changes in the water content of plant canopies, which is suitable for crop agricultural management, forest canopy monitoring, and vegetation stress identification (Hardisky et al. 1983).The SLAVI is proposed to estimate specific leaf area (Lymburner et al. 2000) because it is strongly correlated with the volume and height of vegetation, net photosynthesis, net primary production, and leaf area index (Broge and Leblanc 2001;Yu et al. 2010).The DEM characterizes topographic variations, as being the most distinctive factor affecting the land suitability for coffee plantation allocation.Likewise, the slope (derived from DEM) measures the steepness of a feature, and also characterises slope suitability levels for coffee cultivation, which is identified as a yield limiting factor, especially in undrained or eroded elevated areas (FAO 1996).The LST represents the inherent physical processes of energy and water exchange, and has a strong negative association with NDVI as it affects crop phenology and yields (Jackson et al. 1981;Sandholt et al. 2002;Gusso 2018).

Model performance and classification accuracies
The results of 10-fold cross-validation showed that the MSE values achieved for the 1995, 2005, and 2020 datasets were 0.087, 0.093, and 0.0097, respectively.The larger MSE values were observed for 1995, mainly attributed to mixed-pixel problems caused by the expansion of coffee-planted areas and coffee agroforestry systems into forested areas.The spectral confusion between coffee and other industrial crops could consequently increase MSE values between the estimates and the actual observations.The model coefficients for each year were used to estimate the probability of coffee-planted areas for 1995, 2005, and 2020, respectively.The results for each year were an image scene, where the pixel value is a continuous number between 0 to 1, indicating the proportion of coffee area within a pixel.The thresholds of 0.82, 0.78, and 0.81, which were obtained from the intersection of the producer's and user's accuracies with a searching range from 0.1 to 0.9, were applied to the 1995, 2005, and 2020 prediction maps to convert mixed pixels to pure pixels concerning two desired classes of coffee or non-coffee.
The accuracies of mapping results were verified in two ways.Firstly, the classification maps were assessed using 1,000 random pixels randomly extracted from the ground reference map to compare with those pixels synchronized from the classification maps.The comparison results showed that the overall accuracy and Kappa coefficient achieved for 1995 were respectively 88.5% and 0.77, while the values for 2005 and 2020 were 89.1% and 0.78, and 86.9% and 0.74, respectively (Table 2).Of 1,000 reference pixels used to check the per-class accuracy of classification maps, the producer's and user's accuracies obtained for the 1995 dataset were 88.6% and 88.4%, while the values for the 2005 and 2020 datasets were 89.0% and 89.2%, and 87.0% and 86.8%, respectively.In general, the mapping results were slightly lower for 2020 because this year's coffee acreage was more spatially dispersed throughout the province due to its expansion into other cropland areas and forests.
Secondly, to check for consistency, the government's statistics on coffee-planted areas were compared, district by district, to the satellite-derived coffee areas for 1995, 2005, and 2020.The overall results showed a slight overestimation, in all cases.The relative errors in area, achieved from the area comparisons between these two datasets, were 7.8% for 1995, 1.1% for 2005, and 6.4% for 2020, respectively.Overall, at the district level, the two datasets had a high degree of agreement (Figure 4).The correlation coefficient of determination (R 2 ) values of the fitted models generally explained more than 98% of the variance for 1995, 2005, and 2020 datasets (P-value < 0.001).The RMSPE and MAPE values for 1995 were 10.2% and 8.4%, while values for 2005 and 2020 were respectively 10.7% and 7.9%, and 9.9% and 7.6%, indicating a good consistency between the satellitederived coffee area and the government's statistics.The cause of some outliers in scatterplots indicating greater or decreased discrepancies was mainly attributed to mixed-pixel issues, especially in districts, such as M' Drak, Lak, and Krong Bong, where coffee was less cultivated, and coffee fields were relatively small and highly fragmented.
Overall, this research's methodology is concentrated on large-scale and crop-specific mapping of coffee plantations as a complicated target using the moderate-resolution Landsat archives to generate multitemporal thematic information on coffee-growing areas over a heterogeneous landscape.In contrast to most earlier satellite-based studies, which focused on typical food crops, such as wheat, corn, and rice, our research burrows into a large-scale classification and decadal change detection of a perennial coffee crop, primarily investigated by a few studies using Landsat and Sentinel images in Central America and Southeast Asian countries (Cordero-Sancho and Sader 2007;Ortega-Huerta et al. 2012;Kawakubo and Machado 2016;Tridawati et al. 2020;Maskell et al. 2021;Le et al. 2022).With the recent release of higher spatiotemporal resolution Sentinel-1 Synthetic Aperture Radar (SAR) and Sentinel-2 optical data since 2014, we see the integration of these two datasets for more accurate and scalable mapping of coffee-growing areas and monitoring short-term changes in cultivation areas and plant phenology over a large area, especially in cloudy tropical regions, has become feasible (Jia et al. 2014;Escobar-L opez et al. 2022).However, this is problematic because the unavailability of long-term historical Sentinel-1 and 2 archives places limitations on an analysis of multi-decadal changes in agricultural practices.
To the best of our knowledge, the methods given in this work represent an initial attempt at regional delimitation of coffee-producing areas and decennial changes in agricultural practices through the utilization of a large number of indexes derived from DEM and cloud-free annual Landsat composite datasets with the aid of machine-learning techniques.The merit of image composition is that it chooses the optimal approximation for each pixel in spectral space, thus minimizing the impact of outlier pixels while retaining the spectral relationship between bands.Our research findings also demonstrate the effectiveness of the RF algorithm for classification of coffee plantations from the satellite composite image using only the 10 most important factors (i.e.vegetation, LST, and topographic indices).These results have coincided with previous studies that demonstrate the value of using remotely-sensed vegetation and surface temperature variables in reliefconditioned locations for improved classification of specific crops, including coffee (Langford and Bell 1997;Sesnie et al. 2008;Schmitt-Harsh 2013;Schmitt-Harsh et al. 2013;Chemura et al. 2017a, Cordero-Sancho andSader 2007;Hunt et al. 2020;Miranda et al. 2020).
Compared to recent research employing moderate remote sensing data (e.g.Landsat and SPOT) that have mistakenly identified coffee-planted areas as natural forests, orchards, and perennial trees (Cayuela et al. 2006;Stibig et al. 2007;Ortega-Huerta et al. 2012;Kawakubo and P erez Machado 2016;Kelley et al. 2018;Hunt et al. 2020), our methods yield reasonably good results, but also reveal a challenge for accurately identifying small patches of coffee plantations, mostly due to topographic effects, low spatial resolution of satellite images, and complicated coffee-farming systems in which coffee fields are often intercropped with other perennials (e.g.nut trees, tea, cacao, mulberries, and peppers).This diversification in terms of coffee height strata and crown closure influences the spectral coffee signature in a variety of ways, subsequently misclassifying coffee fields as other woodland trees due to similar crowns.On the other hand, there are coffee parcels that have not been properly updated with ground reference data can also contribute to the increased mapping error as indicated by lower producer's and user's accuracies in 2020.Another issue with our approach is that it cannot accurately depict complex coffeegrowing landscapes using a composite image since coffee plants' spectral signatures alter with phenological cycle and age.Although attempts have been made to separate the coffee from other crops and perennials in the study region using the time-series vegetation data (derived from SPOT-VGT), the results reveal consistently unsuccessful isolating of coffee plantations as a result of temporal confusion among cropping patterns (Stibig et al. 2007;Nguyen et al. 2020).Given these difficulties in classification of coffee plantations in the study region, our efforts are still being made to improve the model accuracy by exploiting multiple biophysical variables acquired from multi-date imageries of different sensors.

Geographic distribution and decadal changes of coffee planted areas
The classification results showed the spatial distributions of areas planted with coffee in 1995, 2005, and 2020 (Figure 5), which are spatially comparable to those on the government's reference map.In general, coffee has been grown in districts that have an altitude above 300 m with mild weather conditions suitable for growing coffee.Areas cultivated for coffee were relatively dispersed in 1995 (Figure 5a), but much more intensively in 2005 and 2020 (Figures 5b, c).Mapping results for these years were compared to examine decadal changes in areas planted with coffee between these periods (i.e.1995À2005, 2005À2020, and 1995À2020).The results showed that the coffee area remained broadly unchanged over the period 1995À2020, amounting to approximately 90,749 ha (or 6.9% of the total area of 13,315,199.3 ha), whilst the values for the periods 1995À2005 and 2005À2020 were 6.2% and 10.3%, respectively (Figure 6).Between 1995 and 2020, the total area under coffee grew by approximately 135,520 ha (or 10.3% of the province's total area).Specifically, newly cultivated coffee areas in the 1995À2005 and 2005À2020 periods were 7% (or 92,508 ha) and 6.7% (or 88,164 ha) respectively, which were mostly observed in districts of the province's central territories as a result of the conversion of forests to coffee plantations.In contrast, a very small proportion of the coffee acreage (2.1%) declined over these 26 years (1995À2020), and approximately 2.8% and 2.9% for 1995-2005 and 2005-2020, respectively.This decline is primarily attributable to the fact that old coffee trees are grown in steep areas (>15%), where soils are unsuitable and irrigation is inadequate for coffee production, such that the coffee had been cleared for other purposes.These declining areas were dispersed across the study area, particularly in the southern and eastern parts of the study area.On the whole, the main driving force of such an expansion of coffee acreage has been attributed to the conversion of forests into coffee plantations due to the suitability of basalt soils for large-scale coffee production and the high profitability of coffee grains for export in the global market (Koninck 1999;D'haeze et al. 2005;Meyfroidt et al. 2013;Dang and Kawasaki 2017;Kissinger 2020).The region's total coffee production had nearly doubled over the course of 26 years from 107,735 ha in 1995 to 209,955 ha in 2020, in which a higher proportion of the increase in coffee acreage was observed for the period 1995À2005 (1.6%), compared to that of 1.2% for the 2005À2020 period (DLSO 2020).In addition, the expansion of coffee-producing areas may also be caused by the quickening rate of population growth resulting from uncontrolled immigration from other provinces, particularly the central and northern parts of the country, boosted by high coffee prices on the world's coffee market over the past two decades (DLSO 2020; GSO 2020).

Conclusions
The findings of this study confirmed the potential application of our mapping approach to spatially delineate coffee-growing areas, and evaluate decadal changes in coffee-farming practices from 1995 to 2020 using Landsat and DEM-derived indices.The RF-based mapping results, derived from image classification using 10 variables selected based on the analysis of predictor importance, compared with the ground reference data confirmed strong agreement between the two datasets, with the overall accuracy and Kappa coefficient greater than 88.5% and 0.74, respectively.Such findings have also been reaffirmed by the good agreement between the satellite-based coffee area and the government's coffee area statistics at the district level.From 2000 to 2020, the coffee area increased by approximately 135,520 ha, mostly in the study region's central districts.This expansion was mainly due to the benefits of the price of coffee on the market and to the demographic growth accelerated by the immigration of persons from other regions.This research eventually led to the recognition of the merit of the use of Landsat-based LST and vegetation indices and topographic variables for coffee mapping using the RF machine-learning algorithm in Vietnam's Central Highlands.Such a mapping approach could provide quantitative information on spatial distributions of coffee-producing areas and multidecadal changes in farming practices, which was helpful to agronomists and policymakers in addressing crop management and agricultural planning problems in the region.

Disclosure statement
No potential conflict of interest was reported by all authors.

Figure 1 .
Figure 1.Map of the study region concerning elevation levels, administrative boundaries of districts, and national geography of Vietnam.

Figure 2 .
Figure 2. Map of coffee-cultivated areas in the study region overlaid with layers of the road network and water bodies.Ground reference pixels (dots) were used for the accuracy assessment of the coffee mapping results.

Figure 4 .
Figure 4. Relationship between coffee cultivated areas, obtained from the classification of Landsat and DEM-derived indices, and government's coffee area statistics: (a) 1995, (b) 2005, and (c) 2020.

Figure 6 .
Figure 6.Changed cafe areas between (a) 1995 and 2005, (b) 2005 and 2020, and (c) 1995 and 2020.The percentage refers to the total area of the study region (1,315,199.3ha).

Table 2 .
Accuracy assessment results for the coffee classification maps in 2005 and 2017.