Semi-Automated Sample-Based Forest Degradation Monitoring with Photointerpretation of High-Resolution Imagery

Forest fragmentation and degradation are a problem in many areas of the world and are a cause for concern to land managers. Similarly, countries interested in curtailing climate change have a keen interest in monitoring forest degradation. Traditional methods for measuring forested landscape pattern dynamics with maps made from classified satellite imagery fall short with respect to the compatibility of their forest definitions with information needs. In addition, they are not easily amenable to interpretation using tools like confidence intervals derived from survey sampling theory. In this paper, we described a novel landscape monitoring approach that helps fill these gaps. In it, a grid of photo plots is efficiently created and overlaid on high-resolution imagery, points are labeled with respect to their land-use by a human interpreter, and mean values and their variance are calculated for a suite of point-based fragmentation metrics related to forest degradation. We presented three case studies employing this approach from the US states of Maryland and Pennsylvania, highlighted different survey sampling paradigms, and discussed the strengths and weaknesses of the method relative to traditional, satellite imagery-based approaches. Results indicate that the scale of forest fragmentation in Maryland is between 250 and 1000 m, and this agrees with compatible estimates derived from raster analytical methods. There is a positive relationship between an index of housing construction and change in forest aggregation as measured by our metrics, and strong agreement between metric values collected by human interpretation of imagery and those obtained from a land cover map from the same period. We showed how the metrics respond to simulated degradation, and offered suggestions for practitioners interested in leveraging rapid photointerpretation for forest degradation monitoring.


Introduction
Links between forest fragmentation, forest loss, and forest degradation have drawn interest in recent years due to increased focus on carbon monitoring and climate change mitigation strategies. For example, the United Nations' Reducing Emissions from Deforestation and Degradation (REDD) program [1] requires the measurement and monitoring of forest cover in exchange for donor support for various forest conservation activities. Forest fragmentation has been proposed as one indicator of forest degradation [2], and thus its measurement and monitoring show promise for providing the information needed for REDD reporting [3,4]. For MD and PG, the sampling design we used was chosen as part of a larger study of land-use in Maryland [19]. We used 5046 plots (described below) for the MD study and 3465 plots for the PG study. Both samples were chosen in a spatially-balanced manner, using the sample selection procedure described in Lister and Scott [39]. This procedure uses a fractal (a space-filling curve) that covers the entire study area to create a tessellation that partitions the sampling frame into units from which samples are drawn in a spatially uniform, equal probability manner; this is analogous to a systematic design but without the same degree of regularity. For SC, the sampling design was chosen as part of a forest inventory precision improvement study by Westfall et al. [23]. We used 2099 plots (described below), again chosen in the manner described by Lister and Scott [39].

Plot Design
The plot design of the MD sample was chosen as part of a larger land-use change study [19]. For that study, a cluster design was chosen, with four points separated by 500 m arranged in a square pattern (Figure 2a). This configuration minimized cost while maximizing land-use change estimate precision.
For PG, a plot design selection procedure was implemented, similar in principle to the concept of cyclic sampling described in Burrows et al. [40]. Based on PI heuristics established by Lister et al. [19], it was determined that a cluster plot with an array of 10 subplots (points) would conservatively allow us to finish all plots within the time and budget allotted for the study. Since the fragmentation metric used in this study, described below, relies upon grouping points into separation distance classes and summarizing land-use class adjacency frequencies, our goal was to select a 10-point configuration that would provide a roughly equal count of pairs of points per separation distance class. This would allow for the calculation of the metric with similar numbers of points in each class. To achieve this, we implemented a procedure in which a random set of 10 points was drawn, with replacement, 50,000 times from a grid of 25 equidistant (50 m spacing) points distributed on a square, 200 m × 200 m grid (Figure 2b). For each set of points drawn, we calculated the separation distance between each point and every other point in the 10-point sample. Next, we counted the number of pairs of points in each separation distance class and stored the result. We then randomly chose the final plot design from For MD and PG, the sampling design we used was chosen as part of a larger study of land-use in Maryland [19]. We used 5046 plots (described below) for the MD study and 3465 plots for the PG study. Both samples were chosen in a spatially-balanced manner, using the sample selection procedure described in Lister and Scott [39]. This procedure uses a fractal (a space-filling curve) that covers the entire study area to create a tessellation that partitions the sampling frame into units from which samples are drawn in a spatially uniform, equal probability manner; this is analogous to a systematic design but without the same degree of regularity. For SC, the sampling design was chosen as part of a forest inventory precision improvement study by Westfall et al. [23]. We used 2099 plots (described below), again chosen in the manner described by Lister and Scott [39].

Plot Design
The plot design of the MD sample was chosen as part of a larger land-use change study [19]. For that study, a cluster design was chosen, with four points separated by 500 m arranged in a square pattern (Figure 2a). This configuration minimized cost while maximizing land-use change estimate precision.
For PG, a plot design selection procedure was implemented, similar in principle to the concept of cyclic sampling described in Burrows et al. [40]. Based on PI heuristics established by Lister et al. [19], it was determined that a cluster plot with an array of 10 subplots (points) would conservatively allow us to finish all plots within the time and budget allotted for the study. Since the fragmentation metric used in this study, described below, relies upon grouping points into separation distance classes and summarizing land-use class adjacency frequencies, our goal was to select a 10-point configuration that would provide a roughly equal count of pairs of points per separation distance class. This would allow for the calculation of the metric with similar numbers of points in each class. To achieve this, we implemented a procedure in which a random set of 10 points was drawn, with replacement, 50,000 times from a grid of 25 equidistant (50 m spacing) points distributed on a square, 200 m × 200 m grid (Figure 2b). For each set of points drawn, we calculated the separation distance between each point and every other point in the 10-point sample. Next, we counted the number of pairs of points in each separation distance class and stored the result. We then randomly chose the final plot design from among the handful of configurations that yielded a uniform count of points per distance class, which, in our case, was nine pairs per plot (Figure 2b). among the handful of configurations that yielded a uniform count of points per distance class, which, in our case, was nine pairs per plot (Figure 2b). For SC, the plot consisted of an array of 52 points arranged uniformly within a circular area of 5810 square meters. This area was chosen because a circle of this size circumscribes a single FIA cluster plot, which was the focus of the study from which the data came [23]. Interpoint spacing was approximately 12 m. The design for PG, which was chosen in order to achieve a roughly equal number of pairs of points per separation distance class. (c) The design for SC, which was chosen as part of a larger study on forest inventory precision.

Fragmentation and Degradation Indicators
For MD and PG, the forest fragmentation metric we chose to illustrate the point-based approach, what we call Point Aggregation Index (PAI) (Table 1), is based on the probabilities used to calculate Shannon and Weaver's [41] conditional entropy. The metric uses class adjacency probabilities as its foundation, like entropy-based metrics, such as Q, the nearest neighbor probability metric [42], Cluster/Interspersion (CL) index [43], and the similar Aggregation Index (AI) [44]. PAI can be used as a component of several contagion indices calculated from the entropy of attribute adjacency [45], only using points instead of raster cells as the element with which probabilities are calculated. The minimum value of PAI is 0 (there are no adjacent points of the class of interest), and the maximum value is 1 (the entire landscape is made up of one class). One key feature of calculating contagion indices based on PAI is that it can be done with a userspecified definition of adjacency. With pixel-based approaches, adjacency is typically limited to different definitions of pixel contiguity (e.g., the 4-or 8 neighbor rule). With points, adjacency can be defined using an annular neighborhood approach-points are considered adjacent if they are found within a chosen distance band of one another. Figure 3d illustrates the calculation of PAI on the scenario depicted in Figure 3a, only using a larger neighborhood to define adjacency. For SC, the plot consisted of an array of 52 points arranged uniformly within a circular area of 5810 square meters. This area was chosen because a circle of this size circumscribes a single FIA cluster plot, which was the focus of the study from which the data came [23]. Interpoint spacing was approximately 12 m.

Fragmentation and Degradation Indicators
For MD and PG, the forest fragmentation metric we chose to illustrate the point-based approach, what we call Point Aggregation Index (PAI) (Table 1), is based on the probabilities used to calculate Shannon and Weaver's [41] conditional entropy. The metric uses class adjacency probabilities as its foundation, like entropy-based metrics, such as Q, the nearest neighbor probability metric [42], Cluster/Interspersion (CL) index [43], and the similar Aggregation Index (AI) [44]. PAI can be used as a component of several contagion indices calculated from the entropy of attribute adjacency [45], only using points instead of raster cells as the element with which probabilities are calculated. The minimum value of PAI is 0 (there are no adjacent points of the class of interest), and the maximum value is 1 (the entire landscape is made up of one class). One key feature of calculating contagion indices based on PAI is that it can be done with a user-specified definition of adjacency. With pixel-based approaches, adjacency is typically limited to different definitions of pixel contiguity (e.g., the 4-or 8 neighbor rule). With points, adjacency can be defined using an annular neighborhood approach-points are considered adjacent if they are found within a chosen distance band of one another. Figure 3d illustrates the calculation of PAI on the scenario depicted in Figure 3a, only using a larger neighborhood to define adjacency.
For SC, several adjacency metrics were adapted for use with points or developed with the goal of creating index values that can be used as plot-based values to use in survey sampling estimation procedures. Furthermore, at the recommendation of Li et al. [46], we focused on commonsense metrics that are related to biological functions forest provide. A list of these metrics is provided in Table 1.  Table 1. For SC, several adjacency metrics were adapted for use with points or developed with the goal of creating index values that can be used as plot-based values to use in survey sampling estimation procedures. Furthermore, at the recommendation of Li et al. [46], we focused on commonsense metrics that are related to biological functions forest provide. A list of these metrics is provided in Table 1.

PI and Subpopulation-Level Estimation Procedures
The land-use category was assessed on each point on every plot by interpreting digital aerial imagery. The 1998 imagery (used as the time 1 dataset for the MD study) consisted of the panchromatic, leaf-on, 2-m-pixel resolution, digital orthophoto quadrangles (DOQs) from a state-level imagery dataset stored locally in an ArcGIS [47] raster catalog. The 2007 imagery (used as the time 2 MD data and for the PG study) consisted of the color infrared, leaf-on, 1-m-pixel resolution, digital images from the National Aerial Imagery Program (NAIP) served over the Internet using a web-mapping service (WMS). The definitions of the forest and nonforest land-use categories used were those used by FIA [8]. A single interpreter was trained and conducted all photo analyses for this study. Quality assurance protocols included a second assessment of a portion of the points by another interpreter and the requirement that the first interpreter re-interpret several hundred randomly-chosen points to assess consistency. In both cases, at least 95% of the interpreter's work met quality standards (i.e., was correct and/or consistent).
To increase PI efficiency, an automation method was developed whereby an ArcGIS tool was used to create subsets of imagery from the WMS to areas encompassing and slightly beyond the extent of the footprint of each plot. In other words, "snapshots" of imagery at a scale of 1:4000 were generated, with each image centered on the plot and containing sufficient detail for the interpreter to assess landuse status. The sets of images for both studies were stored locally and displayed using a Microsoft Access form developed for the project (Figure 4a). A template of this form with all functionality is available in the supplementary material. The form was designed to display the images and allow for data entry in such a way that the number of mouse clicks, wait time for image loading, and data entry times were minimized. In addition, links to both "Google Maps" and "Bing Maps" were enabled for each point so that contextual information, and, where possible, Google's "Street View", could be used to guide the interpreter.

PI and Subpopulation-Level Estimation Procedures
The land-use category was assessed on each point on every plot by interpreting digital aerial imagery. The 1998 imagery (used as the time 1 dataset for the MD study) consisted of the panchromatic, leaf-on, 2-m-pixel resolution, digital orthophoto quadrangles (DOQs) from a state-level imagery dataset stored locally in an ArcGIS [47] raster catalog. The 2007 imagery (used as the time 2 MD data and for the PG study) consisted of the color infrared, leaf-on, 1-m-pixel resolution, digital images from the National Aerial Imagery Program (NAIP) served over the Internet using a web-mapping service (WMS). The definitions of the forest and nonforest land-use categories used were those used by FIA [8]. A single interpreter was trained and conducted all photo analyses for this study. Quality assurance protocols included a second assessment of a portion of the points by another interpreter and the requirement that the first interpreter re-interpret several hundred randomly-chosen points to assess consistency. In both cases, at least 95% of the interpreter's work met quality standards (i.e., was correct and/or consistent).
To increase PI efficiency, an automation method was developed whereby an ArcGIS tool was used to create subsets of imagery from the WMS to areas encompassing and slightly beyond the extent of the footprint of each plot. In other words, "snapshots" of imagery at a scale of 1:4000 were generated, with each image centered on the plot and containing sufficient detail for the interpreter to assess land-use status. The sets of images for both studies were stored locally and displayed using a Microsoft Access form developed for the project (Figure 4a). A template of this form with all functionality is available in the supplementary material. The form was designed to display the images and allow for data entry in such a way that the number of mouse clicks, wait time for image loading, and data entry times were minimized. In addition, links to both "Google Maps" and "Bing Maps" were enabled for each point so that contextual information, and, where possible, Google's "Street View", could be used to guide the interpreter.
For the MD study area, we calculated PAI for five distance classes with cut-offs of 500, 1000, 2000, 3000, 4000, and 5000 m for each of two periods (1998 and 2007), and aggregated values by subpopulations consisting of county groups we created using county-level densities of new housing permits issued between 1990 and 2010 (# housing unit permits issued/area of county) [50]. The county groups, which we hypothesized would correspond with different levels of change in PAI, were formed by implementing Jenks' optimization method [51], which iteratively chooses class breaks, seeking to minimize within-group and maximize the between-group variance of housing permit density ( Figure 5). The goal of these efforts was to characterize the relationships between PAI and scale of analysis (point separation distance) and use them to assess fragmentation at different scales. Table 1. Landscape metrics calculated at the plot-level for the forest fragmentation and degradation analysis conducted in the Schuylkill County (SC) study area. For Maryland (MD) and Prince George's County (PG), only PAI (Point Aggregation Index) was calculated, summarized at the landscape-level. The subset of plots used represents the subset of the total number of plots in the population for which calculation of the metric is feasible. Infeasible calculations can consist of scenarios leading to division by zero, or, for some metrics, distance calculations between forested points on plots with 1 or fewer forested points. Sample sizes (n) used for each metric in the SC study area are given.

Subset of Plots Used Equation Description
Point Aggregation Index (PAI) 1 Plots with at least 1 forested point The number of forest-forest point adjacencies divided by the total number of adjacencies between forest and any other type (including forest) in the estimation area. Calculated with A ij , an i × j (i = j = 1 . . . m points in the estimation area) adjacency matrix with binary elements (1= a forested point adjacent to another forested point, 0 otherwise) and B ij , an i × j adjacency matrix with binary elements (1 = a forested point adjacent to any point type, 0 otherwise). Values can range from 0 (the case where there are no adjacent forest points) to 1 (the case where all points are forest).

TSP distance (TSPd)
Plots with at least 2 forested points (n = 1599) TSPd = i∈D d i The shortest path distance (meters) calculated using an algorithmic solution to what is known as the Traveling Salesperson Problem (TSP), farthest insertion method [48] found within the R package TSP [49]. d i = the i th member of the set D of interpoint distances leading to the shortest overall distance between all forested points on the plot, returning to the first point. Values range from the minimum distance between two adjacent points (the case where there are 2 adjacent forested points on the plot; 12 m in the case of SC) to the sum of the shortest path roundtrip distance between all points on the plot (the case where all points on the plot are forest; approximately 675 m in the case of SC).

Relative TSP distance (rTSPd)
Plots with at least 2 forested points (n = 1599) rTSPd = TSPd k Relative (average per point) distance along the shortest path between all forested points on the plot, returning to point 1. k = the number of forested points in the plot. Values range from 1 2 the distance between 2 adjacent points (the case where there are k = 2 adjacent forested points on the plot) to 1/k th of the sum of the shortest path roundtrip distance between all points on the plot (the case where all points on the plot are forest; in the case of SC, k = 52).

Nearest Neighbour distance (NNd)
Plots with at least 2 forested points (n = 1599) NNd = min C ij, j i Nearest neighbor distance, calculated as the nonzero minimum of the i × j (i = j = 1 . . . k forested points on the plot) distance matrix C with elements comprised of geographic distance between forested points i and j. Values range from the minimum distance between two adjacent points (the case where there are k = 2 adjacent forested points on the plot; 12 m in the case of SC) to the maximum distance between 2 points on the plot (the case where there are k = 2 forest points located at the extreme ends of the plot; approximately 82 m in the case of SC).

Mean Interpoint Distance (MId)
Plots with at least 2 forested points (n = 1599) The average distance between forested points on the plot, calculated as the mean of the i × j (i = j = 1 . . . k forested points on the plot) distance matrix D with elements comprised of geographic distance between forested points i and j. Values range from 1 2 the distance between 2 adjacent points (the case where there are 2 adjacent forested points on the plot) to 1 2 the distance between the 2 points located at extreme ends of the plot (the case where there are 2 forested points located at the extreme ends of the plot).

Number of f-f point adjacencies (NA)
Plots with at least 1 forested point (n = 1602) The number of forest-forest point adjacencies on the plot, summarized from A ij , the i × j (i = j = 1 . . . m points in the plot) adjacency matrix with binary elements (1= a forested point adjacent to another forested point, 0 otherwise). Values range from 0 (the case where there are no forest-forest adjacencies) to the total count of possible point to point adjacencies on the plot (the case where all points are forest; 170 m in the case of SC).

Relative number of adjacencies (rNA)
Plots with at least 1 forested point The number of forest-forest point adjacencies divided by the total number of possible adjacencies of any type on the plot, calculated with NA and D ij , the i × j (i = j = 1 . . . m points in the plot) adjacency matrix obtained when all points are a forest. Values range between 0 (the case where there are no forest-forest adjacencies) and 1 (the case where all points are forest). 1 For MD and PG, the estimation area referred to in the description is the analysis unit as a whole; for SC, the estimation area is individual plots.  ii. An array of points comprising the cluster plot. Points touching forest canopy are "painted" by holding the shift key while passing the mouse over them; iii. Record navigation tools to facilitate switching between plots and selected sets of points.

Mean Forest
For the MD study area, we calculated PAI for five distance classes with cut-offs of 500, 1000, 2000, 3000, 4000, and 5000 m for each of two periods (1998 and 2007), and aggregated values by subpopulations consisting of county groups we created using county-level densities of new housing permits issued between 1990 and 2010 (# housing unit permits issued / area of county) [50]. The county groups, which we hypothesized would correspond with different levels of change in PAI, were formed by implementing Jenks' optimization method [51], which iteratively chooses class breaks, seeking to minimize within-group and maximize the between-group variance of housing permit density ( Figure 5). The goal of these efforts was to characterize the relationships between PAI and scale of analysis (point separation distance) and use them to assess fragmentation at different scales.   ii. An array of points comprising the cluster plot. Points touching forest canopy are "painted" by holding the shift key while passing the mouse over them; iii. Record navigation tools to facilitate switching between plots and selected sets of points.
For the MD study area, we calculated PAI for five distance classes with cut-offs of 500, 1000, 2000, 3000, 4000, and 5000 m for each of two periods (1998 and 2007), and aggregated values by subpopulations consisting of county groups we created using county-level densities of new housing permits issued between 1990 and 2010 (# housing unit permits issued / area of county) [50]. The county groups, which we hypothesized would correspond with different levels of change in PAI, were formed by implementing Jenks' optimization method [51], which iteratively chooses class breaks, seeking to minimize within-group and maximize the between-group variance of housing permit density ( Figure 5). The goal of these efforts was to characterize the relationships between PAI and scale of analysis (point separation distance) and use them to assess fragmentation at different scales.

Variance Estimation for PAI for the MD and PG Study Areas
PAI, as we calculated it for MD and PG, is an index calculated from the set of between-subplot and between-plot values collected on a subpopulation of our PI plots. Because the nature of the sampling distribution of this statistic is unknown, and for illustrative purposes, we chose to use a delete-one jackknife variance estimation approach, described in Thompson [52] and used by Kleinn [34] and Ramezani and Ramezani [38]. In this approach, one cluster plot is deleted from the sample, and PAI is calculated. This is repeated for each plot in succession, and the estimator of the variance of PAI v PAI is calculated as follows: where n = the number of cluster plots used, PAI (i) is PAI generated from the ith deletion, and PAI (.) is the mean of the set of PAI's calculated with n − 1 plots each. An approximate 100(1 − α)% confidence interval for PAI is calculated as where t is the value from the t distribution with n − 1 degrees of freedom.

Comparison with Traditional Raster Data and Methods for the PG Study Area
In order to compare our results with more traditional methods for the PG study area, we used the same sample and plot designs to calculate PAI by using a GIS to intersect the point locations with a raster forest/nonforest map derived from a reclassification of the NLCD 2001 land cover product. The reclassification consisted of grouping all NLCD classes in the 40 range (forest) and 91 (woody wetland) into a forest (f ) category and everything else into the nonforest (nf ) category. In addition to calculating PAI using this f/nf dataset, we also calculated average patch size by first assigning each pixel to a patch based on its contiguity with other pixels of the same class (using an 8 neighbor rule), and then averaging the areas of the patches that were formed.
As an improvement to the traditional, raster-based patch size calculation, we also developed an alternative measure of a property for which patch size is often used as a surrogate-a variable we called mean traversal width-length index (TWL). Mean TWL statistics are calculated by summarizing the area of the class of interest on a per-patch or analysis unit basis using a set of perpendicular transects established through the landscape ( Figure 6): where L i is the length of a transect found within patch i of n patches found on the landscape. The intent is to characterize that property of a patch or set of patches relates to unidirectional flow, as occurs with many materials (wildlife, propagules, wind, water, etc.) that traverse an ecosystem.
We calculated TWL by first establishing what amounts to sets of 30-m wide northerly-and southerly-transects through the 30-m resolution f/nf image by assigning to each pixel the value of its Albers NAD83 Easting and Northing. We then attached to each pixel the f or nf value from the f/nf map and used the ArcGIS 10.0 [47] command "regiongroup" to attach a unique region identifier to each segment of contiguous f pixels formed by both the Easting and Northing transects. This created a set of uniquely-identified, 30-m wide strips of contiguous forest pixels in both the north-south and the east-west directions. The length of each strip (area/30) was then obtained, and a histogram of TWL values was generated.
We acknowledge that TWL is not rotation invariant, but we assumed that the two perpendicular directions would give an adequate characterization of TWL because we saw no evidence of anisotropic pattern on the landscape; if a pattern of north-south or east-west orientation of features is suspected, a random rotation of the land cover image can be easily applied in a GIS prior to TWL generation. each segment of contiguous f pixels formed by both the Easting and Northing transects. This created a set of uniquely-identified, 30-m wide strips of contiguous forest pixels in both the north-south and the east-west directions. The length of each strip (area/30) was then obtained, and a histogram of TWL values was generated.  (3)).
We acknowledge that TWL is not rotation invariant, but we assumed that the two perpendicular directions would give an adequate characterization of TWL because we saw no evidence of anisotropic pattern on the landscape; if a pattern of north-south or east-west orientation of features is suspected, a random rotation of the land cover image can be easily applied in a GIS prior to TWL generation.

Data Collection Methods, Metric Calculation, and Estimation Procedures for the SC Study Area
The land-use category was assessed on each point of each plot by interpreting 1-m-pixel resolution color NAIP imagery from 2015, obtained from a WMS [15]. The definitions of the forest and nonforest categories used were the same as those for the MD and PG study areas. A single interpreter was trained and conducted all photo analyses for this study using the same image preparation procedure and the Microsoft Access form approach described above (Figure 4b). The dataset used and the form are available in the supplementary material. Quality assurance protocols were similar to those described previously for MD and PG.
The metrics shown in Table 1 were calculated using code written in the R Statistical Software language [53], provided in the supplementary material. For SC, the paradigm of estimate calculation is plot-based; i.e., metrics for individual plots were calculated and, thus, treated as plot-level random variables. In order to demonstrate how to estimate change and forest degradation using this set of random variables and estimators from survey sampling, a random number drawn from a uniform distribution ranging from 0 to 1 was assigned to each point where forest cover was detected. Points with a random number value between 0.99 and 1 were removed in order to simulate a random, uniform degradation (loss) of 1% of the forest across the study area; hereafter, we have referred to this simulated degradation dataset as the Time 2 dataset. Metrics in Table 1 were again calculated with the Time 2 data, and estimates of change and their variance were generated. For each estimate, including plot-level change, the following estimators were applied for the (sub)population mean ( ) and its variance ( ( )), assuming simple random sampling:  (3)).

Data Collection Methods, Metric Calculation, and Estimation Procedures for the SC Study Area
The land-use category was assessed on each point of each plot by interpreting 1-m-pixel resolution color NAIP imagery from 2015, obtained from a WMS [15]. The definitions of the forest and nonforest categories used were the same as those for the MD and PG study areas. A single interpreter was trained and conducted all photo analyses for this study using the same image preparation procedure and the Microsoft Access form approach described above (Figure 4b). The dataset used and the form are available in the supplementary material. Quality assurance protocols were similar to those described previously for MD and PG.
The metrics shown in Table 1 were calculated using code written in the R Statistical Software language [53], provided in the supplementary material. For SC, the paradigm of estimate calculation is plot-based; i.e., metrics for individual plots were calculated and, thus, treated as plot-level random variables. In order to demonstrate how to estimate change and forest degradation using this set of random variables and estimators from survey sampling, a random number drawn from a uniform distribution ranging from 0 to 1 was assigned to each point where forest cover was detected. Points with a random number value between 0.99 and 1 were removed in order to simulate a random, uniform degradation (loss) of 1% of the forest across the study area; hereafter, we have referred to this simulated degradation dataset as the Time 2 dataset. Metrics in Table 1 were again calculated with the Time 2 data, and estimates of change and their variance were generated. For each estimate, including plot-level change, the following estimators were applied for the (sub)population mean (Y) and its variance (v Y ), assuming simple random sampling: where y i = the attribute of interest on plot i expressed as the value of the index calculated as described in Table 1; for change estimation, the attribute on each plot is calculated as I T2 -I T1 , where I is the calculated index from Table 1, and T1 and T2 are time periods 1 and 2, respectively, and n = number of plots used in the estimation. For n used in our study, see Table 1. For change estimation, the subset of plots that have per-plot index values calculated at both periods is used. The variance of the (sub)population mean is estimated as: where y i and n are defined as before and s 2 is the sample variance, ignoring the finite population correction factor when sample sizes are sufficiently large.

PAI and NLCD Metrics in the PG Study Area
The plot of PAI vs. separation distance class for both the PI and the NLCD f/nf map is shown in Figure 7. As expected, PAI varied at different scales of observation, or in our parlance, for different definitions of adjacency. The fact that PAI from both PI and NLCD track each other so closely supports our assumption that our human interpreter and PI method is compatible with automated methods and maps derived from image classification. calculated index from Table 1, and T1 and T2 are time periods 1 and 2, respectively, and n = number of plots used in the estimation. For n used in our study, see Table 1. For change estimation, the subset of plots that have per-plot index values calculated at both periods is used.
The variance of the (sub)population mean is estimated as: where yi and n are defined as before and s 2 is the sample variance, ignoring the finite population correction factor when sample sizes are sufficiently large.

PAI and NLCD Metrics in the PG Study Area
The plot of PAI vs. separation distance class for both the PI and the NLCD f/nf map is shown in Figure 7. As expected, PAI varied at different scales of observation, or in our parlance, for different definitions of adjacency. The fact that PAI from both PI and NLCD track each other so closely supports our assumption that our human interpreter and PI method is compatible with automated methods and maps derived from image classification. The pattern of the decay of PAI with distance was expected. The NLCD 2001 map showing the f/nf pattern in Prince George's County indicates a mosaic of irregularly shaped and sized forest patches (e.g., Figure 6). This led us to expect a broad range of values, with smaller separation distances yielding a higher PAI (the existence of larger, aggregated patches would lead to this), and larger separation distances yielding lower values (separation distances beyond the average patch width, among other things, would cause this). The phenomenon that leads to this pattern of decay of PAI with separation distance is related to spatial autocorrelation, a condition whereby phenomena located closer together are more similar than those located farther apart. Spatial autocorrelation is a characteristic of landscape pattern, and its causes and effects are well-described in the literature [54]. The pattern of the decay of PAI with distance was expected. The NLCD 2001 map showing the f/nf pattern in Prince George's County indicates a mosaic of irregularly shaped and sized forest patches (e.g., Figure 6). This led us to expect a broad range of values, with smaller separation distances yielding a higher PAI (the existence of larger, aggregated patches would lead to this), and larger separation distances yielding lower values (separation distances beyond the average patch width, among other things, would cause this). The phenomenon that leads to this pattern of decay of PAI with separation distance is related to spatial autocorrelation, a condition whereby phenomena located closer together are more similar than those located farther apart. Spatial autocorrelation is a characteristic of landscape pattern, and its causes and effects are well-described in the literature [54].
Monitoring changes in the shape of the PAI-separation distance class relationship can offer information relevant to landscape-scale forest dynamics. The form of the relationship between PAI and separation distance indicates that beyond the separation distance range of approximately 250-1000 m, there is a relatively small decrease in the likelihood that two points chosen at random will both be a forest. In other words, the effect of local pattern on PAI diminishes, and the effect of the landscape-scale pattern has begun to take over. This hypothesis is supported by our analysis of the NLCD dataset. The average patch area was 21.2 ha, or, if a square patch is assumed, 460 m on a side (520 m diameter if a circular patch is assumed). The frequency distribution ( Figure 8) and associated summary statistics of the components of the TWL metric, which characterize a property of the landscape related to the linear flow of materials, indicate a relatively similar size range: average TWL was 240 m, and median TWL was 120 m.
Monitoring changes in the shape of the PAI-separation distance class relationship can offer information relevant to landscape-scale forest dynamics. The form of the relationship between PAI and separation distance indicates that beyond the separation distance range of approximately 250-1000 m, there is a relatively small decrease in the likelihood that two points chosen at random will both be a forest. In other words, the effect of local pattern on PAI diminishes, and the effect of the landscape-scale pattern has begun to take over. This hypothesis is supported by our analysis of the NLCD dataset. The average patch area was 21.2 ha, or, if a square patch is assumed, 460 m on a side (520 m diameter if a circular patch is assumed). The frequency distribution ( Figure 8) and associated summary statistics of the components of the TWL metric, which characterize a property of the landscape related to the linear flow of materials, indicate a relatively similar size range: average TWL was 240 m, and median TWL was 120 m.

Temporal Differences in PAI across Different Scales and Historical Housing Permit Densities in the MD Study Area
The relationship between historical housing permit density ( Figure 5), separation distance, and change in PAI over time can be seen in Figure 9. A pattern of decay of PAI with separation distance, similar to that seen in the PG study area, occurs at the state level. Almost all of the distance classes showed a significant (non-overlapping 95% confidence intervals) difference in PAI values between housing unit density levels. For the medium and high housing unit density levels, there was a consistent trend of reductions in PAI between 1998 and 2007 for all distance classes. These findings align with our assumption that PAI differences would be correlated with processes that affect landscape pattern status and dynamics; housing unit density permits can be thought of as a surrogate for development activity, leading to forest degradation and loss. The lower PAI values and greater reduction in PAI over time are found in the counties in the more populous, less-forested portions of the state-those affected most by development during the decade under study ( Figure 5).
It is important to note that for PG and MD, PAI results were calculated as a ratio of counts of point adjacencies at the estimation unit level, not first for individual plots, and then summarized for estimation units. Under this paradigm, our estimates of PAI and the precisions of the estimates are interpreted as characteristics of a subpopulation of unknown size within the analysis unit. The subpopulation must consist of the set of population elements with at least one forested point on them; otherwise, PAI calculation would be impossible because there would be a zero in the denominator (Table 1, PAI). Because of this and the fact that plots contribute unequally to the subpopulation-level, ratio-based estimate, the form of the sampling distribution of PAI is difficult to know. The jackknife

Temporal Differences in PAI across Different Scales and Historical Housing Permit Densities in the MD Study Area
The relationship between historical housing permit density ( Figure 5), separation distance, and change in PAI over time can be seen in Figure 9. A pattern of decay of PAI with separation distance, similar to that seen in the PG study area, occurs at the state level. Almost all of the distance classes showed a significant (non-overlapping 95% confidence intervals) difference in PAI values between housing unit density levels. For the medium and high housing unit density levels, there was a consistent trend of reductions in PAI between 1998 and 2007 for all distance classes. These findings align with our assumption that PAI differences would be correlated with processes that affect landscape pattern status and dynamics; housing unit density permits can be thought of as a surrogate for development activity, leading to forest degradation and loss. The lower PAI values and greater reduction in PAI over time are found in the counties in the more populous, less-forested portions of the state-those affected most by development during the decade under study ( Figure 5).
It is important to note that for PG and MD, PAI results were calculated as a ratio of counts of point adjacencies at the estimation unit level, not first for individual plots, and then summarized for estimation units. Under this paradigm, our estimates of PAI and the precisions of the estimates are interpreted as characteristics of a subpopulation of unknown size within the analysis unit. The subpopulation must consist of the set of population elements with at least one forested point on them; otherwise, PAI calculation would be impossible because there would be a zero in the denominator (Table 1, PAI). Because of this and the fact that plots contribute unequally to the subpopulation-level, ratio-based estimate, the form of the sampling distribution of PAI is difficult to know. The jackknife variance estimator, therefore, appears to be a reasonable approach to addressing this situation, as it does not require distributional assumptions. Both Kleinn [34] and Ramezani and Ramezani [38] used the jackknife variance estimator when calculating landscape metrics based on forest inventory plots.
It must be noted, however, that procedures for comparing means using this method are not straightforward. For example, we inferred significant differences between county-group PAI values for all distance classes because confidence intervals do not overlap. However, inferring the significance of differences in means using confidence interval overlap is not strictly correct as it assumes additivity of standard errors [55]. In our case, since the confidence interval for the difference is inflated due to this unfounded assumption, our conclusion is justified, but clearly, this is not a satisfactory approach to use in all cases, and better ways to compare means with jackknife variance estimators are needed.
Forests 2019, 10, x FOR PEER REVIEW  14 of 19 variance estimator, therefore, appears to be a reasonable approach to addressing this situation, as it does not require distributional assumptions. Both Kleinn [34] and Ramezani and Ramezani [38] used the jackknife variance estimator when calculating landscape metrics based on forest inventory plots. It must be noted, however, that procedures for comparing means using this method are not straightforward. For example, we inferred significant differences between county-group PAI values for all distance classes because confidence intervals do not overlap. However, inferring the significance of differences in means using confidence interval overlap is not strictly correct as it assumes additivity of standard errors [55]. In our case, since the confidence interval for the difference is inflated due to this unfounded assumption, our conclusion is justified, but clearly, this is not a satisfactory approach to use in all cases, and better ways to compare means with jackknife variance estimators are needed.

Sample-Based Estimates of Metrics and Change in the SC Study Area
In contrast to the "single observation unit" paradigm employed in the PG and MD study areas, we chose a traditional survey sampling paradigm for the SC study area. In other words, we performed estimation by first aggregating adjacency counts and other information to the plot level, calculating the metrics in Table 1 by plot, and then calculating means and their variances using this set of plot-level metric values. This does two things-it allows us to use traditional estimators of mean and variance of means, and it shifts the focus to the subpopulation of existing forest and the plot-scale, which in the case of SC is tied to our plot geometry, approximately 1-ha. Under this paradigm, the subpopulation sample size is a random variable, but Thompson [52] (p. 63) showed that for sufficiently large sample sizes, an unbiased estimator of the conditional variance, given the subpopulation sample size, is the same as that shown in equation 5. For forest degradation, estimates of changes in ~hectare-scale metrics calculated only within a specific subpopulation of forested points are arguably more relevant than those calculated at the landscape scale with these metrics.
For the metrics Mean Forest Proportion 18-m window (MFP18), SD Forest Proportion 18-m window (SDFP18), and Proportion Forest (PF), which are calculated using all plots and points, there is no distinction between the estimation paradigms, as parameter estimates are interpreted the same Figure 9. The relationship between three counties grouped by historical housing permit density levels (low, medium, and high), imagery date (1998 and 2007), and the fragmentation metric PAI across different separation distance classes.

Sample-Based Estimates of Metrics and Change in the SC Study Area
In contrast to the "single observation unit" paradigm employed in the PG and MD study areas, we chose a traditional survey sampling paradigm for the SC study area. In other words, we performed estimation by first aggregating adjacency counts and other information to the plot level, calculating the metrics in Table 1 by plot, and then calculating means and their variances using this set of plot-level metric values. This does two things-it allows us to use traditional estimators of mean and variance of means, and it shifts the focus to the subpopulation of existing forest and the plot-scale, which in the case of SC is tied to our plot geometry, approximately 1-ha. Under this paradigm, the subpopulation sample size is a random variable, but Thompson [52] (p. 63) showed that for sufficiently large sample sizes, an unbiased estimator of the conditional variance, given the subpopulation sample size, is the same as that shown in equation 5. For forest degradation, estimates of changes in~hectare-scale metrics calculated only within a specific subpopulation of forested points are arguably more relevant than those calculated at the landscape scale with these metrics.
For the metrics Mean Forest Proportion 18-m window (MFP 18 ), SD Forest Proportion 18-m window (SDFP 18 ), and Proportion Forest (PF), which are calculated using all plots and points, there is no distinction between the estimation paradigms, as parameter estimates are interpreted the same (all plots have equal impact and selection probability, and the population is the entire area of SC). However, with metrics that rely on either one or two forested points per plot, the subpopulation of interest is defined in the somewhat ambiguous way as "the subset of population elements having at least one (or two for some metrics) forest points within a circle of 82 m diameter". The plot-level values, in this case, are considered equal-impact measurements of some local characteristics of forest-related to its spatial configuration. Changes in plot-level mean values of these metrics can thus be seen as a way to track degradation in a manner that more closely addresses the concerns about changes in ecological function that originally motivated interest in degradation monitoring, at the appropriate scale. Table 2 contains results obtained under this paradigm for the PI analysis of the SC study area, alongside results obtained from a degradation simulation that represents 1% forest loss, randomly assigned to the set of all initially-forested points. Sample-based estimates of metric differences (degraded-original), with their standard errors and 95% confidence intervals, are also provided. The trajectory of change in metric values was as expected; when simulated degradation was applied, metrics related to adjacency counts (PAI, NA (Number of forest-forest point adjacencies), and rNA (Relative number of adjacencies)) decreased. This suggests that in the average~1-ha area containing forest, forested points became more dispersed. Table 2. Plot-based mean and precision values of metrics and change described in Table 1 and Equations (4) and (5), calculated for the SC study area. Entries consist of means, absolute and relative (%) standard errors (SE), and, for change, the 95% confidence interval half-widths (SE*1.96). Estimates are from the original photointerpretation (PI) dataset (Original); the same plots, but with 1% random forest loss applied (Degraded); and the mean difference (Degraded − Original). PF, which is not a fragmentation metric included in Similar conclusions can be drawn from changes in metrics related to the intra-plot distance between forested points (traveling salesperson problem distance (TSPd), relative TSP distance (rTSPd), nearest neighbor distance (NNd), and mean interpoint distance (MId)). It is noteworthy that while mean TSPd (the total distance required to visit each forested point on the plot and return to the first point) decreased, the relative value (i.e., the per-point total distance) increased, again suggesting that functions of forest that rely on networks of interconnected forest points at the~1-ha scale, such as providing habitat to forest-dwelling wildlife species, have become degraded. Finally, MFP 18 , SDFP 18 , and PF, which were calculated with all points in SC, show that within the land area of SC, on average, there is not only less forest overall, but points in this area have less forest immediately surrounding them, and the proportion of forest surrounding them is more variable. All differences in means given in Table 2 (except for NNd) are significantly different from zero, given the large sample sizes; this, however, is not relevant if the magnitudes of the differences do not have a practical impact.

The Validity of Fragmentation Metrics and Functional Significance
We would like to emphasize that none of the metrics or change estimators we used could effectively measure what we term functional degradation or the level of forest fragmentation or degradation that has ecological meaning [56,57]. Similarly, many landscape configurations could yield the same value for each of these metrics. Since a full comparison of the performance of our metrics under different scenarios is beyond the scope of the current research, we refer the reader to excellent reviews of this topic [34,46,58,59], and will suffice it to say that in general, the metrics should be seen as components of some multidimensional property of a landscape related to forest quantity, distribution, and patch shape, not necessarily as physical properties with a direct biological meaning. In other words, these types of analyses should serve as guides that will alert managers to the trajectory of forest degradation.

Costs
With respect to costs, our method allowed for the interpretation of roughly one plot (10 points) every 13 seconds for PG, one plot (4 points) every 6 seconds for MD, and one plot (52 points) every 45 seconds for SC. Knowing the time and cost per sampling unit, and the expected frequency of the attribute of interest on the landscape, project managers can very precisely calculate the budgets for PI-based land cover assessments with an expected level of precision [19]. They can also easily conduct quality assurance, and when change is detected, examine and document the details of the change from the imagery. Forest monitoring groups that are required to characterize forest status or degradation in a detailed way with a limited budget can benefit from the PI-based monitoring approach we present. Future work will involve refining the set of landscape metrics presented in Table 1 and assessing their performance under different heterogeneity and change scenarios.

Applications
This research was motivated by the need for information on forest degradation using definitions that correspond precisely with what a human interpreter would choose, as well as interest in leveraging the abundance of georeferenced national forest inventory data [34,38,60] and PI-based forest monitoring information [28] that has become readily available. Potential applications include multi-scale comparisons of changes in metric values or differences between analysis units (e.g., county groups, Figure 5, and related discussion), or use as a covariate in assessments of forest condition and health, as has been done with other fragmentation metrics [61,62]. Finally, results can be used by planners and land managers as a component of a forest monitoring system like that required for participation in REDD+ to help prepare resource assessments and management plans.

Conclusions
We conclude that our approach to using data collected from PI for monitoring landscape pattern status and trends, and inferring degradation, offers many benefits. Results are similar to those obtained from summaries of raster maps, and we showed that the values of the metrics are correlated with the intensity of processes that lead to degradation. Our case studies illustrate two workflows and sampling paradigms and highlight the pros and cons of each. Overall, we find that leveraging existing PI data, like those being collected with tools like Collect Earth Online and other rapid PI methods, has great potential for filling information needs that are not currently met by more traditional degradation monitoring approaches relying on satellite imagery classification.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/10/896/s1: Input datasets for SC study area, R code for generating results tables for SC study area, and form for collecting data in SC with 20 images as an example.