Synthesizing Remote Sensing and Biophysical Measures to Evaluate Human–wildlife Conﬂicts: The Case of Wild Boar Crop Raiding in Rural China

: Crop raiding by wild boars is a growing problem worldwide with potentially damaging consequences for rural dwellers’ cooperation with conservation policies. Still, limited resources inhibit continuous monitoring, and there is uncertainty about the relationship between the biophysical realities of crop raiding and humans’ perceptions and responses. By integrating data from camera traps, remote sensors, and household surveys, this study establishes an empirical model of wild boar population density that can be applied to multiple years to estimate changes in distribution over time. It also correlates historical estimates of boar population distribution with human-reported trends to support the model’s validity and assess local perceptions of crop raiding. Although the model proved useful in coniferous and bamboo forests, it is less useful in mixed broadleaf, evergreen broadleaf, and deciduous forests. Results also show alignment between perceptions of crop raiding and actual boar populations, corroborating farmers’ perceptions which are increasingly dismissed as a less reliable source of information in human–wildlife conﬂict research. The modeling techniques demonstrated here may provide conservation practitioners with a cost-e ﬀ ective way to maintain up-to-date estimates of the spatial distribution of wild boar and resultant crop raiding.


Introduction
Human-wildlife conflict is an increasingly prevalent challenge throughout the world, with potentially severe implications for environmental conservation [1]. These conflicts are centered on competing social and environmental values, and while research cannot point to "correct" solutions, it can help conservation managers predict potential conflicts prior to policy implementation and manage them accordingly. Human-wildlife conflict is a special concern in protected areas and their surroundings. As ecosystems recover, growing populations of some wildlife species can threaten the livelihoods and safety of nearby rural communities [2], especially when these species invade farm fields and devour crops. Not only does crop raiding bring economic hardship to already poor families [3] and contribute to food shortages [4], it can also breed resistance to conservation programs [5] and interfere with their outcomes [3]. Managing wildlife crop raiding is thus socially and environmentally vital. While there have been several studies on crop raiding near protected areas [5][6][7][8], few have examined how crop raiding changes over time, likely because field estimates are expensive to obtain [9].
Crucial to the management of human-wildlife conflict is understanding the habitat and distribution of the wildlife species in question. Remotely sensed imagery, especially those from the Landsat satellite series, have long been used to estimate wildlife distribution and habitat suitability due to their cost-free socioeconomic effects of successful forest conservation programs, help with the development of appropriate mitigation and compensation programs, and further enhance the long-term success of conservation actions.

Study Site
Fanjingshan National Nature Reserve (N 27°44′42′′-28°03′11′′, W 108°34′19′′-108°48′30′′) is located in northeast ern Guizhou Province, southwestern China. Since its establishment in 1978, the reserve has attracted attention from conservationists globally and is inside one of the world's 25 biodiversity hotspots [32]. The mountainous 419-km 2 site is dominated by evergreen, deciduous, and mixed broadleaf forests with tracts of bamboo. Elevation varies widely from 400 m to 2560 m. The site is also home to about 13,000 people, most of whom are subsistence farmers, although some have migrated to cities for work or found employment in the area's burgeoning tourism sector [33]. Human settlements in the study area have historically brought instances of both deforestation and reforestation [34]. As in many other mountainous regions inhabited primarily by subsistence farmers, the area has been involved in China's Grain to Green (GTGP) reforestation program since 2000, which aims to provide participants with monetary and/or in-kind payments in return for replacing croplands on steep hillslopes with forestland or grassland [35]. A total of 774 of the area's approximately 3256 households are enrolled in the program, receiving an average of 230 yuan (currently, 1 yuan = 0.14 USD) per mu (15 mu = 1 ha). Through GTGP, the local government has provided farmers with seedlings of pine, Chinese fir, and other tree species, along with bamboo seedlings. Figure 1 contains a true-color aerial image of the reserve in August 2016.

Empirical Data
A time series of Landsat imagery taken on near-anniversary dates covering the study area was selected (earthexplorer.usgs.gov; Landsat level-2 product): 22 August 1990, 22 August 1996 August 2002, 16 August 2011 (Landsat 4/5 TM), and 29 August 2016 (Landsat 8 OLI). These dates were chosen because they provided mostly cloud-free near-anniversary dates during summer [36] at intervals of approximately five years [37]. Clouds were identified and masked out across the time series using the cloud/quality layer provided within the Landsat Level-2 product. A relative atmospheric correction was applied to the time series imagery using the empirical line correction [38], pseudo-invariant targets (i.e., cells whose land cover remained invariant across the time series), and

Empirical Data
A time series of Landsat imagery taken on near-anniversary dates covering the study area was selected (earthexplorer.usgs.gov; Landsat level-2 product): 22 August 1990, 22 August 1996, 31 August 2002 August 2011 (Landsat 4/5 TM), and 29 August 2016 (Landsat 8 OLI). These dates were chosen because they provided mostly cloud-free near-anniversary dates during summer [36] at intervals of approximately five years [37]. Clouds were identified and masked out across the time series using the cloud/quality layer provided within the Landsat Level-2 product. A relative atmospheric correction was applied to the time series imagery using the empirical line correction [38], pseudo-invariant targets (i.e., cells whose land cover remained invariant across the time series), and the 1990 Landsat TM image Remote Sens. 2020, 12, 618 4 of 18 as the base image. To ensure the integrity of the empirical atmospheric correction, pseudo-invariant targets that significantly diverged from the empirical functions were removed until each function had R 2 values greater than 0.85. Model coefficients obtained for these empirical atmospheric corrections are available in Appendix A. The resulting equations were then applied to the spectral bands of the other images in the time series. Following this radiometric correction, the wide dynamic range vegetation index (WDRVI) was calculated across each image in the time series WDRVI = 0.2×ρ NIR −ρ red 0.2×ρ NIR +ρ red . The WDRVI was chosen because it is sensitive to variation in greenness values when vegetation density is high, whereas other indices saturate [39]. In addition to the Landsat-derived WDRVI, data on elevation were also utilized. These data were derived from a digital elevation model obtained by the Shuttle Radar Topography Mission (earthexplorer.usgs.gov), from which slope and slope aspect were calculated. Elevation, slope, and slope aspect were assumed constant throughout the time series.
Data on wild boar populations were obtained using 69 camera traps (Bushnell Trophy Cam) deployed between April 2015 and November 2016 throughout the study area [40]. Cameras were placed in each of 71 sampling plots (20 m × 20 m), which had been selected to spread across most of the reserve area. In total, 69 of these camera trap placements produced usable observations. While accessibility issues inhibited true random placement, sampling plots were selected based on distance to other plots, elevation, and the advice of reserve staff and local guides. To fit with the anniversary dates of the Landsat time series, which are ideally collected in summer due to the season's phenological stability [41], only camera trap observations from summer were used (1 June-31 August). Boar density at each camera trap location was estimated by dividing the total number of boars sighted by the number of summer days that camera was in operation. To prevent the same boar from being counted more than once in rapid succession, it was assumed that a spurt of snapshots captured within 10 min were of the same boar unless more than one appeared in a single image. Where there was more than one per image, the number of boars was assumed to be the greatest number captured in one image during that cluster. A continuous estimate of this proxy for population density was then derived across the entire study site using the inverse distance weighting interpolation (IDW) procedure, using three different k-values: 1, 2, and 3. This set of k-values is typical in environmental research [42] and allowed us to optimize model calibration by selecting the power that produced the most transferrable results. This produced three separate surfaces, each estimating boar population density across the landscape; k-values refer to the power assigned to each point's distance from the cell in question; when k = 1, the weight would be proportional to the inverse distance; when k = 2, the weight would be proportional to the inverse distance-squared, and so on. Thus, the surface derived with k = 1 gave the most weight to camera traps farther away and those with k = 3 gave the least. Figure 2 illustrates the estimated boar density distribution produced by each k-value. ). The WDRVI was chosen because it is sensitive to variation in greenness values when vegetation density is high, whereas other indices saturate [39]. In addition to the Landsat-derived WDRVI, data on elevation were also utilized. These data were derived from a digital elevation model obtained by the Shuttle Radar Topography Mission (earthexplorer.usgs.gov), from which slope and slope aspect were calculated. Elevation, slope, and slope aspect were assumed constant throughout the time series. Data on wild boar populations were obtained using 69 camera traps (Bushnell Trophy Cam) deployed between April 2015 and November 2016 throughout the study area [40]. Cameras were placed in each of 71 sampling plots (20 m × 20 m), which had been selected to spread across most of the reserve area. In total, 69 of these camera trap placements produced usable observations. While accessibility issues inhibited true random placement, sampling plots were selected based on distance to other plots, elevation, and the advice of reserve staff and local guides. To fit with the anniversary dates of the Landsat time series, which are ideally collected in summer due to the season's phenological stability [41], only camera trap observations from summer were used (1 June-31 August). Boar density at each camera trap location was estimated by dividing the total number of boars sighted by the number of summer days that camera was in operation. To prevent the same boar from being counted more than once in rapid succession, it was assumed that a spurt of snapshots captured within 10 min were of the same boar unless more than one appeared in a single image. Where there was more than one per image, the number of boars was assumed to be the greatest number captured in one image during that cluster. A continuous estimate of this proxy for population density was then derived across the entire study site using the inverse distance weighting interpolation (IDW) procedure, using three different k-values: 1, 2, and 3. This set of k-values is typical in environmental research [42] and allowed us to optimize model calibration by selecting the power that produced the most transferrable results. This produced three separate surfaces, each estimating boar population density across the landscape; k-values refer to the power assigned to each point's distance from the cell in question; when k = 1, the weight would be proportional to the inverse distance; when k = 2, the weight would be proportional to the inverse distance-squared, and so on. Thus, the surface derived with k = 1 gave the most weight to camera traps farther away and those with k = 3 gave the least. Figure 2 illustrates the estimated boar density distribution produced by each k-value. Each camera trap location was also assigned to one of five vegetation classes recorded based on in-situ observation: evergreen broadleaf, bamboo, coniferous, mixed broadleaf, and deciduous. The camera trap locations were converted to Thiessen polygons to estimate vegetation cover across the study area. Camera trap coordinates, dates of operation, and vegetation classes are listed in Appendix B. From here, the estimated boar density and vegetation type were assigned to each of 1527 agricultural plots represented as points-73 of which were removed from analysis because they fell outside the range of the Thiessen polygons and could not be reliably assigned to a forest type. In addition, 156 parcels were removed due to missing household allocations. This left 1298 points (i.e., parcels) in the model. Each point was also assigned the zonal mean WDRVI, elevation, slope, and slope aspect within a 500 m radius; 500 m was chosen to reflect the typical 500-600 m home range of a wild boar [43]. Each camera trap location was also assigned to one of five vegetation classes recorded based on in-situ observation: evergreen broadleaf, bamboo, coniferous, mixed broadleaf, and deciduous. The camera trap locations were converted to Thiessen polygons to estimate vegetation cover across the study area. Camera trap coordinates, dates of operation, and vegetation classes are listed in Appendix B. From here, the estimated boar density and vegetation type were assigned to each of 1527 agricultural plots represented as points-73 of which were removed from analysis because they fell outside the range of the Thiessen polygons and could not be reliably assigned to a forest type. In addition, 156 parcels were removed due to missing household allocations. This left 1298 points (i.e., parcels) in the model. Each point was also assigned the zonal mean WDRVI, elevation, slope, and slope aspect within a 500 m radius; 500 m was chosen to reflect the typical 500-600 m home range of a wild boar [43]. Figure 3 illustrates vegetation distribution and camera trap and parcel locations, with the reserve boundary in white. In-person interviews were conducted with households in a stratified random sample. Based on a 2013 census that had identified 3256 households, 1160 households were selected in hopes of obtaining 650 usable interviews after eliminating households with no knowledgeable member present due to travel or other circumstances. The target of 650 households was selected somewhat arbitrarily to ensure satisfactory statistical power. The 3256 households were divided into 123 sampling units; 58 of these units were selected and assigned to 20 administrative villages in proportion to each village's population size. This produced a slight overrepresentation of smaller villages. In total, 20 households were selected from each administrative village per sampling unit; 1160 households were ultimately selected, and full surveys were completed for 605 households in 2014. In 2015, these 605 households were revisited, and 494 full surveys were completed [44]. Details are available at http://complexities.org/pes/research/recent-updates. The household head was interviewed when possible; otherwise any other available, knowledgeable adult was interviewed. Each interviewee indicated the level of crop raiding over the past twelve months on a five-point Likert scale where 1 = "no crop raiding," 2 = "non-serious crop raiding," 3 = "somewhat serious crop raiding," 4 = "serious crop raiding," and 5 = "very serious crop raiding." The Likert scale is useful for analyzing ordinal data like this when the distances between values, i.e., perceived seriousness, cannot be practically measured [45]. Respondents also answered whether the severity had increased, In-person interviews were conducted with households in a stratified random sample. Based on a 2013 census that had identified 3256 households, 1160 households were selected in hopes of obtaining 650 usable interviews after eliminating households with no knowledgeable member present due to travel or other circumstances. The target of 650 households was selected somewhat arbitrarily to ensure satisfactory statistical power. The 3256 households were divided into 123 sampling units; 58 of these units were selected and assigned to 20 administrative villages in proportion to each village's population size. This produced a slight overrepresentation of smaller villages. In total, 20 households were selected from each administrative village per sampling unit; 1160 households were ultimately selected, and full surveys were completed for 605 households in 2014. In 2015, these 605 households were revisited, and 494 full surveys were completed [44]. Details are available at http://complexities.org/pes/research/recent-updates. The household head was interviewed when possible; otherwise any other available, knowledgeable adult was interviewed. Each interviewee indicated the level of crop raiding over the past twelve months on a five-point Likert scale where 1 = "no crop raiding," 2 = "non-serious crop raiding," 3 = "somewhat serious crop raiding," 4 = "serious crop raiding," and 5 = "very serious crop raiding." The Likert scale is useful for analyzing ordinal data like this when the distances between values, i.e., perceived seriousness, cannot be practically measured [45]. Respondents also answered whether the severity had increased, decreased, or remained constant over the past 10 years. These interviews also covered each household's demographics, work and migration history, participation in GTGP, experience with crop raiding, agricultural holdings, and other topics related to livelihood and lifestyle. This produced a total of 605 household interviews; in 2015, all 605 households were revisited to provide more information on land use and assist in participatory mapping of their agricultural holdings. A total of 494 household interviews were completed in 2015. Each household identified the location of its agricultural parcels on a local map, which were later converted to point coordinates. This provided a dataset with 1298 agricultural fields, each assigned to a household with the household's respective survey responses.

Data Analyses
To evaluate the ability to predict boar density from WDRVI, slope, elevation, and aspect, separate Ordinary Least-Square (OLS) regression analyses were run for each vegetation type and for all five vegetation types combined, using the three interpolated boar densities (obtained using three different k-values in the IDW) as dependent variables. Prior studies have shown OLS regressions are useful for habitat suitability models [46] despite imperfectly linear relationships and may even fit the data better than alternatives like nonlinear quantile regression [47]. While a nonlinear function like a general additive model may have provided slightly more accurate predictions of boar density at the validation parcels based on the calibration parcels, because the OLS stage of this analysis was simply meant to test the usefulness of these variables at each vegetation type before geographically weighted analysis later on, we opted for the simplicity of OLS. Separate regressions for each vegetation type were run to account for differences in how attracted boar are to different types [24]. These OLS models were calibrated using clustered random samples of agricultural fields for each vegetation class; parcels were chosen via a random number generator where the parcel with the unique ID matching that number, and those with the next four unique IDs, were selected. This clustering technique was selected to save time over simple random sampling of one parcel at a time, while the clusters were small enough (5 out of 1298 parcels at a time) to help ensure wide distribution and relative independence among selected parcels. This selection process was repeated until 25% of the parcels had been chosen; these would become the "validation" group and the rest would form the "calibration" group. The linear models derived from each calibration group were then applied to the validation parcels within the same vegetation class, and a new OLS regression was run to relate predicted to observed (i.e., interpolated) boar density. A positive, statistically significant relationship between predicted and observed values indicated that the model was useful in predicting boar density. Closer slope estimates to 1.0, together with intercepts closer to zero, showed that the model exhibited higher accuracy and transferability. After calibrating and validating models for each k-value and each vegetation type, we needed to select the k-value that best reflected the dispersal of wild boar around the camera traps. The best k-value is that for which predicted and measured values are the closest, which can be answered using R 2 [48]. Regressions based on the k = 1 interpolated boar density produced the highest R 2 values in calibration for three of the five forest types and the combined model. It also provided the highest R 2 values when validating the regressions for four of the five forest types and the combined model. Thus, the IDW interpolation using k = 1 was used for the remainder of analysis. Models obtained using the calibration and validation datasets are listed in Table 1.  After investigating whether remotely sensed variables could be used to predict wild boar density, a new series of regression analyses was run to predict boar population density over time. First, the same variables from the OLS models were put into a geographically weighted regression (GWR) model to estimate a spatially variable model of boar density based on vegetation type, WDRVI, elevation, slope, and aspect-derived moisture index in the 500 m surrounding each agricultural field. GWR was selected as it is likely to provide more spatially accurate estimations of the effects of environmental variables than global OLS models could when spatial variability in the research question is of importance [49].

Evergreen broad
(Slope was not included as a coefficient for deciduous forests because its values were too clustered for the software to process without error; we deemed this an acceptable omission given slope's lack of significance in the OLS model for deciduous forests.) These models were calibrated using the WDRVI values from the 2016 Landsat image and the slope, aspect, and elevation values from the DEM, then the point-specific coefficients were applied to the 1990, 1996, 2002, and 2011 Landsat images to estimate boar exposure on each field during those years. The model was applied across multiple years to estimate whether boar populations had increased or decreased overall over time and whether our modeled population densities correlated with crop raiding reported by farmers. These years were selected because they provided mostly cloud-free near-anniversary dates and covered periods before and after the implementation of GTGP in 2000. The boar density value at each parcel for each year was then subtracted from the GWR-predicted value based on the 2016 image, producing four new variables estimating change in boar density between the given year and 2016. This was an appropriate way to estimate historical boar distributions because habitat suitability, which can be modeled through indicators like vegetation index and topography, is a widely used indicator of population distribution [12,16]. Figure 4 illustrates the methodology for estimating boar density in current and past years. This dataset was then joined by household ID with the 2015 household survey dataset (n = 494). These records went into a random effects logistic regression grouped by household, where the dependent variable was whether or not the household had experienced "serious" or "very serious" crop raiding over the past 12 months (scores 4-5 on the Likert scale), and the independent variable was the GWR-predicted boar density from the 2016 image.   Table 2 shows the mean estimated boar density for each vegetation type. Figure 5 shows the percentage of 446 households that rated crop raiding at each level of seriousness based on the Likert scale from "none" to "very serious," as well as the percentages who believed crop raiding had lessened, remained steady, or increased over the past 10 years. 40.2% believed crop raiding had worsened in the past 10 years, 21.4% believed it had remained steady, and 28.4% believed it had lessened. Figure 6 illustrates how these responses were distributed across the reserve. We used a survey conducted on the same sample the previous year (2014) to confirm that wild boars were the most-blamed species for crop raiding, with 93.7% of respondents who experienced crop raiding stating that wild boar was the species that caused the most damage. This confirmed boar population was an appropriate proxy for crop damage risk. "Boar per day" estimates for each camera trap ranged from 0 to 0.18; 48 of the 71 camera traps did not record any boar during summer, and the highest number of boars observed in one day was four.  Table 2 shows the mean estimated boar density for each vegetation type. Figure 5 shows the percentage of 446 households that rated crop raiding at each level of seriousness based on the Likert scale from "none" to "very serious," as well as the percentages who believed crop raiding had lessened, remained steady, or increased over the past 10 years. 40.2% believed crop raiding had worsened in the past 10 years, 21.4% believed it had remained steady, and 28.4% believed it had lessened. Figure 6 illustrates how these responses were distributed across the reserve. We used a survey conducted on the same sample the previous year (2014) to confirm that wild boars were the most-blamed species for crop raiding, with 93.7% of respondents who experienced crop raiding stating that wild boar was the species that caused the most damage. This confirmed boar population was an appropriate proxy for crop damage risk. "Boar per day" estimates for each camera trap ranged from 0 to 0.18; 48 of the 71 camera traps did not record any boar during summer, and the highest number of boars observed in one day was four.  On average, estimated boar density decreased slightly within each 500 m radius around agricultural fields between 1990 and 2016. Across all vegetation types, predicted boar density decreased by a mean of 0.00005 BPD (SD 0.0002). However, trends varied widely by vegetation type. Mean change in boar density by agricultural field since 1990, 1996, 2002, and 2011, along with the maximum increase and decrease, are listed by vegetation type in Table 3. Predictive model accuracies also varied by vegetation type. When a single model was applied to all forest types at once, it performed respectably (R 2 = 0.45) and generally followed the correlations found in separate models. Across the entire study area, boar density was positively explained by vegetation index (+0.048 BPD per unit WDRVI; p < 0.01), positively explained by slope (+0.0005 BPD per degree; p < 0.01), negatively explained by elevation (−0.004 BPD per 100 m; p < 0.01), and negatively explained by aspect-derived moisture index (−0.0010 BPD per point; p < 0.01). This suggests boars overall prefer lower, dryer, steeper slopes with denser vegetation. The model was strongest in conifer forests with high internal consistency (i.e., accuracy of the model within the spatial range of calibration sites) at (R 2 = 0.68) and was well validated (β1 = 0.974; β0 = 0.0007; p < 0.01; R 2 = 0.55). In conifer forests, boar density was positively correlated with vegetation index (+0.014 BPD per one-unit increase in WDRVI; p < 0.01), positively correlated with slope (+0.0001 BPD per degree; p < 0.05), negatively correlated with elevation (−0.003 BPD per 100 m; p < 0.01); and negatively correlated with aspect-derived moisture index (−0.0002 BPD per point; p < 0.01). The bamboo model was moderately strong; there was moderate internal consistency within the model (R 2 = 0.49) and it fit the validation parcels well (β1 = 1.460; β0 = -0.009; p<0.01; R 2 = 0.51). In bamboo forests, boar density was positively correlated with  On average, estimated boar density decreased slightly within each 500 m radius around agricultural fields between 1990 and 2016. Across all vegetation types, predicted boar density decreased by a mean of 0.00005 BPD (SD 0.0002). However, trends varied widely by vegetation type. Mean change in boar density by agricultural field since 1990, 1996, 2002, and 2011, along with the maximum increase and decrease, are listed by vegetation type in Table 3. Predictive model accuracies also varied by vegetation type. When a single model was applied to all forest types at once, it performed respectably (R 2 = 0.45) and generally followed the correlations found in separate models. Across the entire study area, boar density was positively explained by vegetation index (+0.048 BPD per unit WDRVI; p < 0.01), positively explained by slope (+0.0005 BPD per degree; p < 0.01), negatively explained by elevation (−0.004 BPD per 100 m; p < 0.01), and negatively explained by aspect-derived moisture index (−0.0010 BPD per point; p < 0.01). This suggests boars overall prefer lower, dryer, steeper slopes with denser vegetation. The model was strongest in conifer forests with high internal consistency (i.e., accuracy of the model within the spatial range of calibration sites) at (R 2 = 0.68) and was well validated (β1 = 0.974; β0 = 0.0007; p < 0.01; R 2 = 0.55). In conifer forests, boar density was positively correlated with vegetation index (+0.014 BPD per one-unit increase in WDRVI; p < 0.01), positively correlated with slope (+0.0001 BPD per degree; p < 0.05), negatively correlated with elevation (−0.003 BPD per 100 m; p < 0.01); and negatively correlated with aspect-derived moisture index (−0.0002 BPD per point; p < 0.01). The bamboo model was moderately strong; there was moderate internal consistency within the model (R 2 = 0.49) and it fit the validation parcels well (β1 = 1.460; β0 = -0.009; p<0.01; R 2 = 0.51). In bamboo forests, boar density was positively correlated with  Table 3. Predictive model accuracies also varied by vegetation type. When a single model was applied to all forest types at once, it performed respectably (R 2 = 0.45) and generally followed the correlations found in separate models. Across the entire study area, boar density was positively explained by vegetation index (+0.048 BPD per unit WDRVI; p < 0.01), positively explained by slope (+0.0005 BPD per degree; p < 0.01), negatively explained by elevation (−0.004 BPD per 100 m; p < 0.01), and negatively explained by aspect-derived moisture index (−0.0010 BPD per point; p < 0.01). This suggests boars overall prefer lower, dryer, steeper slopes with denser vegetation. The model was strongest in conifer forests with high internal consistency (i.e., accuracy of the model within the spatial range of calibration sites) at (R 2 = 0.68) and was well validated (β 1 = 0.974; β 0 = 0.0007; p < 0.01; R 2 = 0.55). In conifer forests, boar density was positively correlated with vegetation index (+0.014 BPD per one-unit increase in WDRVI; p < 0.01), positively correlated with slope (+0.0001 BPD per degree; p < 0.05), negatively correlated with elevation (−0.003 BPD per 100 m; p < 0.01); and negatively correlated with aspect-derived moisture index (−0.0002 BPD per point; p < 0.01). The bamboo model was moderately strong; there was moderate internal consistency within the model (R 2 = 0.49) and it fit the validation parcels well (β 1 = 1.460; β 0 = −0.009; p < 0.01; R 2 = 0.51). In bamboo forests, boar density was positively correlated with slope (+0.003 BPD per degree; p < 0.01), negatively correlated with elevation (−0.008 BPD per 100 m; p < 0.01), and uncorrelated with vegetation index or aspect. In evergreen broadleaf domains, where the model had a somewhat low internal consistency (R 2 = 0.18), boar density was positively correlated with slope (+0.0004 BPD per degree; p < 0.01), negatively correlated with aspect-derived moisture index (−0.0008 BPD per point; p < 0.01), and uncorrelated with vegetation index and elevation. The validation performed similarly well to the calibration (β 1 = 1.226; β 0 = −0.008; p < 0.01; R 2 = 0.19). Model R 2 for mixed broadleaf forests was also moderately low (R 2 = 0.21), and the model applied less effectively to validation parcels (β 1 = 1.518; β 0 = −0.010; p < 0.01; R 2 = 0.38). Here, boar density was positively correlated with vegetation index (+0.061 BPD per unit WDRVI; p < 0.01), negatively correlated with elevation (−0.005 BPD per 100 m), and uncorrelated with slope and aspect. The deciduous forest model had the highest internal consistency; with an R 2 of 0.98, it showed boar density positively correlated with vegetation index (+0.066 BPD per unit WDRVI; p < 0.01), negatively correlated with elevation (−0.007 BPD per 100 m; p < 0.01), and uncorrelated with slope or aspect. Despite the high internal consistency, model validation failed; the relationship between predicted and interpolated boar density was negative. This was unsurprising given the smaller sample size and dense clustering of most deciduous plots. Overall, estimated boar density decreased slightly across the study site between 1990 and 2016; it decreased in evergreen broadleaf, conifer, and mixed broadleaf forests, but increased in bamboo and deciduous domains. Still, the direction and magnitude of change varied considerably within each vegetation type.   Geographically weighted regressions reflected survey-reported perceptions of crop raiding in some vegetation types. In conifer and deciduous forests, higher estimated boar density correlated with a greater likelihood that a householder had experienced "serious" or "very serious" crop raiding over the previous twelve months. Relationships between predicted boar density and reported crop raiding severity were not statistically significant for evergreen broadleaf, bamboo, or deciduous forests. Probability of reporting crop raiding as "serious" or "very serious" is illustrated in Figure 7. Changes in estimated boar density also correlated with reported changes in severity for some years in evergreen broadleaf, conifer, and deciduous forests. In evergreen broadleaf forests, respondents who had experienced greater decreases in boar density since 1990 or 1996 were more likely to report crop raiding having decreased during the past 10 years. This was also the case for conifer and deciduous forests that had experienced greater decreases in estimated boar density since 1996. Estimated changes in boar density did not correlate with the likelihood of reporting lessened crop raiding in bamboo or mixed broadleaf forests for any year. Figure 8 illustrates the probability a respondent reported crop raiding as decreasing if predicted boar density had decreased by 0.001 BPD between each comparison year and 2016. All logistic regression models are listed in Table 4.

Results
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 17 respondent reported crop raiding as decreasing if predicted boar density had decreased by 0.001 BPD between each comparison year and 2016. All logistic regression models are listed in Table 4. Figure 7. Probability of reporting crop raiding as "serious" or "very serious" by vegetation and two plausible boar densities (* p < 0.05; ** p < 0.01).

Discussion
Our results suggest remotely sensed variables are useful for estimating boar population density and crop raiding severity under certain vegetation types. Bamboo and coniferous forests produced models with moderate to high internal consistency and transferability (moderately high R 2 in the OLS calibration and validation models within each forest type). In these vegetation types, remotely sensed vegetation indices and topography seem to be useful for estimating boar density with minimal costs, and it is also feasible to estimate year-to-year variations in wild boar population density. However, this should be further verified with multiyear ground estimates of boar populations. Although the modeling strategy was less internally consistent and transferrable for mixed and evergreen broadleaf forests than it was for bamboo and coniferous forests, it still may provide some useful information, but results from this model need to be interpreted more cautiously. Furthermore, the model's lack of significant relationship with vegetation index in evergreen broadleaf and bamboo domains suggests this technique would not be useful in providing multi-year estimates there because vegetation index is the only predictive variable that can change substantially from year to year. The model for deciduous forests, although exhibiting extremely high internal consistency, exhibited no transferability as the relationship between predicted and interpolated values was negative and small. This contradiction may have occurred due to the smaller sample size and densely clustered pattern of the deciduous units; most of these fields were densely packed at the south end of the reserve with a few fields scattered along the northern end. The randomly selected calibration parcels all happened to occur at the south end, so it is unsurprising that the model would have poor predictive ability in parcels located more than 30 km away and whose vicinity was not represented in the calibration. With better dispersed survey units, it may be feasible to derive a usable model for deciduous forests using this method. It should also be noted that camera traps were generally placed in more geographically accessible locations due to the region's rough terrain [40]; this may have limited our model's usability in some regions.
While this dataset does not allow us to verify crop damage or historical boar populations through ground estimates, household survey responses offer insights that can be used to corroborate the models. In conifer and deciduous forests, households whose fields showed higher risk of crop raiding by wild boars were more likely to report having experienced "serious" or "very serious" crop raiding in the past 12 months. This suggests that the model constitutes a useful predictor of boar density in these forest types. In conifer, evergreen broadleaf, and deciduous forests, households whose fields had higher estimated decreases in boar exposure since 1996 were more likely to report that crop raiding had decreased in the past 10 years; this was also true for changes since 1990 in evergreen broadleaf forests. While there was no correlation between the estimated boar densities since 2002 or 2011 and the reported crop raiding trajectory, this does not necessarily undermine the hypothesized relationship during the decade in question. It cannot be assumed that respondents confined their thought processes to the 10 years specified in the survey (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), as it is common for respondents to "telescope" prior years into a given period, especially regarding routine, non-landmark events [50]. Their responses may thus reflect changes prior to the decade defined by the survey if changes before 2005 are "telescoped" in. Further, although our results suggest our modeling framework best captures population conditions in conifer forests, it is important to underline that these forests have the lowest boar densities among the vegetation types evaluated despite conifer forests' wide geographic coverage and high sample size (n = 378). While this may reduce the model's utility at this study site, it shows our methods are more useful in conifer-dominated landscapes that experience substantial crop raiding by wild boars.
The results obtained in this study address a critical issue in human-wildlife conflict: that of gaps between human perceptions of wildlife actions and biophysical realities [51]. Human-wildlife conflict is fraught with misconceptions, and farmers often overestimate the threat that wildlife pose to their crops [52]. This can exacerbate the resistance to conservation efforts and environmental gains [3,5]; it may also lead to misdirected efforts to mitigate farmers' losses and reduce tensions. However, this study's correlations between perceived crop raiding occurrence, and its temporal trajectories, with modeled boar population density provides preliminary evidence that perceptions and reality are often in sync, at least within our study area. Under these circumstances, technical assistance to reduce crop raiding may effectively reduce tensions and preserve cooperation with conservation actions. Of course, given that this dataset does not measure actual crop damage or ask landholders for more detailed information on their experiences, this study does not provide sufficient proof on its own. Further research with objective, physical measurements of crop damage, together with more detailed information on landholders' perceptions, are needed to verify correlations among modeled boar density, actual crop raiding, and perceived crop raiding. Nevertheless, this study does provide nascent support for the validity of crop raiding perceptions as well as a methodological skeleton for further investigation with more detailed data.
The authors emphasize this study is exploratory, intended more to test the feasibility of the methodological framework than to provide concrete answers on wild boar and crop raiding at Fanjingshan. Along these lines, it is important to note that forecasting based on population estimates from such a limited timeframe is dubious; this methodology would be better employed with camera trap data from several years, ideally spread out. Adding temporal depth to the camera trapping estimates would also allow for the inclusion of critical climatic data like year-to-year variations in temperature and moisture, which may vastly improve habitat and population predictions. This study indeed demonstrates the utility of camera traps, remotely sensed images, and household surveys in deriving empirical relationships between boar distributions and present-day crop raiding burdens. Meanwhile the correlations between estimated boar distributions over time and landholders' perceptions of worsening or abating crop raiding suggest the modeling techniques described here may be useful for updating habitat maps and planning interventions based on time-variant, remotely sensed variables like vegetation and weather. Creating reliable forecasts of boar population and crop damage using these methods will, of course, call for richer data from which to calibrate the models. Future studies should refine this methodology with multi-year camera trap measurements (or other reliable population estimation techniques), weather data, and physical measurements of crop damage that do not rely solely on landholders' perceptions. After investing in a robust foundation of data, this modeling technique may provide an effective, affordable means of managing wildlife and their damage to crops as populations fluctuate and shift over time.

Conclusions
This study uncovers the potential to improve crop raiding monitoring and management over multi-decade periods at minimal cost after initial ground population estimates are made. This may help design cost-effective, easily updatable compensation schemes and technical interventions to minimize economic burdens to farmers. It also empirically affirms anecdotal understandings of crop raiding's relationship to regional ecological changes in some vegetation domains. Although the model suggests crop raiding is not increasing for all households, it demonstrates fine-scale heterogeneity in both realities and perceptions. While the boar density models obtained for Fanjingshan cannot (and should not) be applied to other areas directly [13], the methods for deriving geographically weighted coefficients described here may allow for cost-effective, long-term monitoring of wild boar populations and crop raiding risks at other sites around the world. Given the expense of continuous ground measurements and the limited resources in many affected communities, this may improve boar management and compensation arrangements at minimal cost by allowing practitioners to update distribution maps using remotely sensed imagery.