1.Question asked and data used: 

How does juniper density vary with slope and aspects characteristics (as indicators of soil moisture)?

This study site for this analysis is the Camp Creek Paired Watershed Study (CCPWS), located in a semiarid region of central Oregon. This research focuses on two watersheds: one in which most juniper was removed in 2005-2006, referred to as the treated WS, and one watershed in which mature juniper is the dominant overstory species, referred to as the untreated WS. Belt transects (3m wide by 30m long) which recorded the number, height, and canopy cover of juniper within the treated WS were used for analysis. The number of juniper found in each transect was represented as a single point at the beginning of each transect. Additionally, NAIP imagery and a 30m DEM were used to assess juniper density and slope and aspect characteristics across both watersheds and in the surrounding areas.

2. Approach and tools used:

Analysis was conducted using R Studio, ArcGIS Pro, and Excel. A Chi-squared test and logistic regression were calculated in R, followed by Geographically Weighted Regression (GWR) in ArcGIS Pro. Excel was used to calculate the percentage of pixels classified as juniper within each buffer.

3. Steps used in analysis:

     A. Data preparation

              1. The following steps were conducted in a previous exercise: 1) a 30m DEM was used to extract slope and aspect characteristics, 2)slope and aspect values were each divided into nine groups, representing characteristics generally associated with increased soil moisture (e.g., flat slopes with northerly aspects were rated highest), 3)using the raster calculator the values of slope and aspect were averaged ((slope category+aspect category)/2). The raster created using this process rated each grid square from 1-9 (henceforth referred to as “rating model”). The rating model is shown below.

             2.Classification of NAIP imagery (also from a previous exercise). The NAIP imagery was classified using Support Vector Machine supervised classification. (Please note: the resolution of the NAIP imagery makes identification of juniper saplings difficult, therefore conclusions from the data should be made with caution.) The NAIP raster was resampled to the same resolution as the 30m DEM. The pixels were reclassified into a binary raster, with each pixel being either “juniper” or “not juniper”. The resampled NAIP imagery is shown below.

            3.The classified NAIP raster and rating model were projected to NAD 1983 UTM zone 10. Both rasters were masked to cover the same extent.

            4. A table was created which represented the number of juniper found per transect as a point. Using the “extract multi values to points” tool in ArcGIS the corresponding value from the rating model was assigned to each transect point.

             5. In a previous exercise, a buffer of 50m was created surrounding each of the belt transects (treated WS only). The percentage of juniper (pixels classified as juniper/total pixels) was calculated in excel and added to the belt transect table created previously.

             6. Using the “create random points” tool in ArcGIS, 96 random points were created within the study area (treated WS, untreated WS, as well as surrounding area). These points were a minimum of 150m apart.

             7. A 50m buffer was created surrounding each of the 96 random points. The percentage of juniper in each buffer was calculated as indicated above. These values were added (using the join function) to a table with the random points.

      B. Analysis

             1. The binary classified NAIP raster and rating model were loaded into R Studio and reviewed (addresses have been abbreviated):

library(raster)
library(sp)
library(rgdal)

apslop<-raster(“C:/…AspSlop_projected_clip.tif”)
juniperonly_NAIP<-raster(“C:/…NAIP_resampled_clip1.tif”)

apslop
juniperonly_NAIP
plot(apslop)
plot(juniperonly_NAIP)

s2<-stack(apslop,juniperonly_NAIP)
v2<-data.frame(na.omit(values(s2)))

             2. Chi-squared test was conducted for the rating model and classified NAIP raster (R):

chisq.test(apslop[],juniperonly_NAIP[])

             3. Logistic model, summary of results and associated plots (R):

glm.fit<-glm(juniperonly_NAIP[]~apslop[],data=v2,family=binomial)
summary(glm.fit)

plot(glm.fit)
exp(coef(glm.fit))

              4. GWR was conducted in ArcGIS for the belt transects (in treated WS), using the juniper density in the 50m buffer as the dependent variable and the rating model as the independent variable. The distance band neighborhood type and golden search neighborhood selection were used.

              5. GWR was conducted in ArcGIS using the 96 random points (in treated WS, untreated WS, and surrounding areas). The same neighborhood parameters described above were used.

              6. Global Moran’s I was calculated using the residuals for both GWR analyses. Inverse distance and the Euclidean distance method were used for analysis.

4. Overview of Results:

      A. The Chi-squared test of the binary NAIP raster (indicating juniper versus non-juniper) and the rating model raster indicated a slightly significant relationship:

Pearson’s Chi-squared test data:

apslop[] and juniperonly_NAIP[]

X-squared = 28.624, df = 15, p-value = 0.01798

      B. The logistic regression model indicated a significant relationship between the aspect and slope classification (from the rating model) and the presence or absence of juniper (dependent variable). However, the difference between the null deviance and the residual deviance suggests that this may not be an appropriate model for predicting the presence or absence of juniper. The exponent of the coefficient indicated that for every increase in aspslop rating the probability of juniper being present increased by 1.06.

Call: glm(formula = juniperonly_NAIP[] ~ apslop[], family = binomial, data = v2) Deviance Residuals:

Min 1Q Median 3Q Max

-0.6934 -0.6576 -0.6317 -0.5984 1.9537

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -1.80732 0.11048 -16.359 < 2e-16 ***

apslop[] 0.05935 0.01948 3.046 0.00232 **

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 7425.5 on 7761 degrees of freedom

Residual deviance: 7416.1 on 7760 degrees of freedom

(656 observations deleted due to missingness)

AIC: 7420.1

Number of Fisher Scoring iterations: 4

(Intercept)      apslop[]

0.164094       1.061147

      C. The GWR of treated WS dataset indicated an adjusted R-squared value of 0.28. Also, the following warning was indicated for this analysis: “At least one of your local regressions had very limited variation after applying the weights. Please use caution when interpreting the results”.  This further suggests that this model may not be indicative of the explanatory factors related to juniper density. The GWR results are provided below.

Golden Search Results

Distance Band     AICc

292.6545      222.3773

884.3284      221.4705

518.6538      216.8374

658.3291      219.2658

432.3298      217.0552

572.0050      217.5695

485.6810      216.6200

465.3026      216.6782

———————

 WARNING 110259: At least one of your local regressions had very limited variation after applying the weights. Please use caution when interpreting the results.

——————————– Analysis Details ——————————-

Number of Features:                                                             33

Dependent Variable:                      BELTTRANSECTS_50M_EXCELTOTABLE.F__JUNIPER

Explanatory Variables:  C250M_ALLSITES_EXCELTOTABLE_PROJECT.ASPSLOP_PROJECTED_CLIP

Distance Band:                                                            465.3026

———————————————————————————

–——– Model Diagnostics ———-

R2:                              0.4918

AdjR2:                           0.2845

AICc:                          216.6782

Sigma-Squared:                  28.4635

Sigma-Squared MLE:              20.4670

Effective Degrees of Freedom:   23.7290

     D. The GWR performed for the larger dataset (points distributed across both watersheds and in the surrounding areas) indicated a similar R-squared value to the GWR based in the treated WS only. A similar warning regarding the residuals was indicated for this GWR analysis. The results of this GWR are displayed below.

Golden Search Results

Distance Band     AICc

1360.4243     -98.1728

3270.5186     -84.8038

2090.0154     -90.0118

2540.9275     -87.1310

1811.3364     -92.9095

1639.1033     -95.1341

1532.6574     -96.4953

1466.8702     -97.2349

1426.2115     -97.6318

1401.0830     -97.8502

1385.5527     -97.9775

1375.9545     -98.0539

———————

WARNING 110259: At least one of your local regressions had very limited variation after applying the weights. Please use caution when interpreting the results.

—————————- Analysis Details —————————

Number of Features:                                                     96

Dependent Variable:     RANDOMPTS_50MBUFFER_JUNIPER_EXCELTOTABLE.F_JUNIPER

Explanatory Variables:            RANDOM_POINTS_105.ASPSLOP_PROJECTED_CLIP

Distance Band:                                                   1375.9545

————————————————————————-

——— Model Diagnostics ———-

R2:                              0.2932

AdjR2:                           0.1943

AICc:                          -98.0539

Sigma-Squared:                   0.0190

Sigma-Squared MLE:               0.0167

Effective Degrees of Freedom:   84.3419

     E.The Global Moran’s I for the GWR residuals from the belt transects dataset (treated WS) indicated clustering (p<0.0001) while the GWR residuals for the larger dataset (both WS with surrounding area) indicated random distribution (p=0.327).

5. Discussion/critique of methods:

The analysis used here suggests that other explanatory factors (outside of the rating model created based slope and aspect) may be in place which influence the density of juniper. The R-squared values of this analysis improved compared to previous exercises but are still relatively low. Additionally it should be noted that the presence of juniper has been correlated to reduced soil moisture in some conditions (Lebron et al., 2007), therefore slope and aspect may not be sufficient representatives of soil moisture in these watersheds.

The intent of this analysis was to develop a potential workflow that can be used with other datasets. Caution should be used in drawing any conclusions from this analysis specifically. Additionally, caution should be taken in interpreting the results of this analysis due to the characteristics of the datasets used. For instance, the resolution of the NAIP imagery makes the detection of juniper saplings difficult. Therefore the classification results may not be an accurate representation of juniper density, particularly in the treated WS. However, this process may be applied to more expansive data to include higher-resolution imagery (such as imagery collected using UAVs, etc.).

Reference:

Lebron, I., Madsen, M. D., Chandler, D. G., Robinson, D. A., Wendroth, O., & Belnap, J. (2007). Ecohydrological controls on soil moisture and hydraulic conductivity within a pinyon-juniper woodland. Water Resources Research, 43(8), W08422. https://doi.org/10.1029/2006WR005398