2.2.1.

Normalized difference vegetation index

The normalized difference vegetation index (NDVI) is an index derived from reflectance values of the red and near-infrared (NIR) regions of the electromagnetic spectrum [Eq. (1)] and is sensitive to various biophysical vegetation characteristics, such as biomass and percent cover^23–25

NDVI was calculated for each of the 13 Landsat 8 scenes and used as inputs to the classification. For Landsat 8, the bands used in NDVI are b4 (red) and b5 (NIR). NDVI values range from $-1$ to 1; nonland surfaces (such as water and snow) typically assume negative values and land surfaces typically assume positive values. As landscapes become more densely vegetated, the NDVI trends to 1.^11,15,26

2.2.2.

Creating layer stacks for different elevation zones

Vegetation communities are generally constrained by elevation gradients and the inclusion of elevation data can greatly increase the accuracy of mapping efforts.²⁷ A 10-m raster DEM from The National Map²⁸ was downloaded, mosaiced, and reprojected (UTM, zone 12, WGS84) using bilinear interpolation, converted into units of feet, degraded in resolution to 30 m, and finally subset to the study area coordinates. The newer 10-m data were used as it is an updated version of the 30 m that corrects for some of the issues with the original 30-m data, such as edge matching/seam issues, while maintaining the potential of being updated from other high-resolution sources (such as LIDAR). A mask was created from the DEM data for elevations above and below 2591 m. Snow contamination affecting the pixels of several image dates in areas above 2591 m was the main reason for splitting the scene using this elevation. Masking scenes using this threshold allowed more dates to be utilized in areas below this elevation because snow contamination in upland areas would have caused some otherwise excellent scenes to be rejected.

Initially, two layer stacks were produced, one for each Landsat path/row. The layer stack for P36/R37 consisted of bands 1 to 7 and NDVI for each of the seven dates followed by the DEM image to create a 57 band image. The layer stack for P35/R37 consisted of bands 1 to 7 and NDVI for each of the six dates, followed by the DEM image to create a 49 band image. The P35R37 layer stack was further split into two separate image stacks, one masked to elevations above 2591 m, and one below that elevation. Only P35R37 contained elevations above this threshold in the study area.

For areas above 2591-m elevation, only four of the six dates were uncontaminated by snow cover; these four scenes and their derived NDVIs and DEM imagery were used to create a new 33 band layer stack for analysis in the high elevation zone. The Landsat path/row boundaries cut diagonally across the study area which resulted in strips of “extraneous” pixels as the individual image bands did not exactly cover the same area. To address this, all the individual image bands in two of the layer stacks affected by this issue (P35R37 below 2591 m and P36R37) were made into a binary mask with 0 representing no data and 1 representing usable image data. The minimum overlap area (intersection of polygons) was then determined and used to mask out the areas of the layer stack that contained any “no data” values. In summary, three layer stacks were made as follows:

2.2.3.

Unsupervised classification

Unsupervised classifications do not require training data like supervised classification and therefore may be more economical in developing land cover maps for large areas. We explored a unique method using iterations of unsupervised classifications and thematic class matching in an attempt to increase the accuracy of the classification.

For portions of the map above 2591 m, a four-date layer stack (with corresponding NDVIs and DEM) was classified to initially create 10 unsupervised classes. Each of those 10 classes were then subset from the same layer stack and further classified into 10 additional classes, resulting in a total of 100 classes. One second-generation class was further subdivided into five classes, so 104 classes in total were created. Below 2591 m (P35R37), a six-date layer stack (with corresponding NDVIs and DEM) was classified to create 20 unsupervised classes. Each of those 20 classes were further subdivided into 20 additional classes for a total of 400 classes. For P36R37, a seven-date layer stack (with corresponding NDVIs and DEM) was iteratively classified ($20\times 20$) to ultimately create 400 unsupervised classes as well. The portion of the mapping area above 2591 m (only in P35R37) was far less than the portion of the image below this elevation, and needed far fewer classes to map.

The ERDAS Imagine (ver. 2014), ISODATA algorithm was used for the unsupervised classifications, and these multiple tiered classifications ultimately resulted in over 900 classes. For each ISODATA run, the parameters were set as follows: initialize from statistics; # of classes from: xx to xx (i.e., from 20 to 20); minimum size (%): 0.001; maximum SD: 3.00; minimum distance: 4.00; max. merges: 1; initialize means along: principle axis; scaling range: automatic; maximum iterations: 100; convergence threshold: 0.999; skip factors: $X$: 1 $Y$: 1; and classify zeros = not checked.

2.2.4.

Vegetation class assignment

To determine the most likely vegetation class of each of the 904 resultant unsupervised classes, they were compared with the original SWReGAP map. Each class was spatially subset from the SWReGAP map, and the top three SWReGAP vegetation class designations were determined for each unsupervised class. If the “top” SWReGAP class, determined by the number of pixels representing that class, was exceeded by summing the second and third place most likely SWReGAP classes, and these classes were to be grouped into a single class (i.e., Dry-Mesic Spruce-Fir, Mesic Spruce-Fir, etc.), the grouped class designation (i.e., Spruce-Fir, etc.) was now considered the most likely vegetation class. In other words, if the second and third place classes were to be grouped into a single class, and the number of pixels representing these summed classes exceeded the original first place selection, it was determined to be the most likely class. See Fig. 2 for further clarification.

Fig. 2

Flowchart of process to determine “most likely” vegetation/land cover class.

After assigning the most likely class to all the unsupervised classes for P35R37, the Ponderosa and Mixed Conifer classes, that straddled the “border” across the 2591-m elevation threshold, were subset out, combined and reclassified (into 30 unsupervised classes) using the same four-date layer stack used for high elevations. Vegetation classes were then assigned using the technique described previously. This resulted in a much “cleaner transition” across this artificial elevation boundary. After the class assignments were made, the below and above 2591-m portions of P35R37 were initially mosaiced together. This was followed by mosaicing (using the overlay function) the classified area that straddled the elevation threshold “over” the original P35R37 mosaic. The resultant P35R37 mosaic was then mosaiced (using the overlay function) “over” the classified P36R37 section of the mapping area. This resulted in roughly the eastern two-thirds of the mapping area being derived from P35R37 imagery, and the remaining western portion being generated from P36R37 imagery. Finally, a $3\times 3\text{pixel}$ majority filter was applied to eliminate single pixels and small pixel groups that are potential “noise” in the image. Urban area boundary shapefiles from 2016 were then downloaded from the census website,²⁹ screen digitized, and added to the map to create an urban class. See Fig. 3 for flowchart of process discussed above.

Fig. 3

Flowchart of process steps in creating the unsupervised map.

For use in accuracy assessment, the SWReGAP vegetation/land cover map was collapsed and recoded from its original 46 classes found in the study area to the final 16 classes used in the unsupervised classification map using the cross walk shown in Table 1. Classes that did not directly translate from the SWReGAP map to a more generalized class of the unsupervised classification map (i.e., recently burned, etc.) were labeled as “unclassified” on this map. None of the “unclassified” group of vegetation/land cover classes came up as a “likely” class in the analysis.

2.2.5.

Supervised classification

A supervised classification was developed to compare with the unsupervised classification. An initial supervised classification was performed using all the field data [specifically vegetation/land cover class and associated global positioning system (GPS) position] in the study area from the 2001 to 2004 SWReGAP and the 2017 field campaign as training sites (1028 points/pixels in all). No attempt was made to verify that the SWReGAP points were indeed still the same vegetation class, nor was burned vegetation removed or reassigned from the dataset. After the SWReGAP vegetation type field data were generalized into one of the 16 final classes, they were combined with the field data points from the 2017 field campaign. The field data points were then grouped into the individual vegetation class shapefiles using ArcMap (v. 10.5). They were then further subdivided by Landsat path/row. These groups of points that represented the same vegetation class (in a specific path/row) were then opened into ERDAS Imagine (v. 2014) and converted into AOI points, before being combined with the respective Landsat 8, NDVI, and DEM image layer files to create representative signature files. The six-date Landsat 8 image stack (with NDVI and DEM) discussed earlier was used to create the signature files for P35R37, whereas the seven-date Landsat 8 image stack (with NDVI and DEM) was used for P36R37. During the supervised classification, “minimum distance” was used as the parametric rule, whereas a nonparametric rule was not employed. Water, agriculture, Riparian Forest/Woodland, mine, and urban were added back in from the unsupervised classification map (non $3\times 3$ spatial filtered version). Then, a $3\times 3$ median filter was run on the image.

2.3.1.

2017 field campaign

The 2017 field campaign was based on a map of semirandom, accessible points in the study area that included a number of points for each of the generalized vegetation/land cover classes. Roads layers from the San Carlos Apache Reservation, Tonto and Apache-Sitgreaves National Forest, and other sources were buffered by 200 m in ArcMap 10.5 (ESRI 2016) to constrain the accuracy assessment points to areas near a road. The map was broken into six regions based on the distribution of major roads and land management/ownership regions. Random points were generated for each vegetation class for each region with a minimum allowed distance between points of 100 m. Some classes, such as Spruce-Fir Forest, are less spatially extensive and therefore not represented in all regions; this resulted in 74 class/region combinations. We attempted to create 10 points in each of the 74 class/region areas; however, because of the small spatial extent of some classes as well as the minimum allowed distance constraint, this was not always possible. Eight of the 74 class/region combinations received $<10$ random sampling points. This resulted in a total of 684 possible field sampling points.

Field work was performed September 18 to 22, 2017, and October 31 to November 1, 2017, by driving along well-maintained roads to predetermined accuracy assessment points. Allowing for a suitable area to pull off, we attempted to locate the point to collect field data (Fig. 4). If a fence line (or a steep drop-off) impeded access, a “subjective point” was placed as near to the actual point as possible. Other subjective points were added along the route, often within broad expanses of a single vegetation type, or in unique vegetation “islands” within a more continuous cover. A total of 95 points were sampled. At each site, a GPS point and photographs in the four cardinal directions (north, east, south, and west) were collected and information on the local vegetative cover were documented, including the major vegetation species and their apparent abundances, and overall vegetation type based on one of the mapped vegetation classes. General vegetation of areas bordering the site was noted as well, as was evidence of fire and tree mortality/regeneration.

Fig. 4

Representative images of some vegetation classes from the study area: (a) Sonoran and Chihuahuan Desertscrub (Sonoran), (b) Plains and Semidesert Grassland, (c) Pinyon-Juniper-Oak Woodland, and (d) Ponderosa Pine Forest.

2.3.2.

2001 to 2004 SWReGAP field campaigns

Over the course of 2−1/2 field seasons, approximately 93,000 reference samples were collected directly, or obtained from other contemporary projects, for the SWReGAP land-cover mapping effort in the southwestern states of Arizona, Colorado, Nevada, New Mexico, and Utah.⁷ Of those, 933 were in the study area (Fig. 5). In the study area, the vast majority of field plots were field collected by traveling easily navigable roads and selecting plots that were at least 1 hectare in size with stand homogeneity.⁷ Ocular estimates included percent cover of dominant species and abiotic site characteristics. GPS coordinates were recorded.

Fig. 5

Sites where vegetation and land cover field data were collected.

2.4.1.

Statistical analysis

The unsupervised classification map was compared with 2017 field data and SWReGAP field data. Additionally, the recoded SWReGAP map was compared with 2017 field campaign data. The SWReGAP map was also compared with the SWReGAP field data to determine an upper boundary for this map’s accuracy over the study area. Finally, the supervised classification map was compared with both the SWReGAP and 2017 field campaign data to estimate upper boundaries for this map’s accuracy. Accuracy assessments were developed in R³³ using the caret package.³⁴ The error matrices for the SWReGAP map to SWReGAP field data, and supervised classification are not shown as the field data (100%) were used as the training data to produce these maps.

3.2.1.

Unsupervised classification

The results of the unsupervised classification had an overall accuracy of 53% when compared with 2017 field data (Table 3). In general, CA was higher than PA. Agriculture (Ag) and Pinyon-Juniper-Evergreen Oak Woodland were the most accurate classes overall with both CA and PA $>0.50$. Madrean Pine-Oak Forest/Woodland was the least accurate class with both CA and PA falling below 0.25.

Table 3

Error matrix comparing the unsupervised classification results (expected) to the 2017 field data (reference).

The unsupervised classification map was assessed using SWReGAP field data (Table 4) and compared with the accuracy assessment using the 2017 field data. Overall accuracy was similar at 52% with mixed results in both CA and PA. Not all classes can be compared between the 2017 field data and the SWReGAP field data since the SWReGAP field data are missing any Agricultural observation points and the 2017 field data is missing Spruce-Fir Forest. All of the classes saw changes in either CA or PA and most classes saw changes in both. Aspen Forest and Pinyon-Juniper-Evergreen Oak Woodland were lower in both CA and PA, while Mixed Conifer Forest, Ponderosa Pine Forest, and Madrean Pine-Oak Forest/Woodland were higher in both CA and PA. Interior chaparral, Montane Grassland, and Plains/Semidesert Grassland had higher PA but lower CA. Sonoran/Chihuahuan Desertscrub had lower PA but higher CA. Riparian Forest/Woodland and Mesquite Upland Scrub maintained their CA but had lower PA. Montane Grassland, Plains/Semidesert Grassland, and Sonoran/Chihuahuan Desertscrub were the most accurate (PA and $\mathrm{CA}>0.50$), whereas Spruce-Fir Forest was the least accurate (PA and $\mathrm{CA}<0.25$).

Table 4

Error matrix comparing the unsupervised classification map (expected) to the SWReGAP field data (reference).

3.2.2.

SWReGAP

The SWReGAP map was assessed using 2017 field data (Table 5). Overall accuracy was 44%, a reduction of 9% from the unsupervised classification map assessed with 2017 field data with mixed results in both CA and PA. Madrean Pine-Oak Forest/Woodland was the only class to improve in both CA and PA. Aspen Forest, Ponderosa Pine Forest, Montane Grassland, and Mesquite Upland Scrub showed decreases in both CA and PA. Pinyon-Juniper-Evergreen Oak Woodland had higher CA but lower PA, whereas Interior Chaparral had lower CA but higher PA. Both Mixed Conifer Forest and Plains/Semidesert Grassland had higher CA, whereas Riparian Forest/Woodland had higher PA. Sonoran/Chihuahuan Desertscrub had lower CA. Just as in the unsupervised classification map, agriculture and Pinyon-Juniper-Evergreen Oak Woodland were the most accurate classes, whereas Madrean Pine-Oak Forest/Woodland was the least accurate.

Table 5

Error matrix comparing the SWREGAP map with generalized classes (expected) to the 2017 field data (reference).

Lowry et al.⁷ reported a validation accuracy of 61% for the entire five state mapping area. This validation accuracy was based on an accuracy assessment of an intermediary product and the accuracy of the final SWReGAP map is likely higher than the validation accuracy. To find an upper boundary for accuracy of the final SWReGAP product over the study area, the SWReGAP map was compared with the SWReGAP field data. The resultant accuracy is overestimated because the test data were taken from the training data³⁰ but it sets the upper boundary of accuracy at 75%. The overall accuracy of the SWReGAP map was roughly considered to be between these boundaries (61% minimum—measured over the five state mapping areas, and 75% maximum—measured over the study area). The error matrix comparing the recoded SWReGAP map to the SWReGAP field data is not presented here, because 100% of the field data were used to create the map, and it does not represent a competent estimate.

3.2.3.

Supervised classification

The results of the supervised classification had an overall accuracy of 43% when compared with 2017 field data, and 47% when compared with the SWReGAP field data. The error matrices are not presented here, because 100% of the field data were used to create the map, and they represent an unrealistic upper limit to the accuracy estimates. Even so, the supervised classification map has very low accuracy and probably requires a significant amount of work with the selection of training sites before being useful as an accurate map product. Moreover, this classification was only meant to explore the initial results of using a supervised classification on the field data to see if that is a reasonable avenue for future mapping.