Canadian Journal of Remote Sensing Journal canadien de télédétection ISSN: (Print) (Online) Journal homepage: https://www.tandfonline.com/loi/ujrs20 Fusing Landsat and SAR Data for Mapping Tropical Deforestation through Machine Learning Classification and the PVts-β Non-Seasonal Detection Approach Yonatan Tarazona, Alaitz Zabala, Xavier Pons, Antoni Broquetas, Jakub Nowosad & Hamdi A. Zurqani To cite this article: Yonatan Tarazona, Alaitz Zabala, Xavier Pons, Antoni Broquetas, Jakub Nowosad & Hamdi A. Zurqani (2021) Fusing Landsat and SAR Data for Mapping Tropical Deforestation through Machine Learning Classification and the PVts-β Non- Seasonal Detection Approach, Canadian Journal of Remote Sensing, 47:5, 677-696, DOI: 10.1080/07038992.2021.1941823 To link to this article: https://doi.org/10.1080/07038992.2021.1941823 © 2021 The Author(s). Published by Informa View supplementary material UK Limited, trading as Taylor & Francis Group. Published online: 26 Jul 2021. Submit your article to this journal Article views: 2989 View related articles View Crossmark data Citing articles: 3 View citing articles Full Terms & Conditions of access and use can be found at https://www.tandfonline.com/action/journalInformation?journalCode=ujrs20 CANADIAN JOURNAL OF REMOTE SENSING 2021, VOL. 47, NO. 5, 677–696 https://doi.org/10.1080/07038992.2021.1941823 Fusing Landsat and SAR Data for Mapping Tropical Deforestation through Machine Learning Classification and the PVts-b Non-Seasonal Detection Approach Fusion des donne es Landsat et RSO pour la cartographie de la deforestation ^ce a tropicale gra  la classification par apprentissage automatique et a  l’approche de detection non saisonnie re PVts-b Yonatan Tarazonaa,f , Alaitz Zabalab , Xavier Ponsb , Antoni Broquetasa , Jakub Nowosadc , and Hamdi A. Zurqanid,e a Department of Signal Theory and Communications, School of Telecommunication Engineering, Universitat Politecnica de Catalunya, 08034 Barcelona, Spain; bGrumets Research Group, Department Geografia, Edifici B, Universitat Autonoma de Barcelona, Bellaterra, 08193 Catalonia, Spain; cInstitute of Geoecology and Geoinformation, Adam Mickiewicz University, Krygowskiego 10, 61-680, Poznan, Poland; dDepartment of Forestry and Environmental Conservation, Clemson University, 261, Lehotsky Hall Box, Clemson 340317, SC, USA; eUniversity of Arkansas Agricultural Experiment Station, Arkansas Forest Resources Center, University of Arkansas at Monticello, Monticello, AR 71655, USA; fAmerican Program in GIS and Remote Sensing, California 15112, Peru ABSTRACT ARTICLE HISTORY This article focuses on mapping tropical deforestation using time series and machine learn- Received 17 December 2020 ing algorithms. Before detecting changes in the time series, we reduced seasonality using Accepted 5 June 2021 Photosynthetic Vegetation (PV) index fractions obtained from Landsat images. Single and multi-temporal filters were used to reduce speckle noise from Synthetic Aperture Radar (SAR) images (i.e., ALOS PALSAR and Sentinel-1B) before fusing them with optical images through Principal Component Analysis (PCA). We detected only one change in the two PV series using a non-seasonal detection approach, as well as in the fused images through five machine learning algorithms that were calibrated with Cross-Validation (CV) and Monte Carlo Cross-Validation (MCCV). In total, four categories were obtained: forest, cropland, bare soil, and water. We then compared the change map obtained with time series and that obtained with the classification algorithms with the best calibration performance, revealing an overall accuracy of 92.91% and 91.82%, respectively. For statistical comparisons, we used deforestation reference data. Finally, we conclude with some discussions and reflections on the advantages and disadvantages of the detections made with time series and machine learning algorithms, as well as the contribution of SAR images to the classifications, among other aspects. RÉSUMÉ Cet article porte sur la cartographie de la de forestation tropicale a l’aide de series chronolo- giques et d’algorithmes d’apprentissage automatique. Avant de detecter les changements dans les series temporelles, nous avons reduit la saisonnalite en utilisant les fractions de l’in- dice de vegetation photosynthetique (PV) obtenues a partir des images Landsat. Des filtres mono et multitemporels ont ete appliques pour reduire le bruit de chatoiement des images RSO (Radar a synthese d’ouverture) (c’est-a-dire ALOS PALSAR et Sentinel-1B) avant de les fusionner avec les images optiques par le biais de l’analyze en composantes principales (ACP). Nous avons detecte un seul changement a la fois dans les series de PV en utilisant une approche de detection non saisonniere, ainsi que dans les images fusionnees, gr^ace a cinq algorithmes d’apprentissage automatique qui ont ete calibres par validation croisee (CV) et validation croisee de Monte-Carlo (MCCV). Au total, quatre categories ont ete obte- es, sol nu et eau. Ensuite, nous avons compare la carte des change- nues: fore^t, terres cultive ments obtenue avec les series temporelles et celle obtenue avec les algorithmes de classification les plus performants en matiere de calibration, revelant une precision globale CONTACT Yonatan Tarazona geoyons@gmail.com Supplemental data for this article can be accessed at https://doi.org/10.1080/07038992.2021.1941823. ß 2021 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http://creativecommons.org/licenses/by- nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way. 678 TARAZONA ET AL. de 92.91% et 91.82% respectivement. Pour les comparaisons statistiques, nous avons utilise des donnees de reference sur la deforestation. Enfin, nous concluons par des discussions et des reflexions sur les avantages et les inconvenients des detections effectuees avec les series temporelles et les algorithmes d’apprentissage automatique, ainsi que sur la contribu- tion des images RSO dans les classifications. Introduction process of effective and rapid monitoring of deforest- ation events; and (ii) they depend heavily on seasonal- Deforestation resulting from anthropogenic activities, ity in the time series. To address these limitations, and defined as the conversion of forest to another Tarazona et al. (2018) proposed a novel detection land use, as well as the long-term reduction of the approach called PVts-b, which has one calibration tree canopy cover below the minimum 10% threshold parameter and detections that are not dependent on (FAO 2018), have had a significant impact on terres- trial biodiversity and productivity all over the World. the seasonal component of the time series. In fact, forest change is one of the major processes of However, in tropical regions cloud cover seriously global land-cover change (Ding et al. 2019). Forests affects the possibilities of adequately mapping defor- play a crucial role in ecosystem services, including estation (Zhu et al. 2019), and some tropical coun- carbon sequestration, climate and water cycle regula- tries have cloud cover that exceeds the long-term tion, and maintenance of biodiversity (Gibson et al. yearly average frequency of 80%. Persistent cloud 2011; Turubanova et al. 2018). In particular, tropical cover inhibits full optical coverage from Landsat-like regions have been undergoing rapid changes in forest sensors even when compositing is performed over cover mainly due to land-use change. For example, 1–2 years (Souza et al. 2013). Synthetic Aperture tropical deforestation accounts for 10–15% of Radar (SAR) can penetrate clouds, and therefore has anthropogenic greenhouse gas emissions (Houghton the potential of complementing optical-based forest 2013). Tropical deforestation has been occurring since monitoring systems (Joshi et al. 2016). Several the 1980s (van der Werf et al. 2009; Hansen et al. authors have used these images to map deforestation 2013; Achard et al. 2014) and threatens biodiversity, (e.g., Park and Chi 2008; Trisasongko 2010; Jia and ecosystem services, and water security on the planet. Lei Wang 2018). Since Advanced Land Observing Remote sensing has played a key role in studying Satellite (ALOS) Phased Arrayed L-band Synthetic the environment, especially since the launch of Aperture Radar (PALSAR) data have become freely Landsat-1 in 1972 (Vidal-Macua et al. 2017). Free and available and with the launch of the Sentinel-1A and open access to the Landsat archive has dramatically 1B C-band SAR satellites in 2014 and 2016 (Torres benefited operational applications, scientific studies, et al. 2012), for the first time, dense SAR time series and discoveries based on analyses of large numbers of are free and openly available for the tropical region. Landsat images (Zhu et al. 2019). For example, annual Nevertheless, shorter wavelength C-band (5.5 cm) forest change has been globally mapped using avail- SAR is generally less useful for monitoring forest able Landsat observations from 2000 to 2012. This change due to the lower penetration depth and rapid global mapping reports net deforestation of 1.5 mil- saturation of the signal over forests (Woodhouse lion km2 worldwide (Hansen et al. 2013). Several 2005). However, the high temporal observation dens- approaches, which can exploit the full temporal detail ity of the Sentinel-1 C-band compensates for the of available archives, have been proposed to map trop- low sensitivity for detecting deforestation, when ical deforestation from local to global scales (e.g., compared to longer wavelength L-band (e.g., ALOS Souza et al. 2005; Asner et al. 2009; Kennedy et al. PALSAR) SAR observations (Reiche et al. 2018). The 2010; Verbesselt et al. 2012; Hansen et al. 2013; Zhu main negative aspect of SAR images is speckle noise, and Woodcock 2014). However, and according to which impoverishes radiometric resolution and Tarazona et al. (2018), the main limitations of most of affects interpretation and classification results. the detection algorithms based on time-series, such as Reducing speckle noise leads to a decrease in spatial Breaks For Additive Season and Trend (BFAST) resolution, and therefore reduces the potential of (Verbesselt et al. 2012) or Continuous Change these images for coverage analyses. Detection and Classification (CCDC) (Zhu and Multi-sensor data fusion approaches (Zhang 2010) Woodcock 2014), are: (i) they have different and combining SAR and Landsat optical sensors have excessive calibration parameters that hinder the demonstrated a clear increase in forest mapping CANADIAN JOURNAL OF REMOTE SENSING 679 accuracy. This is mainly because fusing the different Study area data enhances visual interpretation and also improves Our study area is located in Madre de Dios, a region the performance of quantitative analyses. in the southeast of Peru (see Figure 1). Madre de Dios Furthermore, in recent years there have been signifi- has a flat topography, with a minimum and maximum cant advances in techniques for fusing SAR with altitude of 170 and 360 masl (meters above sea level). optical data, such as the wavelet-merging technique This moderate slope facilitates the co-registration pro- (Hong and Zhang 2008; Lu et al. 2011; Abdikan cess between optical and SAR images. The study area 2016), Principal Component Analysis (PCA) (Walker is included in path/row 003/069 within the World et al. 2010; Pereira et al. 2013; Abdikan 2016) and Reference System-2 (WRS-2) and covers approxi- intensity-hue-saturation (IHS) (Abdikan 2016). mately 1100  1100 Landsat pixels (i.e., 1090 km2). However, there is still no consensus in the scientific Madre de Dios is a hotspot of great biodiversity community on the best method for integrating SAR (Tarazona and Miyasiro-L opez 2020) and 40% of its and optical data for mapping the effects of deforest- area is protected by a set of Protected Natural Areas ation, especially in tropical forests (Pereira et al. and Native Communities. It also has the largest coal 2013). Approaches that combine optical and SAR time reserves in the world (Baccini et al. 2012). However, series imagery to detect deforestation must overcome deforestation that occurs in the region, based mainly various challenges, including accurate co-registration on logging, agriculture, and gold mining (Figure 1c), and speckle noise, among others. is endangering forest ecosystem services and having In this study, we propose detecting deforestation negative effects on the surrounding population and in a tropical region in south-eastern Peru using only the Peruvian state. In fact, the latest reports of defor- optical data (Landsat 5 Thematic Mapper (TM) and estation, made by the Peruvian state through the Landsat 8 Operational Land Imager (OLI)), and national reports of the Ministry of the Environment combining optical data and SAR data, that is, only of Peru (MINAM), indicate an increase in deforest- two dates. An optical and radar image for each date ation in recent years. was used for the data fusion. For the first case, we used long series of Photosynthetic Vegetation (PV) Data and methods index fractions in order to reduce seasonality in the We follow a series of rigorous processes according to time series. These PV fractions were obtained from the premise set out in the introduction (see Figure 2). the Carnegie Landsat Systems Analysis-Lite Section “Satellite image processing” details the pre- (CLASlite) program (Asner et al. 2009) under an processing of optical images with the CLASlite Automated Monte Carlo Unmixing approach software to obtain the PV fractions and a standard (AutoMCU), which provides quantitative analysis of pre-processing of the SAR images that includes co- the fractional or percentage cover (0–100%) of live registration with optical images. Section vegetation. In the second case, the PCA method was “Deforestation reference data” provides a review of used for data fusion, which does not lead to a loss the reference data used to evaluate the obtained maps. in spatial resolution (i.e., it does not perform In Section “Fusion of SAR and optical images through smoothing) and makes it possible to obtain the con- PCA” we conducted a brief and highly focused review tribution (%) of the variables in the fusion. of the optical and SAR data fusion method used. Therefore, our objectives address: (i) to detect Finally, in Section “Mapping deforestation” we include changes in the Peruvian rainforest between 2009 and several techniques and standard guidelines for the 2018 (i.e., obtain a single change map) using the mapping and validation of the changing cartography. PVts-b approach through PV series, (ii) to detect changes using machine learning algorithms with the Satellite image processing fused images for the same period as in the previous Landsat data and fraction indices case, (iii) to evaluate the change maps obtained with In this study, 50 images from Landsat TM and OLI the PVts-b approach and the fused images, and (iv) pairs between 1990 and 2018 with a Tier 1 processing to discuss the main advantages and disadvantages of level (Young et al. 2017) were selected. Exactly 38 detections based on time series and machine learn- images correspond to the TM sensor (1990–2011) and ing, as well as the benefits of the contributions of 12 to the OLI sensor (2013–2018). Photosynthetic SAR images to classifications and change detections. Vegetation fractions were obtained from 50 reflectance 680 TARAZONA ET AL. Figure 1. Political map (panel (a)) and Madre de Dios study area with a Landsat 8 OLI background image in RGB combination: SWIR2, SWIR1 and GREEN bands (panel (b)). Panel (d) shows the Digital Elevation Model (DEM) of the Shuttle Radar Topographic Mission (SRTM), both panel (b) and panel (d) represent the same study area. Panel (c) shows gold mining through PlanetLabs images (at 5 m); in addition, the resolution of this image makes it possible to see the characteristics of the gold mining in the study region (i.e., presence of water bodies). images with the CLASlite software (Asner et al. 2009). cloud-free mosaics for each year, we used the algo- In short, CLASlite uses the physical model Spectral rithm proposed by Tarazona et al. (2018). This algo- Mixture Analysis (SMA) to obtain the PV fractions. rithm consists of three steps: (i) it chooses the image SMA assumes that the energy received within the field less affected by atmospheric noise, retaining the image of vision of the remote sensor can be considered as with the highest number of clear pixels and its acqui- the sum of the energies received from each dominant sition time, (ii) it calculates the distances (i.e., differ- endmember (specifically a Linear Mixing Model). ences in acquisition time) of the image selected in Therefore, we would expect advantages in the detec- step (i) to the remaining images, and (iii) it takes the tions because we are working with changes in photo- clean pixels of the images, starting with the image synthetic activity at the subpixel level. Given the that has the shortest distance to the image selected in potential of CLASlite subpixel analysis, the AutoMCU step (i) and ending with the image most tempor- approach has been widely used to map forest disturb- ally distant. ance and deforestation worldwide (e.g., Asner et al. A total of 28 yearly PV fraction images was created 2005, 2009; Carlson et al. 2012; Allnutt et al. 2013; (1990–2018 time series) (of a total of 50 images). Bryan et al. 2013; Dlamini 2017; Tarazona et al. 2018). ALOS PALSAR and Sentinel-1 data After obtaining the PV fractions we deal with the We used ALOS PALSAR data in Fine Beam Dual atmospheric artifacts using the Function of mask Polarization (FBD) with Terrain Corrected for this (Fmask) (Zhu and Woodcock 2012) to mask clouds study. The image was acquired on August 12th, 2009, and cloud shadows. Subsequently, and to generate with a spatial resolution of 12.5 m, a swath of 70 km, CANADIAN JOURNAL OF REMOTE SENSING 681 Figure 2. Flowchart of the pre-processing and processing of the optical and SAR images is shown. in HH and HV-polarizations. Both polarizations were speckle-noise just as we did with ALOS PALSAR. For acquired in gamma-naught radar backscattering coef- terrain correction, we used the Range-Doppler Terrain ficients (c0 ) (Rosenqvist et al. 2007). We softened the Correction model and a three arc-second DEM salt and pepper effect by using a Simple Speckle Filter obtained from the SRTM and downloaded automatic- (SSF) through the Lee filter (Lee 1980) with a movable ally by the SNAP software. The spatial output reso- 3  3 window. We then converted it into decibels (dB) lution after applying the correction to the terrain was using the SNAP software (Sentinel Application 10 m. The VV and VH radar backscattering coeffi- Platform, see http://step.esa.int/main/download/) of cients (r0 in dB) were obtained from digital number the European Space Agency. (DN) values through Eq. (1): Sentinel-1 is a C-band Synthetic Aperture Radar dBi, j ¼ 10  log10 ðDNi, j 2 =A2 Þ (1) (5.5 cm wavelength) imaging mission (Torres et al. 2012) that supports a wide variety of applications, where, i, j represents the position (row and column), including deforestation assessment. Sentinel-1B with and A represents the calibration factor, which can be VV þ VH polarization (S1B-VV þ VH) data were looked up in the sigma naught values in the lookup acquired in Interferometric Wide Swath mode (IW table (LUT). with a 250 km2 swath). The S1B-VH þ VV image was In addition, accurate multi-temporal co-registration acquired on December 22nd, 2018, and was down- is an essential prerequisite for change detection and loaded in Ground Range Detected (GRD) format. The was also carried out with the SNAP software. The pre-processing was also performed with the SNAP Optical Landsat (see Section “Fusion of SAR and software and consisted of three steps: (i) calibration, optical images through PCA”) image was selected as (ii) terrain correction, and (iii) smoothing of master and the SAR images as slaves, and a Root 682 TARAZONA ET AL. Table 1. Number of tie points, threshold of accuracy and A total of nine variables (Zi ) are considered in the RMSE in pixel units obtained after co-registration between two periods. Specifically, six optical bands for both optical and radar data. 2009 and 2018, two HH and HV polarizations for Fusion Tie points RMSE mean RMSE std 2009, and another two VV and VH polarizations for 2009 SWIR1 (master) þ HH (slave) 457 0.18167 0.09060 SWIR1 (master) þ HV (slave) 489 0.16986 0.07929 2018. In addition, we include the HH/HV and VV/ 2018 SWIR1 (master) þ VV (slave) 568 0.24612 0.11843 VH ratios for 2009 and 2018, respectively. The two SWIR1 (master) þ VH (slave) 566 0.24455 0.11167 std: standard deviation. ratios are effective in detecting deforestation and are correlated with forest phenology (e.g., Zeng et al. 2014 Mean Square Error (RMSE) accuracy threshold of 0.5 and Frison et al. 2018). pixels was applied (i.e., 15 m). The RMSE of the co- The Landsat images used for the fusion were registration between the optical and SAR images for selected to match the acquisition date of the SAR the two periods was less than 0.5 pixels (Table 1). images as closely as possible. Therefore, August 28, 2009 and September 6, 2018 were selected from among the available optical images. To fuze the Deforestation reference data images we used the “fusionRS” function (see In order to validate the detected deforestation events, Appendix B for more details). This function uses the two sources have been used as reference data (i.e., prcomp function embedded in R (R Core Team 2020) deforested areas and intact forest masks). First, we to execute the PCA. The outputs of the “fusionRS” used the official deforestation reports from the function are the eigenvalues, eigenvectors, the correl- Peruvian government carried out by the Ministry of ation of the variables with the PC and the contribu- the Environment (MINAM 2015) (see: http://geobos- tion in percentage of each variable in each PC in a ques.minam.gob.pe/geobosque/view/index.php for user-friendly and straightforward way (all the docu- more details). The methodology of these reports is mentation on how to execute the function and also based on the work of Hansen et al. (2008) and examples can be found in Appendix B). Potapov et al. (2011). These reports show the deforest- Moreover, given the nature of the SAR images and ation between 2001 and 2017. Second, we used the optical images, all variables were centered and scaled global forest change data from Hansen et al. (2013), (i.e., standardized) (Eq. (2)) before computing the which were obtained from time-series analysis of PCA to prevent any variable from having a dispropor- Landsat images including the OLI sensor. These tionate influence on the PCA. reports are open-access and are usually used as xl ground truth (global forest change between 2000 and Z¼ (2) r 2018). We used the Google Earth Engine platform (Gorelick et al. 2017) to download them (see: https:// Where x is the original variable, l is the mean, r the earthenginepartners.appspot.com/science-2013-global- standard deviation of the original variables and Z is forest/download_v1.6.html for more details). the standardized variable. Fusion of SAR and optical images through PCA Mapping deforestation To fuze optical (Landsat imagery) and SAR images Non-seasonal detection approach (ALOS PALSAR and Sentinel-1) we used the PCA Time-series methods for mapping deforestation are method that is widely known by the scientific com- generally classified into seasonal and non-seasonal munity (e.g., Pereira et al. 2013; Werner et al. 2014; approaches. Seasonal approaches are widely used and Abdikan 2016). The goal is to reduce the dimensional- a variety of algorithms have been proposed (e.g., those ity through linear combinations of the original varia- proposed by Verbesselt et al. 2012 or Zhu and bles, obtaining the contribution of the optical bands Woodcock 2014), but all of them are highly depend- and SAR in the principal components (PC). We ent on the seasonal component of the time series. As believe that a description of the mathematical depth a result of this dependence, they have a number of of the PCA is important for understanding the shortcomings that make it difficult to map changes in internal process of fusion of the two images with dif- the forest quickly and accurately. However, the non- ferent observation geometries and which few studies seasonal approach (e.g., PVts-b proposed by Tarazona have been able to develop; therefore, we present the et al. (2018)), does not depend on seasonality and mathematical description in Appendix A. detections may be more accurate than the others. CANADIAN JOURNAL OF REMOTE SENSING 683 Figure 3. Examples of time series of forests that have been deforested (points in blue). The Landsat images in the top row are in combination SWIR2, NIR and RED. The middle row contains the PV fractions scaled between 0 and 100. Finally, the lower row shows the detections made by the PVts-b approach with b ¼ 5. The PV series were smoothed, before detecting deforestation, with the “smootH” function (see Appendix B). Therefore, we used the PVts-b approach proposed Classification methods to detect deforestation between 2009 and 2018. For After applying PCA to fuze the images, five popular this, we used series of PV fractions to minimize the robust machine learning algorithms in the R language effect of seasonality and improve detection accuracy. were used to detect changes. Support Vector Machine This approach consists of a series of processes that are (SVM) was the first algorithm used. The SVM classi- detailed in the work of Tarazona et al. (2018). The fier is a supervised non-parametric statistical learning optimal threshold values (b) (Eq. 3) are 5 and 6 for technique that does not assume a preliminary distri- PV series. In this work we used the value b ¼ 5. bution of input data (Mountrakis et al. 2011; P^ ossa and Maillard 2018). Its discrimination criterion is a Li, j ¼ li, j  b  ri, j (3) hyperplane that separates the classes in the multidi- mensional space in which the samples that have estab- lished the same classes are located, generally some Where, Li, j is the lower limit, li, j is the mean, ri, j is training areas. The SVM classifier has been efficiently the standard deviation, b is the threshold magnitude used for mapping vegetation and forest (Knorn et al. and i, j is the row and column position of a 2009; Duro et al. 2012). The second algorithm used given pixel. was Random Forest (RF). RF is a derivative of For example, Figure 3 shows two time series of the Decision Tree which provides an improvement over PV fraction for pixels where the forest is deforested. DT to overcome the weaknesses of a single DT (Pal Figure 3a shows a deforestation event due to agricul- 2005; He et al. 2017; Vlachopoulos et al. 2020). The tural activity in 2011, and Figure 3b shows another prediction model of the RF classifier only requires deforested area in 2010, caused by gold mining activ- two parameters to be identified: the number of classi- ity, which has a large negative impact on the fication trees desired, known as “ntree,” and the num- environment. ber of prediction variables, known as “mtry,” used in 684 TARAZONA ET AL. each node to make the tree grow (Rodriguez-Galiano Table 2. Number of points for the calibration of the classifiers et al. 2012). The third algorithm used was Decision and also for obtaining the change cartography for the years 2009 and 2018, as well as samples for validating the changes Tree (DT). DT is also a supervised non-parametric obtained with the classifiers and the PVts-b approach. statistical learning technique, where the input data is Calibration divided recursively into branches (Chen et al. 2018) Year 2009 Year 2018 Change validation depending on certain decision thresholds until the Classifiers data are segmented into homogeneous subgroups. Forest 225 225 Training This technique has substantial advantages for remote Cropland 125 125 336 Bare soil 105 105 Testing sensing classification problems due to its flexibility, Water 105 105 224 intuitive simplicity, and computational efficiency Total 560 560 560 Change maps 550 (Friedl and Brodley 1997). The fourth algorithm used PVts-b approach was Naive Bayes (NB). The NB classifier is an effect- ive and simple method for image classification based (10 iteration) and an iteration equal to 100 for CV on probability theory. The NB classifier assumes an and MCCV respectively. In order to select the best underlying probabilistic model and captures the classifier, which will then serve as input for the 2009 uncertainty about the model in a principled way, that and 2018 classifications, it was necessary to minimize is, by calculating the occurrence probabilities of differ- the testing bias (the mean of errors ent attribute values for different classes in a training ei ¼ 1  overall accuracy) (Eq. (4)) and the testing set (Pradhan et al. 2010). The last algorithm used was variance (the standard deviation of errors ei) (Eq. (5)) Neuronal Networks (NN). This classification consists so that the model works well for both the training of a neural network that is organized into several data and the test data. This process and analysis were layers, that is, an input layer of predictor variables, written in R with the help of some libraries, such as one or more layers of hidden nodes, in which each caret and willcoxCV. node represents an activation function acting on a Xn weighted input of the previous layers’ outputs, and an Bias ¼ e =n (4) i¼1 i output layer (Liu et al. 2019). rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The fine-tuning of model hyperparameters can sig- Xn Variance ¼ i¼1 i ðe  eÞ2 =n (5) nificantly increase the target in both training and test- ing; however, it also increases the cost of training Where i represents the position of a specific iteration models and the risk of overfitting. Therefore, in this research, we decided to use the default parameters and n represents the total number of itera- (see S1 in supplementary material) after examining a tions executed. large number of studies (Belgiu and Draguţ 2016; Accuracy assessment and map comparison Tyralis et al. 2019; Zurqani et al. 2019). A total of 1670 training points were taken to calibrate Since our study area is located in the Amazon, obviously the classes to be identified were mainly for- and validate the different cartographic products gener- est, cropland, bare soil, and water. Classification of ated (Table 2). These points were selected through the crop types or differentiating secondary forests from interpretation of the satellite images themselves, truth- virgin forests will not be addressed in this study due terrain reference data elaborated by MINAM (2015) to the nature and objectives of the research. The and, on some occasions, we used high-resolution objective is to separate the forest from other classes. images such as Google Earth. First, we selected 560 points for the calibration of the five classifiers for Calibration each of the two years, 2009 and 2018, through a Undoubtedly, given a large number of machine learn- stratified random sampling approach (Olofsson et al. ing algorithms, it was necessary to select the one with 2014). These same points were also used to obtain the the best performance in the classification, that is, the definitive 2009 and 2018 class maps using the classi- algorithm in which the training and testing data used fier with the best performance. In both periods, the converge the learning iteratively to a solution that 560 points were split into training (60%) and testing appears to be satisfactory. Therefore, calibration at (40%) sets through proportional stratified random this stage was crucial and was carried out through sampling. Finally, to validate the change map gener- Cross-Validation (CV) and Monte Carlo Cross- ated with the classifiers and the changes with the Validation (MCCV). We used k groups equal to 10 PVts-b approach, we selected 550 points using CANADIAN JOURNAL OF REMOTE SENSING 685 Figure 4. Panel (a) and (b) show the combinations NIR, RED and GREEN and SWIR2, SWIR1, GREEN from the year 2009 (Landsat- 5 TM). Panel (c) and (d) show the same region with polarizations HH and HV (ALOS PALSAR—2009) respectively. Panel (e) and (f) show the combinations NIR, RED and GREEN and SWIR2, SWIR1, GREEN from the year 2018 (Landsat-8 OLI). Panel (g) and (h) show the same region with polarizations VV and VH (Sentinel 1B—2018). Finally, panel (i) and (j) show the fused images of 2009 and 2018 respectively, in both cases using a combination with the first three principal components. stratified random sampling (80% for intact forest and Results the rest for deforested areas). Fusion of optical and SAR images To evaluate the accuracy of both the image fusion classifications and the detections of the PVts-b The fusions of optical and SAR data for the years approach, we used indicators such as overall accuracy, 2009 and 2018 are shown in Figure 4a–h. In both user’s accuracy, producer’s accuracy, and kappa index panels, they are shown in a 1-2-3 combination of 9 (Congalton and Green 2008). These indicators have principal components (Figure 4i–j). Only the first been widely discussed and used in a large number of principal components provided the greatest variance papers, such as Schultz et al. (2016) and Duro et al. in explaining the data in both 2009 and 2018 (Figure (2012), among others. 5). At this stage we used the jambu elbow (Tibshirani 686 TARAZONA ET AL. Figure 5. Proportion of variance of each PC and the cumulative variance for the year 2009 is shown in panel (a); panel (b) shows the same for the year 2018. Table 3. Contribution in percentages (%) of the variables to the principal components after the fusion of optical and SAR images for the years 2009 and 2018. Contributions (%) Year 2009 PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 Proportion of variance (%) 48.601 21.013 10.806 8.067 5.890 4.579 0.722 0.286 0.035 BLUE 20.017 0.641 0.053 7.611 2.763 5.132 56.201 0.000 7.945 GREEN 21.138 0.204 0.061 6.675 2.370 2.780 2.486 1.194 63.212 RED 20.618 1.015 0.369 7.844 1.656 0.183 38.885 1.498 28.332 NIR 0.409 34.911 8.292 7.437 1.304 37.021 1.351 8.363 0.344 SWIR1 9.507 17.939 3.343 8.973 5.077 10.057 0.984 42.521 0.058 SWIR2 17.963 1.647 0.208 0.574 4.488 27.996 0.009 46.417 0.106 HH 1.725 18.939 35.887 15.133 27.235 2.118 0.002 0.001 0.000 HV 5.402 23.393 2.858 5.995 49.887 14.039 0.078 0.000 0.000 HH/HV 3.216 1.305 48.925 39.755 5.216 0.670 0.000 0.002 0.000 Year 2018 Proportion of variance (%) 57.723 21.485 9.650 6.618 3.964 0.339 0.103 0.089 0.029 BLUE 16.698 0.772 1.568 8.781 9.807 53.780 0.200 0.001 8.389 GREEN 17.618 0.375 1.715 5.779 7.148 2.609 0.018 1.728 63.006 RED 17.816 0.125 2.151 2.780 6.491 41.863 0.163 0.547 28.058 NIR 0.382 28.121 17.701 40.529 11.477 0.547 0.117 1.100 0.021 SWIR1 7.373 22.559 4.034 22.293 2.661 0.388 3.596 36.929 0.163 SWIR2 12.684 11.431 0.357 18.788 0.110 0.008 3.052 53.207 0.358 VV 11.513 16.100 5.491 0.220 10.170 0.094 52.850 3.555 0.001 VH 12.173 5.548 9.535 0.500 48.260 0.590 21.554 1.835 0.001 VV/VH 3.737 14.965 57.444 0.325 3.871 0.116 18.446 1.093 0.000 Lower contributions in each PC are shown in red and blue for optical and radar bands respectively. The last PC that will not be used in the classification are shown in gray. et al. 2001) to decide which components to take out On the other hand, Table 3 shows the contribu- of the equation due to the poor contribution to the tions (%) of the original variables to the principal variance. In fact, for 2009 and 2018, only the six and components. In general, for the year 2009 almost all five principal components had an accumulated vari- optical bands had a significant contribution in PC1 ance of 98.96% and 99.44% respectively (Figure 5). In (except NIR band), the NIR and SWIR1 bands in addition, increasing the number of principal compo- PC2:3, all bands in PC4 (except SWIR2), all (except nents would bring us closer to a context of using RED and NIR bands) in PC5 and the NIR, SWIR1 almost all the variables involved and an increase in and SWIR2 bands in PC6. While all polarizations had noise, which could ultimately hinder the classification a strong contribution in all PC1:6 (except HH in PC1, process through automatic learning algorithms. HH/HV in PC2 and PC6, HV in PC3 and HH in CANADIAN JOURNAL OF REMOTE SENSING 687 Figure 6. The calibrations were performed through CV and MCCV, panel (a) is for the year 2009 and panel (b) is for the year 2018. The minimum, maximum and average of the overall accuracies for each algorithm are shown. In addition, we show the vari- ance obtained for the k ¼ 10 groups and for 100 iterations. PC6). In addition, it is important to note that, of all However, all classifiers had an acceptable overall the SAR bands, the HV polarization had the largest accuracy greater than 90%. In both periods, the classi- contribution in PC1, followed by HH/HV. For 2018, fier with the lowest overall average accuracy was DT. the behavior is quite similar. Like 2009, almost all optical bands made a large contribution in PC1 (except NIR band), NIR, SWIR1 and SWIR2 bands in Visual examination of thematic maps PC2, RED, NIR and SWIR1 in PC3, all bands in PC4 Since the objective is to quantify deforestation, it is and all bands (except SWIR1 and SWIR2) in PC5. very important that the algorithm optimally discrimi- While all polarizations, especially VH and VV, had a nates the forest from other types of cover. Figures 7 significant contribution in PC1:5, with the exception and 8 show the thematic maps obtained by the five of PC4. classifiers in 2009 and 2018 respectively. In both peri- ods, especially in 2009, the visual interpretation Algorithm calibrations between the thematic maps obtained by the five classi- From a series of five automatic learning algorithms, fiers revealed an erroneous amount of the cropland the idea was to select the classifier that meets two class represented in the extreme north and southwest conditions: (i) high overall accuracy and (ii) best per- of the study area. In fact, a visual inspection of the formance in predicting the forest and cropland cate- Landsat images and high-resolution PlanetLabs images gories within the validation data. revealed that these areas are predominantly covered The bias and variance of the calibration of the clas- by forests. However, there is also a certain level of sifiers had notable differences with CV and MCCV deforestation in these areas, which could have led to (Figure 6). That is, calibrations with CV had a slightly confusion between forest and cropland classes. greater variance than calibrations with MCCV, mainly In general, all algorithms produced reasonably because CV only explores some of the possible ways accurate land cover maps in both study periods. in which data could have been partitioned (e.g., test- Nevertheless, the results of the RF algorithm in 2009 ing and training). In fact, MCCV allowed 90 times and NN in 2018 were of special interest as they show more exploration of possible partitions from a total of a better performance in separating forest and cropland 280 C560  1: In short, the calibrations indicated that the cover with low omission and commission levels, and algorithms with the best performance in overall accur- with an overall accuracy of 95.09% in both years acy and the lowest variance were SVM, NN, and RF (Table 4). A more detailed analysis showed that the (overall average accuracy of 97.50% and 98.00% with average commission for the forest class in 2009 was CV and MCCV respectively) for 2009 and NN (over- 0.33% higher than in 2018, and the average omission all average accuracy of 93.93% and 94.80% with CV for the same class in 2009 was 2.86% lower than in and MCCV, respectively) for 2018 (Figure 6). 2018. However, the “salt and pepper” effect was 688 TARAZONA ET AL. Figure 7. Classifications using the six principal components of the fused image in the year 2009. Panels (a), (b), (c), (d) and (e) represent the classifications obtained with SVM, DT, RF (the highest overall accuracy), NB and NN, respectively. present in both periods, especially in 2018 for the 92.91%. This was higher by 1.09% than the overall classifiers NB, DT, and NN (Figures 7 and 8). accuracy obtained with the RF-NN classifiers (Table 5). In addition, the lowest omissions and commissions were obtained with the PVts-b approach and RF-NN Change validation and statistical comparisons classifiers with 9.85% and 9.98% respectively. Visual Figure 9 shows deforestation areas obtained with the analysis revealed that the commission of the PVts-b PVts-b approach, the RF-NN classifiers (i.e., changes approach was mainly in gold mining areas, while the obtained using the RF and NN classifiers) and the ref- omission was in agricultural and gold mining areas, erence data. Visually, the detections of the PVts-b which were not optimally detected. Moreover, for a approach (Figure 9a) and the RF-NN algorithms rigorous statistical comparative analysis, a forest/non- (Figure 9b) reveal a lot of similarity with the reference forest mask of the reference data from Hansen et al. data of Hansen et al. (2013) (Figure 9c) and the (2013) was used to eliminate commissions from detec- MINAM reports (Figure 9d), although a closer inspec- tion approaches. The statistics revealed that there was tion showed some visual differences between the two no significant difference between the mapped areas detection methods, mainly south of the study region with respect to the reference data. In fact, both the (i.e., in elliptical white line). That is, a compact and PVts-b approach and the RF-NN algorithms detected continuous quantification should be expected. a total of 18142.74 ha and 16276.05 ha respectively, However, only the change obtained by the RF-NN while the reference data from Hansen et al. (2013) classifiers shows something more detectable and recorded 19708.35 ha between 2009 and 2018. For the logical (i.e., detection of deforestation in areas of gold period 2009–2017, the PVts-b approach detected activity should include not only the bare soil but also 15998.31 ha, while Hansen et al. (2013) and MINAM the water bodies in the areas). reference data reported 16903.89 ha and 17877.15 ha However, the change maps showed that the overall respectively. Finally, from the changes obtained with accuracy obtained with the PVts-b approach was the RF-NN classifiers, it is important to mention that CANADIAN JOURNAL OF REMOTE SENSING 689 Figure 8. Classifications using the five principal components of the fused image in the year 2018. Panels (a), (b), (c), (d) and (e) represent the classifications obtained with SVM, DT, RF, NB, and NN (the highest overall accuracy) respectively. Class “Difference” represents the difference map of the same classifier between 2009 and 2018 with ground truth (Hansen et al. 2013). between 2009 and 2018 the forests mainly changed to possible to map the forest class optimally in 2009 (see cropland 7169.49 ha, bare soil 4507.56 ha, and omission in Table 4). water 4599.00 ha. To explore the performance of image fusion in detections, we detected changes only with optical images. These detections were underestimated by Discussion 2660 ha with respect to the fused images, which is shown in S3 in the supplementary material. We As expected, image fusion showed that the HV and believe that the difference between what is detected VH cross-polarizations had a slightly greater contribu- with fused images and optical images is relatively tion than the HH and VV polarizations in PC1 (Table small but it is an indicator that the fusion of optical 3). In general, for the year 2009, it is observed that and radar data contributes to the mapping of almost all the radar bands made a significant contri- deforestation. bution to PC2:5, while for the year 2018, the highest To analyze whether speckle noise and feature reso- radar contributions were observed in PC1:5 (except in lution had any influence on change detection, we PC4). In addition, the NIR band (a representative decided to apply a Multi-temporal Speckle Filter optical band for vegetation mapping) made a lower (MSF) (Quegan and Yu 2001) to the SAR images contribution than any other optical band in PC1 before fusing with the optical data. For the year 2009, (even lower than the radar bands) in both study peri- we obtained HH, HV, and HH/HV polarizations from ods. However, this was not the case with the RED a total of 7 ALOS PALSAR images; while for the year band. On the other hand, the vegetation structure 2018, we obtained VV, VH, and VV/VH polarizations (Shimada et al. 2014) and the detection of non-photo- from a total of 11 S1B images (see S4 in the synthetic activity (e.g., fallen leaves) due to the pene- Supplementary material). The change map showed tration of ALOS PALSAR images, probably made it deforestation of 16132.77 ha between 2009 and 2018 690 TARAZONA ET AL. Table 4. Overall accuracy and Kappa index for all classification algorithms and both dates. Water Forest Cropland Bare soil Overall Accuracy Kappa Index OM CO OM CO OM CO OM CO 2009 SVM 95.09 93.07 0.00 0.00 0.00 9.00 22.22 2.78 2.44 2.44 DT 93.30 90.55 2.13 0.00 1.10 10.0 26.67 8.33 2.44 4.76 RF 95.09 93.11 2.13 0.00 1.10 6.25 17.78 5.13 2.44 6.98 NB 93.30 90.58 2.13 0.00 2.20 9.18 24.44 10.53 2.44 4.76 NN 95.09 93.09 0.00 0.00 1.09 8.16 20.00 5.26 2.44 2.44 2018 SVM 94.20 91.94 0.00 4.65 2.33 9.68 20.37 4.44 0.00 0.00 DT 92.86 90.12 0.00 10.87 3.49 6.74 20.37 10.42 4.65 0.00 RF 94.64 92.58 0.00 0.00 4.65 7.87 12.96 9.62 2.33 0.00 NB 93.75 91.38 0.00 4.65 6.98 6.98 14.82 11.54 0.00 0.00 NN 95.09 93.18 0.00 0.00 2.33 9.68 16.67 4.25 0.00 0.00 Omission errors (OM) and commission errors (CO) are also shown. All values are in percentages (%). For overall accuracy, the highest values for each date are shown in bold, as well as lower commission and omission errors in the forest and cropland categories. The green background indicates the best classifier for each year, selected to maximize overall accuracy as well as minimize omission and commission errors for forest and cropland areas. See S2 in the supplementary material for more details. Figure 9. Panel (a) and (b) show deforestation detections (2009–2018) obtained from the PVts-b approach and the RF-NN classi- fiers respectively. Panel (c) and (d) show the reference data from Hansen et al. (2013) and MINAM respectively. All the panels have the SWIR1, NIR and RED combination from the year 2018 (Landsat 8 OLI) as the background. Table 5. Overall accuracy, kappa index, and the average com- mission and omission errors (average commissions and omis- (in both periods, we used PC1:6  99% of cumulative sions from deforestation classes and intact forest) are shown variance) using the algorithms that performed best in for the PVts-b approach, and for the change obtained with the calibrations, that is, RF in 2009 as well as 2018. the RF-NN classifiers. Although MSF resulted in sharper and clearer boun- Overall accuracy Kappa index Commission Omission daries (Figure 10a–c in yellow and red lines), SSF PVts-b approach 92.91 78.36 11.67 9.85 together with optical data quantified the deforestation Change with RF-NN 91.82 72.26 9.98 16.7 better. Similar results were obtained in the study by All values are in percentages (%). CANADIAN JOURNAL OF REMOTE SENSING 691 Figure 10. Column (a) shows SAR images without a speckle filter, column (b) images with SSF, column (c) shows a Multi-temporal Speckle Filter (MSF). To evaluate the reduction of speckle noise we used three indicators: (i) visual interpretation, (ii) Signal Level Ratio (SLR) (White et al. 2020), and (iii) the standard deviation of the intensity before and after applying filters for areas with higher signal echo intensity (P3) and areas with lower signal echo intensity (P4) (column d). The parameters for the MSF filter were used by default: (i) filter: filter Lee, (ii) Number of looks: 1, (iii) Window size: 7  7, and (iv) Sigma: 0.9. Table 6. SLR is a ratio of the unfiltered image to the filtered deforestation better, as can be seen in section S5 on image averaged over a region (R1 and R2). the supplementary material. Mean SLR Another aspect that needs to be discussed is the RAW SSF MSF SSF MSF presence of “salt and pepper” in all the classifications R1 obtained, although much less in the period 2009. In HH 0.2139 0.2140 0.2130 0.9995 1.0042 HV 0.0642 0.0642 0.0640 1.0000 1.0031 general, SVM revealed a lower proportion of this type R2 of noise in the two study periods. It is not clear VV 0.1417 0.1414 0.1774 1.0021 0.7988 VH 0.0331 0.0331 0.0390 1.0000 0.8487 whether the speckle noise of the SAR images was the A mean of one indicates that the filtered image preserves the signal level main origin of the “salt and pepper” since Duro et al. accurately. If the mean varies from one, then there is a bias in the filter- (2012), who analyzed optical images, also obtained ing method. similar results to ours with classifications contami- Mirelva and Nagasawa (2017) (although these authors nated by this type of noise in several algorithms of only worked with radar data). Overall, no substantial machine learning, although in a smaller proportion, differences were found between the two filter types by classifiers like SVM and RF. In addition, the study (Figure 10d, panels 3 and 4, Table 6), although for the area is a region with very particular characteristics areas with higher signal echo intensity the SSF filter and it is a challenge to map deforestation since the performed better. In addition, the MSF seems to have areas deforested by gold mining trigger not only the an effect on the contribution of the SAR bands to the presence of bare soil (i.e., sand) but also water bodies first PC, since lower contributions were observed with (Figure 1c). This is a product of the activity itself respect to SAR bands with SSF (see Table 3 and Table which hinders certain pre-processing methods such as S4 for more details). masking of atmospheric noise in optical images Although the classifications obtained with the five (Tarazona et al. 2018). That is, the masking of atmos- classifiers performed well in overall accuracy, they all pheric noise could mask areas that are deforested by had confusions between forest and cropland. This is mining activity. Therefore, the quantification of defor- mainly due to the similarity of the spectral response estation would reveal a certain degree of omission. of the two coverages, predominantly in the optical This is clearly seen in the quantifications (yellow electromagnetic spectrum. In these situations an lines) carried out with the PVts-b approach and the increase in the number of variables, for example, reference data, but not so much with the quantifica- vegetation indices or other images from different tions carried out with the RF-NN classifiers (see dates, improves the separability of certain relatively Figure 9). This difficulty is not present in the SAR confusing covers and slightly improves the overall images (this can be seen in the detection of fusion), accuracy of the classification, but does not detect where water bodies resulting from gold mining are 692 TARAZONA ET AL. recorded with low backscatter values (see white lines deforestation compared to using only optical data. in Figure 9b). The contribution of SAR bands, especially HV, HH/ Another important point of discussion is the detec- HV, VH, and VV, was substantially important in both tion threshold of the PVts-b approach, since the omis- study periods. We believe that the benefits can be sion of the detections occurred mainly in agricultural greater if the SAR images have been acquired after the coverages. To avoid this deficiency, it is possible to optical images since in areas of high precipitation lower the threshold to b ¼ 3; however, this could trig- cloud-free optical images are usually only available ger a slight increase in the commission. In any case, between June and September. Therefore, when only for future work, it will be important to assess to what optical images are used there is a void of three extent a threshold of b ¼ 3 can increase commission months with no record of the changes in the forest. levels, although the latter is possible to combat with a Regarding speckle noise, although MSF better pre- forest/non-forest mask. serves the characteristics of the entities (better “spatial Finally, according to this work’s context, it is resolution”), it seems to decrease the contributions important to make an analysis of the detection (unlike SSF) of the SAR bands in the data fusion and approaches based on time series and those based on underestimates the detection of deforestation. In add- machine learning. Each approach has advantages and ition, there were no significant differences between disadvantages in certain aspects. For example, it is the backscatter values obtained after applying MSF clear that a machine learning approach, in contrast to and SSF filters. However, the calibration of the a time series approach, is unlikely to detect slight machine learning algorithms revealed that all classi- changes over time in photosynthetic activity. fiers had an acceptable overall accuracy. In fact, all However, if the atmospheric noise and seasonality of classifications had an overall accuracy of more the time series are not eliminated or minimized cor- than 90%. rectly, it could make any further analysis in the time The overall accuracy of the PVts-b and RF-NN series difficult. On the other hand, the machine learn- detections was very similar, although the highest com- ing approach has the advantage from a computational mission and lowest omission were obtained with the and execution point of view (see S6 in the supplemen- PVts-b approach. We believe that this approach is a tary material). Furthermore, the main global deforest- simple and intuitive detection method. Moreover, it ation mapping methodologies are based on can be applied to images that come from different classifications such as Hansen et al. (2013), Hansen sensors with dissimilar periodicity because it does not et al. (2016), and Asner et al. 2009 (CLASlite), depend on the seasonality of the data (i.e., correlated although lately, thanks to the GEE platform (Gorelick data). In addition, it is important that the data set is et al. 2017), there have been efforts to implement large enough (>5 data) to avoid bias in the calcula- detections based on time series, such as the tion of the mean and standard deviation. LandTrendr method (Kennedy et al. 2010), which is Finally, given the increased interest and activity in embedded in GEE (https://emapr.github.io/LT-GEE/). simple, effective, and low-cost forest monitoring However, most detection methods based on time ser- methods, PVts-b is clearly a strategic tool. Therefore, ies have: (i) excessive calibration parameters, (ii) are the PVts-b approach, the algorithm for fusing optical based strongly on seasonality, and (iii) are complex to and radar data and other functions is implemented in execute for a standard user. On the other hand, the an open-source package called “ForesToolboxRS” (see PVts-b approach is a method that: (i) is simple and Appendix B for more details). intuitive and which does not model seasonality, (ii) has only one calibration parameter to detect deforest- Acknowledgments ation, and iii) can be easily implemented by any The work of Yonatan Tarazona Coronel has been partially standard user. funded by American Program in GIS and Remote Sensing and National Program of Scholarships and Educational Credit (PRONABEC–Peru) through RJ: N 4276-2018- Conclusions MINEDU/VMGI-PRONABEC-OBE and RJ: N 942-2019- Fusing optical and SAR data has been studied at MINEDU/VMGI-PRONABEC-OBE. Xavier Pons is the recipient of an ICREA Academia Excellence in Research length in the past, but the contribution (%) of each of Grant. The GRUMETS Research Group is partially sup- the original variables within the fusion has not been ported by the Catalan Government under Grant SGR2017- addressed optimally. This research shows the benefits 1690. This work has also been supported by the Spanish of combining optical and SAR data for detecting MCIU Ministry through the NEWFORLAND project CANADIAN JOURNAL OF REMOTE SENSING 693 (RTI2018-099397-B-C21 (MCIU/AEI/ERDF, EU)). In add- Carlson, K.M., Curran, L.M., Ratnasari, D., Pittman, A.M., ition, we would especially like to thank Dr. Xin Miao for Soares-Filho, B.S., Asner, G.P., Trigg, S.N., Gaveau, D.A., his kind and appropriate suggestions and comments, as well Lawrence, D., and Rodrigues, H.O. 2012. “Committed as different anonymous reviewers for their timely comments carbon emissions, deforestation, and community land that improved the content of the manuscript. Finally, we conversion from oil palm plantation expansion in West thank the General Direction of Environmental Territorial Kalimantan, Indonesia.” Proceedings of the National Planning of the Ministry of the Environment of Peru for Academy of Sciences of the United States of America, Vol. providing us with some data used in this work. 109(No. 19): 7559–7564. doi:10.1073/pnas.1200452109. Chen, Y., Niu, Z., Johnston, C.A., and Hu, S. 2018. “A uni- fying approach to classifying wetlands in the Ontonagon Disclosure statement River Basin, Michigan, using multi-temporal Landsat-8 OLI imagery.” Canadian Journal of Remote Sensing, Vol. The authors report that there is no conflict of interest. 44(No. 4): 373–389. doi:10.1080/07038992.2018.1526065. Congalton, R.G., and Green, K. 2008. Assessing the Accuracy ORCID of Remotely Sensed Data, Principles and Practices. Boca Raton, FL: CRC Press. Yonatan Tarazona http://orcid.org/0000-0002-5208-1004 Ding, Y., Zang, R., Lu, X., and Huang, J. 2019. “Functional Alaitz Zabala http://orcid.org/0000-0002-3931-4221 features of tropical montane rain forests along a logging Xavier Pons http://orcid.org/0000-0002-6924-1641 intensity.” Ecological Indicators, Vol. 97: 311–318. doi:10. Antoni Broquetas http://orcid.org/0000-0001-9801-9145 1016/j.ecolind.2018.10.030. Jakub Nowosad http://orcid.org/0000-0002-1057-3721 Dlamini, W.M. 2017. “Mapping forest and woodland loss in Hamdi A. Zurqani http://orcid.org/0000-0002-9250-4534 Swaziland: 1990–2015.” Remote Sensing Applications: Society and Environment, Vol. 5: 45–53. doi:10.1016/j. rsase.2017.01.004. References Duro, D.C., Franklin, S.E., and Dube, M.G. 2012. “A com- parison of pixel-based and object-based image analysis Abdikan, S. 2016. “Exploring image fusion of ALOS/ with selected machine learning algorithms for the classifi- PALSAR data and LANDSAT data to differentiate forest cation of agricultural landscapes using SPOT-5 HRG area.” Geocarto International, Vol. 33(No. 1): 21–37. doi: imagery.” Remote Sensing of Environment, Vol. 118: 10.1080/10106049.2016.1222635. 259–272. doi:10.1016/j.rse.2011.11.020. Achard, F., Beuchle, R., Mayaux, P., Stibig, H.J., Bodart, C., FAO. 2018. Global Forest Resources Assessment 2020. Rome, Brink, A., Carboni, S., et al. 2014. “Determination of Italy: Food and Agricultural Organization. tropical deforestation rates and related carbon losses Friedl, M.A., and Brodley, C.E. 1997. “Decision tree classifi- from 1990 to 2010.” Global Change Biology, Vol. 20(No. cation of land cover form remotely sensed data.” Remote 8): 2540–2554. doi:10.1111/gcb.12605. Sensing of Environment, Vol. 61(No. 3): 399–409. doi:10. Allnutt, T.F., Asner, G.P., Golden, C.D., and Powell, G.V.N. 1016/S0034-4257(97)00049-7. 2013. “Mapping recent deforestation and forest disturb- Frison, P.-L., Fruneau, B., Kmiha, S., Soudani, K., Dufr^ene, ance in northeastern Madagascar.” Tropical Conservation E., Le Toan, T., Koleck, T., Villard, L., Mougin, E., and Science, Vol. 6(No. 1): 1–15. doi:10.1177/ Rudant, J.-P. 2018. “Potential of Sentinel-1 data for mon- 194008291300600101. itoring temperate mixed forest phenology.” Remote Asner, G.P., Knapp, D.E., Broadbent, E.N., Oliveira, P.J.C., Sensing, Vol. 10(No. 12): 2049. doi:10.3390/rs10122049. Keller, M., and Silva, J.N. 2005. “Selective logging in the Gibson, L., Lee, T.M., Koh, L.P., Brook, B.W., Gardner, Brazilian Amazon.” Science, Vol. 5747: 480–482. doi:10. T.A., Barlow, J., Peres, C.A., et al. 2011. “Primary forests 1126/science.1118051. are irreplaceable for sustaining tropical biodiversity.” Asner, G.P., Knapp, D.E., Balaji, A., and Paez-Acosta, G. Nature, Vol. 478(No. 7369): 378–381. doi:10.1038/ 2009. “Automated mapping of tropical deforestation and nature10425. forest degradation: CLASlite.” Journal of Applied Remote Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Sensing, Vol. 3(No. 1): 033543. doi:10.1117/1.3223675. Thau, D., and Moore, R. 2017. “Google Earth Engine: Baccini, A., Goetz, S.J., Walker, W.S., Laporte, N.T., Sun, Planetary-scale geospatial analysis for everyone.” Remote M., Sulla-Menashe, D., Hackler, J., et al. 2012. “Estimated Sensing of Environment, Vol. 202: 18–27. doi:10.1016/j. carbon dioxide emissions from tropical deforestation rse.2017.06.031. improved by carbondensity maps.” Nature Climate Hansen, M.C., Krylov, A., Tyukavina, A., Potapov, P.V., Change, Vol. 2: 1–4. doi:10.1038/nclimate1354. Turubanova, S., Zutta, B., Ifo, S., Margono, B., Stolle, F., Belgiu, M., and Draguţ, L. 2016. “Random forest in remote and Moore, R. 2016. “Humid tropical forest disturbance sensing: A review of applications and future directions.” alerts using Landsat data.” Environmental Research ISPRS Journal of Photogrammetry and Remote Sensing, Letters, Vol. 11(No. 3): 034008. doi:10.1088/1748-9326/ Vol. 114: 24–31. doi:10.1016/j.isprsjprs.2016.01.011. 11/3/034008. Bryan, J.E., Shearman, P.L., Asner, G.P., Knapp, D.E., Aoro, Hansen, M.C., Potapov, P.V., Moore, R., Hancher, M., G., and Lokes, B. 2013. “Extreme differences in forest Turubanova, S.A., Tyukavina, A., Thau, D., et al. 2013. degradation in Borneo: Comparing practices in Sarawak, “High resolution global maps of 21st-century forest cover Sabah, and Brunei.” PLoS One, Vol. 8(No. 7): e69679. change.” Science, Vol. 342(No. 6160): 850–853. doi:10. doi:10.1371/journal.pone.0069679. 1126/science.1244693. 694 TARAZONA ET AL. Hansen, M.C., Stehman, S.V., Potapov, P.V., Loveland, T.R., MINAM. 2015. “Protocolo de clasificaci on de Perdida de Townshend, J.R.G., DeFries, R.S., Pittman, K.W., et al. Cobertura en los Bosques H umedos Amaz onicos entre 2008. “Humid tropical forest clearing from 2000 to 2005 los a~nos 2000–2011.” Resolucion Ministerial Per u, Lima, quantified by using multitemporal and multiresolution http://www.bosques.gob.pe/archivo/files/pdf/protocolo remotely sensed data.” PNAS, Vol. 105: 9439–9444. Mirelva, P.R., and Nagasawa, R. 2017. “Single and multi- He, Y., Lee, E., and Warner, T.A. 2017. “A time series of temporal filtering comparison on synthetic aperture radar annual land use and land cover maps of China from data for agriculture area classification.” Paper presented 1982 to 2013 generated using AVHRR GIMMS NDVI3g at ICISPC 2017: Proceedings of the International data.” Remote Sensing of Environment, Vol. 199: 201–217. Conference on Imaging, Signal Processing and doi:10.1016/j.rse.2017.07.010. Communication, Penang, Penang, July 2017. Hong, G., and Zhang, Y. 2008. “Comparison and improve- Mountrakis, G., Im, J., and Ogole, C. 2011. “Support vector ment of wavelet-based image fusion.” International machines in remote sensing: A review.” ISPRS Journal of Journal of Remote Sensing, Vol. 29(No. 3): 673–691. doi: Photogrammetry and Remote Sensing, Vol. 66(No. 3): 10.1080/01431160701313826. 247–259. doi:10.1016/j.isprsjprs.2010.11.001. Houghton, R.A. 2013. “The emissions of carbon from defor- Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., estation and degradation in the tropics: Past trends and Woodcock, C.E., and Wulder, M.A. 2014. “Good practi- future potentia.” Carbon Management, Vol. 4(No. 5): ces for estimating area and assessing accuracy of land 539–546. doi:10.4155/cmt.13.41. change.” Remote Sensing of Environment, Vol. 148: James, G., Witten, D., Hastie, T., and Tibshirani, R. 2013. 42–57. doi:10.1016/j.rse.2014.02.015. An Introduction to Statistical Learning - With Pal, M. 2005. “Random forest classifier for remote sensing Applications in R. New York, NY: Springer. classification.” International Journal of Remote Sensing, Jia, M., and Lei Wang, L. 2018. “Novel class-relativity non- Vol. 26(No. 1): 217–222. doi:10.1080/0143116041 23312 local means with principal component analysis for multi- 69698. temporal SAR image change detection.” International Park, N.W., and Chi, K.H. 2008. “Integration of multitem- Journal of Remote Sensing, Vol. 39(No. 4): 1068–1091. poral/polarization C-band SAR data sets for land-cover doi:10.1080/01431161.2017.1395966. classification.” International Journal of Remote Sensing, Jolliffe, I.T. 2002. Principal Component Analysis (Second Vol. 29(No. 16): 4667–4688. doi:10.1080/0143116080 Edition ed.). New York, NY: Springer. 1947341. Joshi, N., Baumann, M., Ehammer, A., Fensholt, R., Pereira, L.O., Freitas, C.C., Anna, S.J.S.S., Lu, D., and Grogan, K., Hostert, P., Jepsen, M.R., et al. 2016. “A Moran, E.F. 2013. “Optical and radar data integration for review of the application of optical and radar remote land use and land cover mapping in the Brazilian sensing data fusion to land use mapping and mon- Amazon.” GIScience & Remote Sensing, Vol. 50(No. 3): itoring.” Remote Sensing, Vol. 8(No. 1): 70. doi:10.3390/ 301–321. doi:10.1080/15481603.2013.805589. rs8010070. P^ossa, E.M., and Maillard, P. 2018. “Precise delineation of Kennedy, R.E., Yang, Z., and Cohen, W.B. 2010. “Detecting small water bodies from Sentinel-1 data using support trends in forest disturbance and recovery using yearly vector machine classification.” Canadian Journal of Landsat time series: 1. LandTrendr — temporal Remote Sensing, Vol. 44(No. 3): 179–190. doi:10.1080/ segmentation.” Remote Sensing of Environment, Vol. 07038992.2018.1478723. 114(No. 12): 2897–2910. doi:10.1016/j.rse.2010.07.008. Potapov, P., Turubanova, S., and Hansen, M.C. 2011. Knorn, J., Rabe, A., Radeloff, V.C., Kuemmerle, T., Kozak, “Regional scale boreal forest cover and change mapping J., and Hostert, P. 2009. “Land cover mapping of large using landsat data composites for European Russia.” areas using chain classification of neighboring Landsat Remote Sensing of Environment, Vol. 115(No. 2): satellite images.” Remote Sensing of Environment, Vol. 548–561. doi:10.1016/j.rse.2010.10.001. 113(No. 5): 957–964. doi:10.1016/j.rse.2009.01.010. Pradhan, R., Ghose, M.K., and Jeyaram, A. 2010. “Land Lee, J.S. 1980. “Digital image enhancement and noise filter- cover classification of remotely sensed satellite data using ing by use of local statistics.” IEEE Transactions on Bayesian and hybrid classifier.” International Journal of Pattern Analysis and Machine Intelligence, Vol. 2(No. 2): Computer Applications, Vol. 7(No. 11): 1–4. doi:10.5120/ 165–168. doi:10.1109/TPAMI.1980.4766994. 1295-1783. Liu, X., Wang, C., Meng, Y., Wang, H., Fu, M., and Quegan, S., and Yu, J.J. 2001. “Filtering of multichannel Bourennane, S. 2019. “Target detection based on spectral SAR images.” IEEE Transactions on Geoscience and derivation in HSI shadow region classified by convolu- Remote Sensing, Vol. 39(No. 11): 2373–2379. doi:10.1109/ tional neural networks.” Canadian Journal of Remote 36.964973. Sensing, Vol. 45(No. 6): 782–794. doi:10.1080/07038992. R Core Team. 2020. R: A Language and Environment for 2019.1697221. Statistical Computing. Vienna, Austria: R Foundation for Lu, D., Li, G., Moran, E., Batistella, M., and Freitas, C.C. Statistical Computing. 2011. “Mapping impervious surfaces with the integrated Reiche, J., Hamunyela, E., Verbesselt, J., Hoekman, D., and use of Landsat Thematic Mapper and radar data: A case Herold, M. 2018. “Improving near-real time deforestation study in an urban–rural landscape in the Brazilian monitoring in tropical dry forests by combining dense Amazon.” ISPRS Journal of Photogrammetry and Remote Sentinel-1 time series with Landsat and ALOS-2 Sensing, Vol. 66(No. 6): 798–808. doi:10.1016/j.isprsjprs. PALSAR-2.” Remote Sensing of Environment, Vol. 204: 2011.08.004. 147–161. doi:10.1016/j.rse.2017.10.034. CANADIAN JOURNAL OF REMOTE SENSING 695 Rodriguez-Galiano, V.F., Ghimire, B., Rogan, J., Chica- Tyralis, H., Papacharalampous, G., and Langousis, A. 2019. Olmo, M., and Rigol-Sanchez, J.P. 2012. “An assessment A brief review of random forests for water scientists and of the effectiveness of a random forest classifier for land- practitioners and their recent history in water resources. cover classification.” ISPRS Journal of Photogrammetry Water, Vol. 11(No. 5): 910. and Remote Sensing, Vol. 67: 93–104. doi:10.1016/j. van der Werf, G.R., Morton, D.C., DeFries, R.S., Olivier, isprsjprs.2011.11.002. J.G.J., Kasibhatla, P.S., Jackson, R.B., Collatz, G.J., and Rosenqvist, A., Shimada, M., Ito, N., and Watanabe, M. Randerson, J.T. 2009. “CO2 emissions from forest loss.” 2007. “ALOS PALSAR: A pathfinder mission for global- Nature Geoscience, Vol. 2(No. 11): 737–738. doi:10.1038/ scale monitoring of the environment.” IEEE Transactions ngeo671. on Geoscience and Remote Sensing, Vol. 45(No. 11): Verbesselt, J., Zeileis, A., and Herold, M. 2012. “Near real- 3307–3316. doi:10.1109/TGRS.2007.901027. time disturbance detection using satellite image time Schultz, M., Clevers, J.G.P.W., Carter, S., Verbesselt, J., series.” Remote Sensing of Environment, Vol. 123: 98–108. Avitabile, V., Quang, H.V., and Herold, M. 2016. doi:10.1016/j.rse.2012.02.022. “Performance of vegetation indices from Landsat time Vidal-Macua, J.J., Ninyerola, M., Zabala, A., Domingo- series in deforestation monitoring.” International Journal Marimon, C., and Pons, X. 2017. “Factors affecting forest of Applied Earth Observation and Geoinformation, Vol. dynamics in the Iberian Peninsula from 1987 to 2012. 52: 318–327. doi:10.1016/j.jag.2016.06.020. The role of topography and drought.” Forest Ecology and Shimada, M., Itoh, T., Motooka, T., Watanabe, M., Management, Vol. 406: 290–306. doi:10.1016/j.foreco. Shiraishi, T., Thapa, R., and Lucas, R. 2014. “New global 2017.10.011. forest/non-forest maps from ALOS PALSAR data Vlachopoulos, O., Leblon, B., Wang, J., Haddadi, A., (2007–2010).” Remote Sensing of Environment, Vol. 155: LaRocque, A., and Patterson, G. 2020. “Delineation of 13–31. doi:10.1016/j.rse.2014.04.014. bare soil field areas from unmanned aircraft system Souza, C., Roberts, D.A., and Cochrane, M.A. 2005. imagery with the mean shift unsupervised clustering and “Combining spectral and spatial information to map can- the random forest supervised classification.” Canadian opy damage from selective logging and forest fires.” Journal of Remote Sensing, Vol. 45: 489–500. doi:10.1080/ Remote Sensing of Environment, Vol. 98(No. 2–3): 07038992.2019.1635877. Walker, W.S., Stickler, C.M., Kellndorfer, J.M., Kirsch, 329–343. doi:10.1016/j.rse.2005.07.013. K.M., and Nepstad, D.C. 2010. “Large area classification Souza, C.M., Siqueira, J.V., Sales, M.H., Fonseca, A.V., and mapping of forest and land cover in the Brazilian Ribeiro, J.G., Numata, I., Cochrane, M.A., Barber, C.P., Amazon: A comparative analysis of ALOS/PALSAR and Roberts, D.A., and Barlow, J. 2013. “Ten-year Landsat Landsat data sources.” IEEE Journal of Selected Topics in classification of deforestation and forest degradation in Applied Earth Observations and Remote Sensing, Vol. the Brazilian Amazon.” Remote Sensing, Vol. 5(No. 11): 3(No. 4): 594–604. doi:10.1109/JSTARS.2010.2076398. 5493–5513. doi:10.3390/rs5115493. Werner, A., Storie, C.D., and Storie, J. 2014. “Evaluating Tarazona, Y., Mantas, V.M., and Pereira, A.J.S.C. 2018. SAR-optical image fusions for urban LULC classification “Improving tropical deforestation detection through using in Vancouver Canada.” Canadian Journal of Remote photosynthetic vegetation time series – (PVts-b).” Sensing, Vol. 40(No. 4): 278–290. doi:10.1080/07038992. Ecological Indicators, Vol. 94: 367–379. doi:10.1016/j.eco- 2014.976700. lind.2018.07.012. White, L., Brisco, B., Murnaghan, K., Pasher, J., and Duffe, Tarazona, Y., and Miyasiro-L opez, M. 2020. “Monitoring J. 2020. “Temporal filters for mapping phragmites with tropical forest degradation using remote sensing. C-HH SAR data.” Canadian Journal of Remote Sensing, Challenges and opportunities in the Madre de Dios Vol. 46(No. 3): 376–385. doi:10.1080/07038992.2020. region, Peru.” Remote Sensing Applications: Society and 1799770. Environment, Vol. 19: 100337. doi:10.1016/j.rsase.2020. Woodhouse, I.H. 2005. Introduction to Microwave Remote 100337. Sensing. Boca Raton, FL: CRC Press. Tibshirani, R., Walther, G., and Hastie, T. 2001. “Estimating Young, N.E., Anderson, R.S., Chignell, S.M., Vorster, A.G., the number of clusters in a data set via the gap statistic.” Lawrence, R., and Evangelista, P.H. 2017. “A survival Journal of the Royal Statistical Society: Series B, Vol. guide to Landsat preprocessing.” Ecology, Vol. 98(No. 4): 63(No. 2): 411–423. doi:10.1111/1467-9868.00293. 920–932. doi:10.1002/ecy.1730. Torres, R., Snoeij, P., Geudtner, D., Bibby, D., Davidson, Zeng, T., Dong, X., Quegan, S., Hu, C., and Uryu, Y. 2014. M., Attema, E., Potin, P., et al. 2012. “GMES Sentinel-1 “Regional tropical deforestation detection using ALOS mission.” Remote Sensing of Environment, Vol. 120: 9–24. PALSAR 50 m mosaics in Riau province, Indonesia.” doi:10.1016/j.rse.2011.05.028. Electronics Letters, Vol. 50(No. 7): 547–549. doi:10.1049/ Trisasongko, B.H. 2010. “The use of polarimetric SAR data el.2013.4254. for forest disturbance monitoring.” Sensing and Imaging: Zhang, J. 2010. “Multi-source remote sensing data fusion: An International Journal, Vol. 11(No. 1): 1–13. doi:10. Status and trends.” International Journal of Image and 1007/s11220-010-0048-8. Data Fusion, Vol. 1(No. 1): 5–24. doi:10.1080/1947983 Turubanova, S., Potapov, P.V., Tyukavina, A., and Hansen, 0903561035. M.C. 2018. “Ongoing primary forest loss in Brazil, Zhu, Z., and Woodcock, C.E. 2012. “Object-based cloud Democratic Republic of the Congo, and Indonesia.” and cloud shadow detection in Landsat imagery.” Remote Environmental Research Letters, Vol. 13(No. 7): 074028. Sensing of Environment, Vol. 118: 83–94. doi:10.1016/j. doi:10.1088/1748-9326/aacd1c. rse.2011.10.028. 696 TARAZONA ET AL. Zhu, Z., and Woodcock, C.E. 2014. “Continuous change combination) we will apply the Lagrange multiplier: detection and classification of land cover using all avail-   able Landsat data.” Remote Sensing of Environment, Vol. Lðu1 Þ ¼ ut1 Ru1  k ut1 u1  1 (4) 144: 152–171. doi:10.1016/j.rse.2014.01.011. oL and then ou1(derive with respect to the incognita u1 ) and Zhu, Z., Wulder, M.A., Roy, D.P., Woodcock, C.E., Hansen, equal to zero: M.C., Radeloff, V.C., Healey, S.P., et al. 2019. “Benefits of the free and open Landsat data policy.” Remote Sensing of oL ¼ 2Ru1  2R ¼ 0 ) ðR  kI Þu1 ¼ 0 (5) Environment, Vol. 224: 382–385. doi:10.1016/j.rse.2019. ou1 02.016. Eq. (5) can be deduced that jR  kI j ¼ 0 and with roots Zurqani, H.A., Post, C.J., Mikhailova, E.A., and Allen, J.S. k1 , k2 , . . . , kn , ordered from highest to lowest k1 > k2 > 2019. “Mapping urbanization trends in a forested land- . . . > kn : On the other hand, if we multiply to the left by scape using Google Earth Engine.” Remote Sensing in ut1 , we conclude that: Earth Systems Sciences, Vol. 2: 173–182. doi:10.1007/ s41976-019-00020-y. ut1 ðR  kI Þu1 ¼ 0 ) ut1 Ru1 ¼ k ) varðY1 Þ ¼ ut1 Ru1 ¼k Appendix A To maximize the variance of Y1 , we have to take the largest eigenvalues k1 and the corresponding eigenvectors All the variables (optical and SAR) are centered and scaled, and it is a question of calculating a new set of bands (Y1 , u1 , so that the first component is defined as Y1 ¼ Zu1 : Y2 , … , Yp ) uncorrelated among themselves, whose varian- Therefore, all components are expressed as the product of ces progressively decrease. Each component Yj the matrix of eigenvectors multiplied by the vector Z con- (j ¼ 1, . . . , p) is a linear combination of the original taining the original variables: variables: Y ¼ uZ Y1 ¼ u11 Z1i þ u12 Z2i þ    þ u1p Zpi (1) Being, In abbreviated notation we will say Y1 ¼ Zu1 . So the variance of the first component will be: Y ¼ ðY1 Y2 .. Yp Þ u 1 Xn 1 ¼ ðu11 u12    u1p u21 u22    u2p             up1 up2    upp ÞZ varðY1 Þ ¼ Y 2 ¼ Y1t Y1 i¼1 1i (2) n n . Eq. (1) in Eq. (2) ¼ ðZ1 Z2 .. Zp Þand fvarðY1 Þ ¼ k1 varðY2 Þ   . 1 t t 1 u1 Z Zu1 ¼ ut1 Zt Z u1 ¼ ut1 Ru1 (3) ¼ k2 .. varðYp Þ ¼ kp n n The expression n1 Zt Z is the matrix of correlations (R) on With the same logic it is possible to find the eigenvalues standardized variables (i.e., centered and scaled) (Jolliffe and eigenvectors for the remaining components. 2002; James et al. 2013). The first component Y1 is obtained so that its variance is maximum, that is, it is about finding Pp Y1 maximizing Appendix B. Code varðY1 Þ ¼ ut1 Ru1 with the restriction j¼1 ut1 u1 ¼ 1: Since the objective is to maximize a function of several The source code and full instructions on how to execute the variables subject to a constraint and the incognita being u1 functions can be found through GitHub (please see https:// (the unknown vector that gives the optimal linear github.com/ytarazona/ForesToolboxRS for more details).