Predicting Pinus monophylla forest cover in the Baja California Desert by remote sensing

Study area

Sierra La Asamblea is located in Baja California’s central desert (−114°9′W 29°19′N, elevation range 280–1,662 m, Fig. 1). The climate in the area is arid, with maximum temperatures of 40 °C in the summer (García, 1998). The sierra is steeper on the western slopes, with an average incline of 35°, and with numerous canyons with occasional springs and oases. Valleys and plateaus are common in the proximity of the Gulf of California. Granite rocks occur south of the sierra and meta-sedimentary rocks along the north and southeast of the slopes. The predominant type of vegetation is xerophilous scrub, which is distributed at elevations ranging from 200 to 1,000 m. Chaparral begins at an altitude of 800 m, and representative specimens of Adenostoma fasciculatum, Ambrosia ambrosioides, Dalea bicolor orcuttiana, Quercus tuberculata, Juniperus california, and P. monophylla are also present at elevations above 1,000 m. Populations of the endemic palm tree Brahea armata also occur in the lower parts of the canyons with superficial water flow and through the rocky granite slopes (Bullock & Heath, 2006).

Datasets

Pixel-based classification

Classification method

Pixel-based classification was carried out in order to predict four different types of land cover in the study area (P. monophylla, scrub, chaparral, and no apparent vegetation). A supervised classification approach with a backpropagation-type artificial neural network (BPNN) (Tan & Smeins, 1996) was applied. BPNN is widely used because of its structural simplicity and robustness in modeling non-linear relationships. In this study, the BPNN comprises a set of three layers (raster): an input layer, a hidden layer, and an output layer (Richards, 1999). Each layer consists of a series of parallel processing elements (neurons or nodes). Each node in a layer is linked to all nodes in the next layer (Guo et al., 2013).

The first step in BPNN supervised classification is to enter the input layer, which in this study corresponded to the values of the pixels of 10 Sentinel-2 bands and of the NDVI image. Weights were then assigned to the BPNN to produce analytical predictions from the input values. These data were contrasted with the category to which each training pixel belongs, corresponding to Georeferenced sites (Datum WGS-84, 11N) obtained in the field in October 2014 and October 2015.

A stratified random sampling method (Olofsson et al., 2013) was used to generate the reference data in QGIS software (QGIS Development Team, 2016). A total of 2,143 random points were sampled, with at least 400 points for each class (Goodchild, 1994). The following classes were considered: (i) P. monophylla, 536 sites, (ii) scrub, 764 sites, (iii) chaparral, 405 sites, and (iv) no apparent vegetation, 438 sites (INEGI, 2013). Class discrimination processes occurred in the hidden layer and the synapses between the layers were estimated by an activation function. We used a logistic function and training rate of 0.20, previously applied to land cover classification (Hepner et al., 1990; Richards, 1999; Braspenning & Thuijsman, 1995). Learning occurs by adjusting the weights in the node to minimize the difference between the output node activation, and BPNN then calculates the error at each iteration with root mean square (RMS) error. The output layer comprised four neurons representing the four target classes of land cover (P. monophylla, Scrub, Chaparral, and no apparent vegetation). Average spectral signatures for the four different types of land cover are shown in Fig. 2.

Relationship between presence of P. monophylla and environmental variables

To model and test the association between presence/absence of P. monophylla in the study area and topographical or climate variables, a Kruskal–Wallis test was used to estimate the difference in the median values in relation to presence and absence of P. monophylla. All variables for which no significant difference between the median values was predicted after Bonferroni correction (α = 0.0005) were excluded from further analysis. The collinearity between the variables with a significant difference between the medians of presence and absence was estimated using the Spearman correlation coefficient (r_s). When the r_s value for the difference between two variables was greater than 0.7, only the variable with the lowest p value in the Kruskal–Wallis test was used in the models (as reported by Salas, Valdez & Michel, 2017 and Shirk et al., 2018). Finally, stepwise multiple linear binominal logistical regression and Random Forest classification including cross valuation (10-fold) were used to model the associations between presence/absence of P. monophylla and the most important topographical and climate variables (Shirk et al., 2018).

Regression and classification including cross-validations were carried out using the trainControl, train, glm (family = “binomial”), and rf functions, as well as the “randomForest” and “caret” packages (Venables & Ripley, 2002) in R (version 3.3.2) (R Development Core Team, 2017). The goodness-of-fit of the logistical regression model was evaluated using the Akaike information criterion, root mean square error, and residual deviance. Validation of the randomForest model was performed using under the curve (AUC; Fawcett, 2006), true skill statistic (TSS; Allouche, Tsoar & Kadmon, 2006), Kappa (Abraira, 2001), specificity, and sensitivity.

Pixel-based classification

We estimated the area of P. monophylla cover with a margin of error (at approximate 95% confidence interval) of 6,653 ± 319 (standard error) ha in Sierra de la Asamblea, Baja California, Mexico (Fig. 3). The confusion matrix of the accuracy assessment is listed in Table 3 including user’s and producer’s accuracy for each class. The supervised classification with BPNN yielded predictions with an overall accuracy of identification of 87.74% (Table 4). This level of accuracy was estimated in the 32 interactions with 0.04 RMS training. The proportion of omission errors in the P. monophylla class was only 2.62%, i.e., 97.39% of the pixels were correctly classified. The shrub class had the larger proportion of omission errors (18.98%). The value of NDVI in the P. monophylla forest fluctuated between 0.30 and 0.41, and in chaparral between 0.24 and 0.28. The smallest values of NDVI corresponded to scrub vegetation, with values between 0.10 and 0.15.

Relationship between presence of P. monophylla and environmental variables

The Kruskal–Wallis test indicated that the median values for ruggedness TRI (p < 2.1 × 10⁻¹⁶), slope (p < 2.2 × 10⁻¹⁶), ruggedness VRM (p = 4.9 × 10⁻⁹), mean temperature in the warmest month (MTWM) (p = 0.000014), Maximum temperature in the warmest month (MMAX) (p = 0.000048), and SPRP (p = 0.00037) were most variable between sites in which P. monophylla was present and absent. The variable slope was closely correlated with ruggedness as well as with MMAX and MTWM (r_s > 0.7). The p_slope of the Kruskal–Wallis test was larger than p_ruggedness and p_MMAX was larger than p_MTWM. Slope and MMAX were therefore excluded from the model analysis. The stepwise multiple linear binominal logistical and Random Forest models showed that the “presence of P. monophylla” model included the independent variables ruggedness, ruggedness VRM, and average temperature in the warmest month (MTWM) (Table 5).

The ruggedness factor was the most influential predictor variable and indicated that the probability of P. monophylla occurrence was larger than 50% when the degree of ruggedness TRI was greater than 17.5 m (Table 5). The ruggedness VRM also indicated that a minimum change in roughness increases the probability of presence of the pine. The probability of occurrence of P. monophylla decreased when MTWM increased from 23.5 to 25.2 °C (Table 5). After cross validation (10-fold), the Random Forest model revealed that the variables ruggedness TRI, ruggedness VRM, and MTWM yielded a greater correlation for their ability to predict presence of the P. monophylla (AUC = 0.920, TSS = 0.690, Kappa = 0.691). The sensitivity was 0.812 and specificity was 0.878.

Pixel-based classification

Predicting the presence of pine forest by using BPNN proved feasible. The NDVI was one of the variables that contributed to the prediction and clearly separated forest cover (NDVI > 0.35) from the other types of vegetation cover (NDVI < 0.20). The overall accuracy of classification (K = 0.87) was similar to that reported in other studies using Sentinel-2A MSI images. For example, Immitzer, Vuolo & Atzberger (2016) reported a K of 0.85 for tree prediction in Europe by using five classes and a random forest classifier. Vieira et al. (2003) reported a K = 0.77 in eastern Amazon using seven classes and 1999 Landsat 7 ETM imagery. However, Sothe et al. (2017) reported K values of 0.98 and 0.90 for, respectively, three successional forest stages and field in a subtropical forest in Southern Brazil by using Sentinel-2 and Landsat-8 data associated with the support vector machine algorithm. Kun et al. (2014) estimated K values of 0.70 to 0.85 for land-use type prediction (including forest) in China by using the support vector machine algorithm classifier and Landsat-8 images of rougher spatial resolution than Sentinel images. The very greater accuracy of predictions by Kun et al. (2014) was probably due to the large-scale of the study and the clearly differentiated types of land considered.

Relationship between presence of P. monophylla and environmental variables

Ruggedness of the terrain was the most important topographic variable, significantly explaining the presence of pines in Sierra La Asamblea (Table 5). Ruggedness, which is strongly positively correlated with slope, may reduce solar radiation, air temperature, and evapotranspiration due to increased shading (Tsujino et al., 2006; Bullock et al., 2008; Allen et al., 2010). The ruggedness indicated by the TRI index explains the presence of the pines because Sierra La Asamblea is heterogeneous in terms of elevation. The VRM index was less important partly because the index is strongly dependent on the vector aspect (Gisbert & Martí, 2010) and in the case of Sierra Asamblea the aspect is very homogeneous and the index values therefore tend to be very low (Table 5), as also reported by Wu et al. (2018). The pines were expected to colonize north facing slopes, which are exposed to less solar radiation than slopes facing other directions. However, the topographical variable aspect was not important in determining the presence of P. monophylla var. californiarum in the study site, possibly because of physiological adaptations regarding water-use efficiency and photosynthetic nitrogen-use efficiency (DeLucia & Schlesinger, 1991), as reported for the P. monophylla, P. halepensis, P. edulis, and P. remota in arid zones (Lanner & Van Devender, 2000; Helman et al., 2017). The Mediterranean climate, with wet winters and dry summers, is another characteristic factor in this mountain range. In the winter in this part of the northern hemisphere, the sun (which is in a lower position and usually affects the southern aspect by radiation) is masked by clouds, rainfall and occasional snowfall (León-Portilla, 1988). During the summer, the solar radiation is more intense, but similar in all directions because the sun is closest to its highest point (Stage & Salas, 2007).

The above-mentioned finding contrasts with those of other studies reporting that north-eastern facing slopes in the northern hemisphere receive less direct solar radiation, thus providing more favorable microclimatic conditions (air temperature, soil temperature, soil moisture) for forest development, permanence, and productivity than southwest-facing sites (Åström et al., 2007; Stage & Salas, 2007; Huang et al., 2009; Marston, 2010; Klein et al., 2014). DeLucia & Schlesinger (1991) reported that P. monophylla populations in the Great Basin California desert with summer rainfall (monsoon) preferred an east-southeast aspect with less intense solar radiation and evapotranspiration.

The probability of occurrence of P. monophylla was also related to the climatic variable MTWM. In Sierra La Asamblea, this pine species was found in a narrow range of MTWM of between 23.5° and 25.2° (Table 1), which, however, is a smaller range than reported for the other pine species (Tapias et al., 2004; Roberts & Ezcurra, 2012). Therefore, this species should adapt well to greater temperatures in the summer (Lanner & Van Devender, 2000), which is usually a very dry period in the study site (León-Portilla, 1988). However, the probability of occurrence was greatest for an MTWM of 23.5 °C (Table 5), which occurred at the top of Sierra La Asamblea, at an elevation of about 1,660 m). We therefore conclude that this species can also grow well when the MTWM is below 23.5 °C. On the other hand, considering MTWM as factor yielded a probability of occurrence of 25%–80%. The spatial resolution of the climatic data by the national database run by the University of Idaho is probably not adequate for describing the microhabitat of P. monophylla (Rehfeldt et al., 2006; Marston, 2010).

Identification of the P. monophylla stands in Sierra La Asamblea as the most southern populations represents an opportunity for research on climatic tolerance and community responses to climatic variation and change.