Predictive Mapping of Low-Density Juniper Stands in Prairie Landscapes of the Northern Great Plains

ABSTRACT Expanding distributions of native juniper species have had significant ecological and economic impacts on prairie ecosystems of the Great Plains. Juniper encroachment reduces rangeland production by decreasing herbaceous biomass and affecting natural ecosystem functions as it alters other native plant communities, microclimates, and soils. Juniper distribution maps are needed to support proactive management, but they often underestimate the extent of low-density juniper stands. Our objectives were to extend a previous juniper mapping study by 1) fitting a predictive ecological model for low-density (< 15% fractional cover) juniper stands and assessing the classification accuracy, 2) determining the habitat variables that had the strongest associations with low-density juniper, and 3) applying the model to map low-density juniper stands, where proactive management has the greatest potential for stopping further juniper encroachment. The study area included counties bordering the Missouri River in southeastern South Dakota and northeastern Nebraska covering approximately 23 000 km2. Environmental predictors included seed source distance and density, as well as topography, climate, soils, and land use variables. Areas of low-density juniper were identified by visual interpretation of sample plots from digital aerial photography. We used a machine-learning approach to classify low-density juniper with the random forests algorithm. Model accuracy was high with an area under the receiver operating characteristic curve of 0.884. Variables related to seed sources were the most important predictors, and precipitation, slope angle, and the local intensity of human land use also had substantial influences. A previous map based on Landsat imagery identified 209 968 acres (84 971 ha) as juniper with in the study area, and this study found an additional 430 648 acres (174 277 ha) classified as low-density juniper stands. These results can provide agencies and land managers with more accurate information about the distribution of juniper, and the underlying techniques can be extended to map woody plant encroachment in other areas.


Introduction
In many parts of the North American Great Plains, expansion of woody plants threatens the natural functions of grassland ecosystems ( Knapp et al., 2008 ;Van Auken 2009 ;Ratajczak et al., 2012 ). Notably, eastern redcedar (Juniperus virginiana) and Rocky Mountain juniper (Juniperus scopulorum) have affected carbon storage, soil characteristics, and plant communities within the prairie ecosystems of the Great Plains ( Norris et al., 2001 ;McKinley and Blair 2008 ;Pierce and Reich 2010 ). In the central United States, the encroachment of these species (hereafter referred to collectively as "juniper") resulted in the conversion of 205 0 0 0 acres (82 Management of encroaching woody plants is crucial, especially when attempting to control a quickly expanding juniper footprint. Proactive and reactive management are two approaches for controlling juniper ( Simonsen et al., 2015 ). Proactive management measures are implemented before juniper has established or is in vulnerable seedling and sapling stages and can include grazing, haying, and low-intensity prescribed burning ( Wilson and Schmidt 1990 ;Smith 2011 ;Simonsen et al., 2015 ). In contrast, reactive management practices are implemented where juniper is already well established. Mechanical removal by timber cutting, herbicides, and high-intensity prescribed burning are commonly used methods for reducing juniper on established sites ( Wilson and Schmidt 1990 ;Smith 2011 ;Simonsen et al., 2015 ). Reactive management of juniper is costly and time consuming, and it becomes less effective as stand density and tree size increase ( Buehring et al., 1971 ;Ortmann et al., 1998 ;Bidwell et al., 2002 ). Therefore, proper planning is essential so that proactive management methods can be used in a timely manner.
Juniper maps can aid in targeting areas for proactive management, particularly if they can identify low-density sites before juniper is fully established. Newly established juniper can mature and begin producing seeds within 6 yr ( Twidwell et al., 2021 ) and can reach 15% canopy cover within 10 −15 yr . Remote sensing is a common way of obtaining spatial information about forest distribution and conditions that can be applied for management purposes ( Franklin 2001 ;Giri 2012 ). Juniper distributions have been mapped with various remotely sensed data sources including very high spatial resolution (VHSR) aerial imagery ( Anderson and Cobb 2004 ;Poznanovic et al., 2014 ), multispectral imagery ( Sankey and Germino 2008 ;Kaskie et al., 2019 ), hyperspectral imagery ( Wylie et al.,20 0 0 ), and multisource fusion of data from active and passive sensors ( Sankey et al., 2010 ;Wang et al., 2017 ). Although these studies have successfully mapped mature juniper stands, mapping areas with sparse coverage during the establishment phase has been more challenging. Wang et al., (2017) found that < 30% of juniper sites in Oklahoma with 10 −20% tree cover was correctly classified using multitemporal synthetic aperture radar (SAR) and Landsat imagery. Kaskie et al., (2019) mapped junipers in South Dakota and Nebraska and similarly found that < 50% of juniper sites with 10 −20% tree cover and < 15% with 1 −10% cover were correctly classified using snow-covered winter Landsat imagery. Because of these challenges, tree cover estimates in rangelands have been largely absent from operational, continental-scale product such as the National Land Cover Dataset ( Kaskie et al., 2019 ). More work is needed to improve the detection of areas in the early stage of juniper establishment where proactive management aimed at preventing juniper encroachment is most effective ( Simberloff 2003 ;Yokomizo et al., 2009 ).
Several approaches have been used to map low-density tree cover in rangeland ecosystems. Techniques that use very-highresolution imagery from digital aerial photographs ( ≤ 1-m pixel size) combined with object-based feature detection methods such as spatial wavelet analysis can detect individual trees ( Strand et al., 2006 ;Falkowski et al., 2017 ). However, high data volumes and requirements for specialized software can make these methods challenging to apply. When used with machine learning techniques and large training datasets, Landsat can be used to predict fractional tree cover as a continuous variable ( Allred et al., 2021 ). Landsat data are widely used for vegetation mapping because of their free availability, global coverage, and large historical archive. However, detecting low-density juniper stands with Landsat is difficult because the pixel size of Landsat (30 m) typically encompasses a mixture of juniper and other background materials. Alternatively, predictive vegetation mapping ( Franklin 1995 ) can be used to supplement remotely sensed maps by identifying areas suitable for juniper establishment based on environmental variables. Predictive modeling of juniper encroachment in riparian habitats has found that juniper establishment is influenced by soils, floods, and the historical effects of flow regulation ( Greene and Knox 2014 ;Illeperuma and Dixon 2021 ). However, riparian areas are only one of the many habitats where juniper occurs ( Lawson 1990 ;Noble 1990 ), and larger-scale predictive models are needed for broader applications.
The goal of this study was to predict the probability of lowdensity juniper ( < 15% cover) in areas of southeastern South Dakota and northeastern Nebraska that were not identified as containing juniper in a previous Landsat-based classification ( Kaskie et al., 2019 ). This product reliably classified juniper pixels with fractional cover > 15% (89 −95% overall accuracy) but had very low detection probabilities for sparser juniper cover. Thus, an additional modeling study focused on low-density juniper was undertaken to better characterize the extent of juniper encroachment across the study area. Specific objectives were to 1) fit a predictive ecological model for the probability of low-density juniper occurrence and assess the resulting classification accuracy, 2) determine habitat variables that have the strongest associations with low-density juniper stands, and 3) extrapolate predictions across the study area to identify areas where low-density juniper is most abundant.

Study area
Our study area covered 14 contiguous counties (nine counties in southeastern South Dakota and five counties in northeastern Nebraska [ Fig. 1 ] all bordering the Missouri River. This area has a Köppen climate classification of Dfa; humid continental climate consisting of warm to hot summers and cold winters ( Kottek et al., 2006 ), with an annual temperature range of 6 −11 °C and an annual average precipitation of 498 −796 mm. Native vegetation consists of mixed-grass prairie species such as little bluestem (Schizachyrium scoparium), big bluestem (Andropogon gerardii), western wheatgrass (Pascopyrum smithii), sideoats grama (Bouteloua curtipendula), and green needlegrass (Nassella viridula). Woodlands are primarily found near drainages and riparian lowlands with the exception of small groves scattered across the prairie uplands. The most common deciduous species includes the plains cottonwood (Populus deltoids) with the occasional green ash (Fraxinus pennsylvanica) and American elm (Ulmus americana). Juniper species such as Rocky Mountain juniper (Juniperus scopulorum) and eastern redcedar (Juniperus virginiana) are also common ( Barker and Whitman 1988 ). Steeply sloped drainages disrupt a flat to rolling topography composed largely of agriculture (48%) and herbaceous grasslands (39%), producing a fragmented landscape ( Wimberly et al., 2018 ). Primary land uses include the agricultural production of corn, soybeans, and wheat, as well as cattle ranching ( Wimberly et al., 2017 ).

Juniper training and validation data
We used a stratified random sampling design to collect data from photo-interpreted juniper plots over a range of juniper densities. We allocated four strata, which included closed canopy woodlands, buffered closed canopy woodlands, planted shelterbelts (narrow rows of trees planted for wind protection), and nonwoodland areas. We digitized the closed canopy woodlands and planted shelterbelts in ArcGIS 10.5 following the guidelines outlined in Bauman et al., (2016) using 60-cm very high spatial resolution (VHSR) data from the National Agricultural Imagery Program collected in 2014 and 2016. To obtain samples of intermediate juniper densities, we placed a 90-m buffer around the digitized closed canopy stratum. The remaining nonwoodland areas had relatively few junipers.
We generated 4 308 random points including 1 282 points in closed-canopy woodlands, 671 points in shelterbelts, 1 210 points in buffered closed-canopy woodlands, and 1 145 points in nonwoodland areas. Each random point was referenced to a Landsat pixel by converting it to a 30 × 30 m polygon and snapping to the Landsat 8 pixel grid. We visually estimated juniper cover in each polygon using a combination of VHSR imagery, which included National Agricultural Imagery Program 2016 and other sources of winter imagery accessed through Google Earth from 2013 to 2017. All estimates were made by the same interpreter for consistency. These data were originally used for accuracy assessment of the remotely sensed juniper map developed by Kaskie et al., (2019) . To increase the size of the dataset for training the juniper encroachment model, we sampled an additional 500 points in closed canopy woodlands and 500 points in buffered closed canopy woodlands. The final training dataset included all points with juniper cover ≤ 15% or no juniper cover. There were 865 lowdensity juniper points (1 −15% cover) and 2 468 juniper absence points ( Fig. 2 ).

Ecological predictor variables
We examined 15 predictor variables ( Table 1 ) that were hypothesized to be associated with juniper encroachment. These variables measured topography (slope and aspect); climate (mean annual temperature and total annual precipitation); soils (per-cent sand, percent silt, percent clay, soil available water storage, depth to restricted layer, root zone depth, and soil drainage class); land use in the surrounding landscape (percent agriculture); seed sources (distance to juniper and percent juniper in the surrounding landscape); and flooding potential (distance to surface water). We processed these variables in Ar-cGIS 10.5 and resampled them all to a 30-m spatial resolution raster.
We used the National Elevation Dataset with 30-m spatial resolution to derive topographic indices. The National Elevation Dataset, created by the US Geological Survey, was accessed through the US Department of Agriculture Geospatial Data Gateway. Percent slope represented the physical gradient of the landscape and aspect indicated the direction the slope faced. We reclassified aspect into nine categories including flat (0% slope); cardinal directions (N, E, S, W); and ordinal directions (NE, SE, SW, NW).
We used the PRISM (Parameter-elevation Relationships on Independent Slopes Model) dataset to extract 30-yr normals (1981 −2010) of climate variables, including average annual minimum temperature, average annual maximum temperature, and average annual precipitation. We used the minimum and maximum temperatures of each pixel to derive the mean annual temperature.
We used the gridded Soil Survey Geographic (gSSURGO) database obtained through the Geospatial Data Gateway to extract multiple soil characteristics. Gridded SSURGO contains the same data provided in the US Department of Agriculture Natural Resources Conservation Service SSURGO database in a 10-m spatial resolution raster format. Soil texture is represented as three sepa-  rate factors including percent sand, percent silt, and percent clay ( Wang et al., 2018 ). A combination of these three factors can be used to define the soil classification within the soil profile. We used the available water storage estimate, which represented the volume (mm) of plant available water the soil can store within a 0to 150-cm soil profile. We used depth to restricted layer as a measure of the distance (cm) within the soil profile showing any restricting features that may constrain root growth or the movement of water and air. We used the root zone depth as a measure to the depth at which plants, root systems can effectively obtain nutrients and water. We also used soil drainage class, which is characterized into seven classes (excessively drained to very poorly drained) and reflects the natural frequency and duration of wet periods for the soil.
We obtained land cover data at a 30-m spatial resolution from the 2016 Cropland Data Layer (CDL). The CDL is produced annually and contains a remote-sensing based classification of different crop types. We used the 2016 CDL dataset to extract the percent agricul-ture within a 5 × 5-pixel window surrounding an individual 30-m pixel. We denoted agriculture as any human use of the landscape for cultivated crops, alfalfa, or other hay fields.
We used a 30-m classified juniper distribution map from a previous study ( Kaskie et al., 2019 ) to estimate distance to seed sources. This map was based on two snow-covered Landsat 8 images that encompassed the study area. Juniper pixels were classified using a matched filtering approach in ENVI version 5.4 software (Exelis Visual Information Solutions, Boulder, CO). Matched filtering produces continuous values that indicate the similarity of each pixel to the spectral signature of pure juniper, and these values were thresholded to produce a classification of juniper presence/absence. The resulting map was validated using canopy cover data obtained from visual interpretation of very high resolution imagery. When juniper cover was above 50%, the detection probability was ≥ 90%. However, when juniper cover was lower than 20%, the true positive rate was < 50%. We used this juniper presence-absence map to compute Euclidian distance to the near- Distance to surface water was measured using the 1:24 0 0 0 scale US Geological Survey National Hydrography Dataset. We used the National Hydrography Dataset to map all surface water features within the study area. The Euclidean distance from each pixel to the closest surface water feature was then calculated.

Juniper encroachment model
The random forests model is an ensemble-learning algorithm ( Breiman 2001 ). This approach has been used in developing susceptibility models ( Ismail et al., 2010 ;Wang et al., 2015 ;Youssef et al., 2016 ;Chen et al., 2017 ) and has been applied in numerous predictive mapping applications and across multiple disciplines ( Prasad et al., 2006 ;Belgiu and Dr ȃgu ¸t 2016 ;Biau and Scornet 2016 ). We used random forests to predict the probability of lowdensity juniper presence given the set of observed environmental variables.
During each iteration of the random forests algorithm, the data are separated into two groups. The in-bag dataset consists of a bootstrap sample with approximately two thirds of the observations in the entire dataset and is used to train the algorithm. The in-bag samples are used to train a decision tree without pruning, where each split is based on one of a randomly selected set of predictor variables (Mtry). The out-of-bag (OOB) dataset contains the remaining observations and is used to generate predictions and estimate prediction error. Because the OOB predictions and associated error estimates are independent of the training dataset, they produce results that are similar to a leave-one-out cross-validation ( Hastie et al., 2008 ). This process is repeated multiple times until a user-defined number of decision trees (Ntree) is reached. Each decision tree rule is used to cast a vote for a predicted class with the maximum number of votes becoming the final classification ( Breiman 2001 ). The probability associated with a class is calculated as the total proportion of votes across all trees.
We developed random forests models using R statistical software version 3.4.1 ( R Core Team 2021 ) and the randomForest package ( Liaw and Wiener 2002 ). We set the Ntree parameter to 500, as trial runs with this value showed a stabilization of the errors before this maximum tree number was attained ( Fig. 3 ). Increasing Ntree to 10 0 0 or more had minimal effects on the results, so the parameter value of 500 was retained for computational efficiency. The Mtry parameter was set as the default, which is the square root of the total number of the predictor variables used in the model.

Model assessment
The accuracy of model predictions was assessed by comparing observations with OOB predictions that were independent of the observations, providing results similar to a leave-one-out crossvalidation ( Hastie et al., 2008 ). We also computed accuracy statistics using a data-splitting approach, where 70% of points were used as training data and 30% were reserved as a validation dataset.
We used the Hosmer-Lemeshow goodness-of-fit test to evaluate model calibration. This test is commonly used in risk and susceptibility modeling ( Bai et al., 2010 ;Catry et al., 2010 ;Fang et al., 2013 ). It categorizes subgroups (referred to as deciles of risk) and performs a Pearson chi-square statistic on the expected and observed frequencies. A close relation of the expected and observed frequencies reflects a model with good fit ( Hosmer and Lemeshow 20 0 0 ).
We evaluated model predictions using receiver operating characteristics (ROC). ROC examines multiple classification cutpoints by plotting the probability of detecting a true positive (sensitivity) against a false positive (1-specificity) over a range of cutoffs for classifying low-density juniper based on the proportion of votes from the random forests ensemble ( Hosmer and Lemeshow 20 0 0 ). The area under the curve (AUC) measures discrimination or the likelihood the model will predict target versus nontarget sites. Hosmer and Lemeshow (20 0 0) suggest that AUC values represent discrimination as being none (0.5), acceptable (0.7 −0.8), excellent (0.8 −0.9), or outstanding (0.9 −1.0). We investigated the model performance of a full model and compared it with the performance of models with manually removed variables to eliminate unnecessary variables while improving our final model computation time ( Plant 2012 ). We began by removing strongly correlated predictor variables followed by a trial-and-error process in which we independently excluded the remaining variables. We found the final model to have as much predictive power as our full model with a minimal effect on the modeling error.

Model accuracy
The full model with all predictor variables achieved an out-ofbag (OOB) error rate of 15.3% (84.7% accuracy). After removing correlated and redundant variables, the final model contained eight predictor variables: distance to juniper, percent juniper in the surrounding landscape, total annual precipitation, percent slope, percent agriculture, soil available water storage, mean annual temperature, and aspect. The OOB estimate of error rate for the final model was 15.2% (84.8% accuracy), showing that the reduced model had as much predictive accuracy as our full model. The model had an OOB class error rate of 7.2% (92.8% accuracy) for nonjuniper (juniper absence) and 36.9% (63.1% accuracy) for lowdensity juniper (juniper presence).
We used the Hosmer-Lemeshow goodness-of-fit test to evaluate the calibration of our final model. The P value for the Hosmer-Lemeshow test was > 0.05 ( χ² = 5.8565, df = 8, P value = 0.6633), indicating that there was not a statistically significant difference between the expected and observed values and the final model was well calibrated. We also used the ROC to evaluate model predictions, and accuracy was high with an AUC of 0.884 ( Fig. 4 ). The Hosmer-Lemeshow test and AUC results indicated that the final low-density juniper model was a good predictor of low-density juniper distributions.
Using the data-splitting validation approach, the overall error rate was 14.6% (85.4% accuracy). The class error rates were 6.3% (93.7% accuracy) for nonjuniper (juniper absence) and 38.6% (61.4%   5. Mean decrease accuracy of final eight low-density conditioning factors (listed in descending order) assigned by our random forests model. The eight predictor variables listed are distance to juniper, percent juniper in the surrounding landscape, total annual precipitation, percent slope, percent agriculture, soil available water storage, mean annual temperature, and aspect. Variable abbreviations are defined in Table 1 . accuracy) for low-density juniper (juniper presence). The AUC was 0.872. These results were similar to those based on the OOB predictions, confirming that accuracy remained high even when predicting observations that were outside of the training dataset.

Environmental predictors
The importance of the eight predictor variables is shown in Fig. 5 by descending order of mean decrease in accuracy. Distance to juniper was observed to have the highest conditional importance (51.71) with a strong negative relationship between lowdensity juniper and the distance to a seed source ( Fig. 6 a). The highest probabilities of low-density juniper ( > 0.25) were within 100 m of a mature juniper seed source, and the probability decreased rapidly out to approximately 300 m. In the field, we observed that low-density junipers located longer distances ( > 1 0 0 0 m) from a seed source were typically planted immature junipers. Percent juniper within a 5 × 5-pixel window (36.81) had the second highest importance with a strong positive relationship (see Fig. 6 b), indicating that there was an increase in low-density juniper where multiple mature seed sources were in close proximity. Total annual precipitation (28.99) was observed to have a negative relationship (see Fig. 6 c) and was followed by percent slope (26.40) with a positive relationship (see Fig. 5 d) and percent agriculture (25.11) with a negative relationship (see Fig. 6 e). The three variables with lowest importance included soil available water storage (15.67) with a negative relationship (see Fig. 6 f), mean annual temperature (13.14) with a positive relationship (see Fig. 6 g), and aspect with the lowest importance (6.86) and no clear relationship (see Fig. 6 h).

Juniper encroachment maps
We produced a final low-density juniper map by grouping the continuous probability values into five Jenks natural breaks groups: 0.0 0 0 −0.075, 0.075 −0.176, 0.176 −0.380, 0.380 −0.694, and 0.694 −1.0 0 0 ( Fig. 7 ). By adjusting the probability thresholds for classifying low-density juniper, map users can control the tradeoff between detection probability (sensitivity) and correctly predicting juniper absence (specificity) as shown by the ROC curve in Fig. 4 . For example, using probabilities > 0.38 (the two highest probability classes in Fig. 7 ) as a cutoff for low-density juniper prediction resulted in a sensitivity of 0.71 and a specificity of 0.88.
Visualizing the probabilities as a continuum provides information about the relative densities of low-density juniper pixels across the landscape. The probability map was able to capture areas of low-density, encroaching juniper that were not classified as juniper in the Kaskie et al., (2019) Landsat-based map ( Fig. 8 ), with few low-density juniper pixels where probabilities were < 0.176 and increasing juniper density over probabilities ranging from 0.176 to 1. The summaries in Table 2 show that 430 648 acres (174 277 ha) had a predicted juniper probability > 0.38 in addition to the 209 968 acres (84 971 ha) of juniper originally classified in Kaskie et al., (2019) . Thus, the footprint of juniper encroachment was considerably larger once the extent of low-density juniper stands was incorporated.

Discussion
Satellite remote sensing has been shown to accurately predict the spatial distribution of relatively dense, mature juniper stands but is less effective at predicting locations where juniper is still establishing and cover is relatively low Kaskie et al., 2019 ). In this study, we showed that it is possible to map lowdensity stands with juniper cover < 15%) using a model with predictor variables for topography, climate, soils, land use, and seed source availability. The model had a relatively high prediction accuracy (AUC = 0.88) and highlighted areas with low-density juniper that were not identified by our previous Landsat-derived map of juniper presence/absence. The predictive map of low-density juniper highlights locations where encroaching junipers are still primarily in seedling and sapling stages and are potentially suitable for proactive management, such as low-intensity burning, haying, and grazing, which can prevent the establishment of mature juniper that will be more difficult to remove.
Although the original Landsat-based map ( Kaskie et al., 2019 ) was not effective at predicting low-density juniper stands, it allowed us to extract predictor variables that greatly influenced the juniper encroachment model. This was expected, as low-density ju- Fig. 6. a-h, Partial dependence plots for eight predictor variables used in the final random forests model for predicting low-density juniper. Variable abbreviations are defined in Table 1 . niper stands are usually associated with seed dispersal from established juniper sites ( Holthuijzen et al., 1987 ;Yao et al., 1999 ). Holthuijzen and Sharik (1985) showed that within abandoned fields, eastern redcedar density decreased as the distance from a seed source increased. This result supports our finding that the probability of observing low-density juniper was highest < 100 m from established juniper and decreased with increasing distance from established juniper. Additionally, Owensby et al., (1973) found highest establishment of eastern redcedars within fields that were already heavily invaded by junipers. We similarly found that an increase in the distribution of established juniper within a 5 × 5pixel window ( ≈60-m radius) increased the probability of lowdensity juniper occurrence. The fact that both variables were important in our model emphasizes that seed source density and distance to seed source are important predictors of juniper establishment.
Predictor variables related to juniper cover had the highest conditional importance within our model, but climatic and topographical factors were influential as well. The 30-yr normal of total annual precipitation had the third highest importance. We    Owensby et al., (1973) similarly found that for every additional inch (2.54 cm) of precipitation, the juniper invasion rate decreased by 0.2 trees per acre (0.49 trees per ha). We also found a positive relationship with slope angle, in agreement with previous findings that Eastern redcedar can be associated with moderate to steep slopes ( Anderson 2003 ), which have shallow soils and provide less competition from other plants and more protection from fire ( Bryant 1989 ;Pierce and Reich 2010 ). The negative relationship with agricultural land use in a 5 × 5-pixel window suggests that activities such as active cultivation, haying, and grazing prevent juniper establishment and reduce the potential for spread into the surrounding landscape. Fire is another type of disturbance that can greatly reduce juniper abundance, particularly at the seedling and sapling stages. We did not have reliable maps of fire history for our study area, but where these data are available they could be incorporated as an additional predictor variable.
Variable importance ranking for our remaining variables fell considerably below the highest-ranked variables, yet removal resulted in a slight reduction in model performance and they were therefore retained in the final model. Low-density juniper occurrence decreased with increasing soil water storage but then increased again at water storage values > 30%. Similar to the effects of slope, very low and high water storage values may indicate sites where competition from other plant species is minimized. Low-density juniper was also increased slightly at warmer temperatures, although there was only a small variation in mean annual temperature (46 −50 °F, 7.8 −10 °C) across the study area. Aspect was the least important environmental factor and did not have a clear relationship with low-density juniper. Some studies have shown aspect to influence juniper establishment ( Lawson 1990 ;Schmidt and Stubbendieck 1993 ), whereas others have found no effects of aspect on juniper establishment ( Tunnell et al., 2004 ).

Implications
Previously developed juniper classification maps based on Landsat imagery allowed us to characterize current distributions of established juniper ( Kaskie et al., 2019 ). The predictive model of low-density juniper developed in this study provided additional information about locations with that were not detectable with Landsat imagery. These low-density juniper stands were abundant within the study area. Adding the locations where low-density juniper was predicted with high probability ( > 0.38) to our previous Landsat-based map more than tripled the prediction of juniper extent from 84 971 to 259 248 ha. If these additional effort s were not made to map low-density juniper stands, the extent of juniper encroachment would be greatly underestimated.
The greater time and cost of managing juniper encroachment with increasing size and density emphasizes the importance of detecting and mapping these low-density areas ( Ortmann et al., 1998 ;Bidwell et al., 2002 ). Emerging management frameworks for reducing woody plant encroachment are inherently spatial ( Twidwell et al., 2021 ), and early detection and rapid response following woody plant recruitment has been identified as a core management philosophy. Predictive maps such as the one generated in this study can support these effort s by highlighting places in the landscape that contain low-density juniper abundance to prioritize them for proactive management. Estimates of the total area of low-density juniper at a particular location can be valuable for budgeting and determining equipment and personnel needed for management. Knowing the locations with low-density juniper can help to direct on-the-ground surveillance effort s and determine the management activities that will be most suitable at particular sites. Because seed sources are one of the most important constraints to juniper establishment, maps of the current juniper distribution are essential for understanding where dispersal and recruitment will occur in the future.
The techniques that we have used are distinctive from, but complementary to, other studies that have mapped rangeland tree cover as continuous variables ( Falkowski et al., 2017 ;Allred et al., 2021 ) and detected rangeland regime shifts through spatial analysis of multitemporal land cover data ( Uden et al., 2019 ). All the geospatial data that we used to develop the low-density juniper map are freely available for the entire United States and many other parts of the world. In addition, the modeling techniques used in the study have been implemented in a variety of freely accessible software platforms such as R ( R Core Team 2021 ) and Google Earth Engine ( Gorelick et al., 2017 ). Therefore, the approach that we developed and tested can be used by a wide variety of stakeholders, easily updated in the future, and extended to other rangeland systems where encroachment of juniper and other woody species are management concerns.