Abstract
The illegal ivory market poses a serious threat of extinction to the African elephant (Loxodonta africana). Conservation efforts to protect elephants are challenging due to their vast habitat range. Conservation efforts focused on high risk areas provide a more efficient management method than randomly patrolling entire habitat areas. The use of unmanned aerial vehicles is having a profound effect on conservationists' abilities to detect and stop poachers. This study provides methods for identifying high risk elephant poaching areas and for modeling drone surveillance capabilities. Point pattern analyses were conducted to identify spatial distribution patterns of elephant poaching incidents in the Tsavo National Parks area in Kenya. We performed geospatial analyses on the physical environment to create a risk map based on how roads, water features, and land cover correlate to poaching incidents. We also modeled drone flight paths for surveillance of high risk areas based on aircraft flight characteristics and the horizontal field of view for a selected infrared camera. We conclude that poaching incidents were geographically clustered and followed a deterministic (predictive) process. Poaching incidents correlated with close proximity to roads and water features, and were predominately upon specific land cover types. Drone flight modeling determined the number of flights and duration necessary to cover a specific high risk area. The procedures discussed can be applied using a combination of both GIS software and freely available statistical analysis software. The analysis and risk identification will enable conservation groups with limited budgets to improve the efficiency of their anti-poaching efforts.
Introduction
World elephant populations are declining rapidly due to illegal poaching and trade in ivory. This illegal ivory trade is fueled mostly by East Asia's demand for ivory trinkets and statues, which are a symbol of wealth and social status. High demand has driven black market values for raw elephant tusks to approximately US$1,700 per pound while finished ivory art pieces can sell for as high as US$830,000 each [1]. Global trade for illegal wildlife generates approximately US$10 billion per year, making it the fifth most profitable elicit trade [2]. In addition to threatening the existence of elephants, profits from the illegal ivory trade have been linked to funding crime syndicates and terrorist groups in Africa [3]. The African elephant (Loxodonta africana) is being illegally killed at unsustainable rates, which resulted in a 64 % reduction of the Central African elephant population from 2002 to 2012 [4]. While other human activities also negatively affect the elephant population, such as retaliatory killings by farmers suffering damaged crops, and habitat loss from development [5], poaching for ivory is currently the greatest threat.
Preventing poaching is difficult due to the vast expanse of habitat elephants traverse. African elephants have been estimated to routinely travel between 30 km and 60 km in a single day [6]. Tracking wildlife and intercepting poachers is more challenging when the majority of the area is covered in dense vegetation and response forces are limited to ground vehicles. The advancement of unmanned aerial vehicles (UAV), better known as drones, is revolutionizing the ability to track and identify poachers. Many conservation planners view drones as the future of conservation, and estimate the value of one drone to that of 50 rangers [7]. Equipped with infrared sensors and high definition cameras, drones can detect poachers and wildlife from greater distances and at night. Data they capture are then relayed to ground personnel to assist in directing their movements [8]. A drone's presence also acts to deter poachers as word spreads that flying machines, which can see at night, are reporting their location to rangers [9].
Conventional aerial surveillance of wildlife can cost up to US$50,000 dollars a week, whereas a drone costs as little as US$3,000 per week [10]. This significantly lower cost allows small conservation groups with limited budgets the opportunity to have an aerial surveillance program. Conservation groups with larger budgets also benefit through being able to expand their remote surveillance programs, and reduce the risks to human pilots when flying over dangerous areas.
Although unmanned aerial surveillance programs drastically increase the ability to monitor and protect wildlife, it is still nearly impossible to effectively monitor large elephant ranges within conservation areas. This highlights the challenge to monitoring and response efforts, which need to be focused on areas that have a higher probability of poaching incidents. A recent study in central Africa revealed that prioritizing ranger patrols to focus on the most threatened conservation areas can reduce ranger patrol costs by 63 % [11]. Plumptre et al. [11] analyzed environmental factors (such as land cover, climate variables, elevation data, soil types and distance to rivers) along with human factors (such as distance to roads, distance to villages and presence of agricultural land) to model wildlife distributions and suitable wildlife habitats. They also analyzed the occurrence of illegal activities and ranger patrol areas to identify specific regions that were more susceptible to these activities. Their analysis revealed that the majority of ranger patrols occurred within 3 km of patrol stations, which left 40 % of the conservation areas unpatrolled. Unpatrolled areas contained ideal wildlife habitats and high occurrences of illegal activity. By reassessing patrol stations and routes to focus on areas of suitable wildlife habitats and high threat values, the study reported that the 13,800 km2 study area could be patrolled effectively for 63 % of the cost of patrolling the entire area [11].
For this study we explored the advantages of statistical analysis and Geographic Information Systems (GIS) to support anti-poaching elephant conservation efforts. Our objectives were to: (1) determine the characteristics and predictability of elephant poaching incidents through point pattern analysis, (2) identify high risk poaching areas based on a geospatial analysis of the physical environment, (3) demonstrate the aerial surveillance that could be achieved through drone aircraft flight path modeling, and (4) locate optimal guard station sites by analyzing vehicle travel time over the conservation landscape.
Methods
Study area
Our analyses focused on the protected area complex in the southeastern region of Kenya (Fig. 1) hereafter referred to as the Tsavo ecosystem. The Tsavo ecosystem is composed of the Tsavo East National Park, Tsavo West National Park, Chyulu Hills National Park, Ngai Ndethya National Reserve, Taita Hills Wildlife Sanctuary, and the Lumo Wildlife Sanctuary. The entire ecosystem contains 21,128 km2 of mostly dry, flat, savannah vegetation [12]. The Tsavo ecosystem is home to a wide variety of wildlife including elephant, rhino, buffalo, lion, giraffe, leopard, hippopotamus, crocodile, waterbucks, and over 500 species of birds [12]. Critical to these wildlife populations are the rivers and streams of the protected area complex, including the perennial Galana River and the seasonal Tiva River.
The Tsavo ecosystem in southeastern Kenya. The analyses performed in this study focused on the area within the ecosystem with the exception of the density analyses that included the surrounding Study Area to eliminate edge effects.
The combined Tsavo East and Tsavo West National Parks exceed more than 4 % of Kenya's total land area, and together comprise one of the largest park complexes in the world. The Tsavo is also home to the largest population of Kenya's elephants and is one of four “Monitoring of Illegal Killing of Elephants” (MIKE) sites in Kenya [13]; MIKE sites were created to gather information about threats to elephants to assist in refining management strategies as part of the Convention of International Trade in Endangered Species (CITES) [14].
Our study area boundary (Fig. 1) encompassed the entire Tsavo ecosystem and nearby surrounding area to mitigate for potential edge effects [15]. The mitigation of edge effects was necessary for our point pattern analyses and kernel density estimation of all poaching incidents, as we expected our findings would be significantly different if we excluded poaching incidents outside yet close to the boundary of the protected areas.
Data analysis
We used ESRI's ArcGIS Geographic Information System (GIS) software and the Spatial Analyst Extension to perform the majority of analyses included herein [16]. GIS tools were used (1) to determine the relationship between poaching incidents and physical environmental characteristics in our study area; (2) to model drone surveillance areas; and (3) to determine the strategic placement of ranger guard posts in response to identified high risk areas. We also used open source statistical software, including R-STAT [17] and GeoDa [18], to perform the point pattern analysis and spatial autocorrelation analysis respectively.
Point pattern analysis on elephant poaching incidents
Spatial statistical analysis offers tools for identifying the patterns and characteristics of geographic point data. This type of analysis can reveal whether point distributions are connected to underlying properties of the physical environment (first order effects) and whether a point has a positive or negative effect on the presence of additional points (second order effects) [15]. Spatial statistical analysis can also determine if the point pattern being analyzed is a realization of a deterministic process or a stochastic (random) process. This is a crucial aspect of point pattern analysis given that a stochastic (random) process would not permit predictive modeling. Each of these analyses was performed on the 156 elephant poaching incident points provided by a Kenyan conservation group, and are further described in the following sections.
First order effects – kernel density analysis
As described by O'Sullivan and Unwin [15], first order effects reveal how properties of the physical environment influence the distribution of points at a regional level, and are directly related to density based measurements. The ArcMap Kernel Density tool [16] was used to create an area-wide incident density map to identify specific areas that have a higher concentration of poaching incidents, as well as to determine which, if any, physical properties in the environment correlate with the measured density. The bandwidth, or search radius, for the kernel density tool is user-specified and determines the area that a single data point can influence when calculating density. Large bandwidths allow data points to have more influence on one another over larger areas and are used to analyze data at regional levels. A bandwidth of 25 km was selected to achieve this regional level density analysis.
Second order effects – G and K function distance based measurements
Second order effects describe how individual points in a pattern influence the placement of additional points. Strong second order effects can be interpreted as a dense clustering of points, or the opposite, as equally spaced points that inhibit points from existing near one another [15]. Second order effects are useful in identifying the presence of hotspots, or of concentrated areas of a given activity or data points, and yield insights into the probability of the placement of future data points.
Second order effects are analyzed through distance based measurements. The G function, also referred to as the refined nearest neighbor measurement, measures the distance from each point to its next closest point and then displays the nearest neighbor distance for each point as a cumulative frequency distribution [15]. This resulting graph of a G function analysis displays the fraction of points (G) that have a nearest neighbor distance less than the given distance (r). We performed a G function analysis using R-Stat [17]. To reduce edge effects resulting from not analyzing nearby points outside the study area, we used three separate edge-corrected outputs, which included the Kaplan-Meier, border, and Hanisch corrected estimates [19]. The G function output also includes a random Poisson process for the same area. Areas where the G function outputs for the poaching incident data differ from the random Poisson output indicate possible clustering [19], while steep rises over short distances are similarly indicative of point clustering [15].
We also calculated the K function of poaching incident data using R-Stat [17]. The K function is an additional distance-based measurement useful in determining the specific properties of a point pattern at different scales. It differs from the G function in that it measures the distances from each point to all other points instead of just to the nearest point. This is accomplished by creating a circle of radial distance (r) around each point and tallying other points within each circle. The average number of these points is then divided by total point density for the entire study area to give the K(r) value [15]. This process is then repeated for progressively larger circles. Similar to the G function, R-Stat created three edge corrected outputs for the K function: the Ripley isotropic, translation, and border corrected estimates [20]. The K function output also included a completely random Poisson process for the same area. Areas where the outputs for poaching incident data differ from the Poisson output indicate possible clustering [20].
Deterministic versus stochastic process – The Monte Carlo Procedure and IRP/CSR
An important goal of spatial statistical analysis is to determine if a point pattern results from a deterministic or stochastic process. A stochastic process, also referred to as an independent random process (IRP), or a process of complete spatial randomness (CSR), has an equal probability of a point occurring in any location, and the location of each point occurs independently to the location of other points [15]. These characteristics do not allow the prediction of additional points, as they arise from a stochastic process.
To eliminate the possibility that our data resulted from a stochastic process, we performed 100 Monte Carlo simulations [15] on both the G and K function outputs for the poaching incident data. This procedure generates random point pattern simulations within the study area boundary using the same number of points as the poaching incident data. The procedure then measures the nearest neighbor point distance values of the randomly generated patterns using a selected function, such as the G and K functions. The total function measurement outputs for all simulations are then combined to create an envelope of function output values that would be considered part of an IRP/CSR process, to which the function output of the actual point pattern being analyzed is compared in order to determine if there are distance ranges with more or less clustering than would be expected [15].
Spatial Auto Correlation – Moran's Index and local indicators of spatial association
We used the Moran's Index (or Moran's I) analysis performed with the GeoDa software [18] to determine the degree of spatial autocorrelation for the poaching incident data. The Moran's Index, a technique used for measuring spatial autocorrelation across an entire study area [15], calculates how different a value is at a specific point compared to the mean attribute value of all points, and the similarity of surrounding point attribute values in order to determine if the data is positively or negatively correlated.
To complete this analysis, poaching incident counts needed to be contained within individual sample units to facilitate comparisons. Sampling units were created using the ArcGIS Create Fishnet tool [16] to generate a lattice of grid cells (10 km × 10 km) that served as the sampling units covering the entire Tsavo ecosystem area. Totaling the number of poaching incidents per cell enabled the completion of the Moran's I analysis.
The resulting Moran's I scatter plot points, with higher positive values on the x-axis, are associated with grid cells containing poaching incident counts that were higher than average. Higher positive values on the y-axis indicated that surrounding grid cell poaching incident counts were also above the average, whereas negative axis values indicated cells with below average sums and below average adjacent cells. A scatter plot containing a majority of points in the first quadrant, with positive x and y-axis values, indicated positive spatial autocorrelation. Moran's Index analysis also yields a Moran's I correlation value; values above the 0.3 threshold indicate positive spatial autocorrelation [15].
To assess spatial autocorrelation for specific locations within the study area, we performed a Local Indicators of Spatial Association (LISA) analysis with GeoDa [18], which analyzes individual grid cells and their immediate neighbors. LISA yields a cluster map depicting each of the lattice grid cells based on comparisons of their own incident counts to their neighboring cells' incident counts. Areas that have a higher than average incident rate and have neighboring units that also have higher than average incident rates (“High-High”) contribute to positive spatial autocorrelation.
The LISA procedure also creates a significance map to further expand the cluster map and determine which areas of clustering are statistically significant. Monte Carlo simulations were used to assess the statistical significance of the relationship between a single lattice grid cell and its surrounding neighbors. The Monte Carlo simulation preserves the attribute value of the lattice grid cell being analyzed while running random simulations of the surrounding lattice cells. If the relationship to neighboring values is outside the normal distribution of random values, this indicates a higher significance value for the analyzed lattice grid cell.
Identifying poaching high risk areas based on physical environment and incident intensity
We looked for a positive correlation between physical characteristics of the environment in the Tsavo ecosystem and the occurrence of poaching incidents. In particular, we focused on land cover types, proximity to water sources, and proximity to roads. Poaching risk values were assigned to each physical environment feature based on the total number of poaching incidents that occurred within that feature. These features were selected for analysis based on observations regarding preferred elephant habitat and general poaching behavior.
Land cover type
We analyzed poaching incidents to determine if poachers favored areas of a specific land type according to observed elephant behavior. Evans and Harris [21] found that male African elephants have seasonal preferences for specific land cover types, while female elephant habitat use was observed to be proportional to the total area of each land cover type. As the poaching incident data used in our study does not contain temporal attributes that allow for the analysis of seasonal behavior, land cover preference assumptions were based on female elephant observations.
We identified fifteen land cover types within the Tsavo ecosystem boundary. We calculated the total number of poaching incidents for each land cover type and divided it by the total number of incidents for the entire study area. The percentage of incidents for each land cover type was then converted into a risk number, and assigned to the land cover type for use in the total poaching risk calculation. We assumed that poaching incidents occurring in proportion to the total area of each land cover type would indicate that poachers were specifically targeting areas that have higher elephant counts. Identifying land cover types with low total area but high numbers of poaching incidents would indicate areas preferred by poachers.
Proximity to roads and water features
We analyzed the number of poaching incidents occurring at varying distances to roads to test if poaching behavior is effected by their proximity in the Tsavo ecosystem area. Haines et al. [22] found that poaching incidents are more likely to occur within close proximity to roads due to the easy access to wildlife they provide, and to facilitate escape after poaching has occurred.
In addition, Cushman et al. [23] revealed that the probability of an elephant being in a specific area was found to decrease as a function of the square root of the distance of the elephant to a water source. We therefore also assessed poaching incident data to determine if poachers were exploiting this observed elephant behavior by analyzing the number of poaching incidents that occurred within varying distances from water features.
We used the ArcGIS Euclidian Distance tool [16] to create multiple buffer zones around water features and roads. Buffer distance intervals were set to a value of 2,000 m and covered the entire ecosystem area. The ArcGIS Spatial Join tool [16] was used to determine the number of incidents within each buffer zone. Similar to the land cover type risk values, we calculated the percentage of poaching incidents within each buffer zone and used these to assign each a poaching risk value.
Small bandwidth density analysis
The density of poaching incidents also contributed to determining the overall risk values. Kernel density estimates (KDE) are often used in crime modeling to determine specific areas with high crime intensities [24]. This analysis relied on the ArcGIS Kernel Density Estimate tool [16] using a small bandwidth (search radius) of 7 km which emphasized specific areas of high poaching intensity. Larger risk values were assigned to high poaching density areas, and were used to calculate the total poaching risk.
Calculating total risk
A map of total risk was created to identify varying levels of poaching risk throughout the study area based on the sum of risk values associated with land cover type, proximity to water features, proximity to roads, and incident density using the ArcGIS Raster Calculator tool [16].
Surveillance drone flight paths
Low cost drones are being used successfully in conservation efforts for a wide variety of tasks including wildlife population monitoring and land use change [25]. The increasing popularity of drone use has led to the creation of the free open source flight planning software, ArduPilot Mega Planner, which assists users in determining flight mission capabilities [25]. We used standard factors to model the drone flight paths required to cover a single poaching area that had been determined high risk. The flight path model takes into consideration several key factors, including cruising altitude, cruising speed, and flight time at cruising speed, to determine the total distance a drone could fly on one fueling (battery or gas tank). The surveillance area covered by the drone was determined by the field of view width from the attached camera used for spotting poachers and the line distance the drone is able to fly.
For flight modeling, we used the Hawkeye drone model RQ-84Z AeroHawk [26]. The AeroHawk's cruising altitude can reach a maximum of 91 m, a cruising speed of 59.5 km/h, and a total flight time of 90 minutes at cruising speed, permitting it to cover a total linear distance of 89.3 km. The field of view used in this analysis was based on the FLIR Tau 2 640 thermal imaging camera with a wide field of view 7.5 mm lens. Based on the camera and lens specifications, the horizontal field of view at an altitude of 91 m would be approximately 187 m.
We modeled drone flight paths for a 23.2 km2 high risk region near the center of Tsavo East National Park. To assist in the placement of the flight lines and to ensure total coverage without overlap, we used the ArcGIS Create Fishnet tool [16] to create a sampling grid that covered the selected high risk area. The width and height of each grid cell was set to the total horizontal field of view for the AeroHawk (187 m). In addition to the grid, the Create Fishnet tool also creates a point feature at the center of each grid cell. Each point feature was used as a guide for the point placement of a new flight path. The flight path followed the grid format, and the drone returned to the launch site for refueling before the maximum distance of 89.3 km was reached.
Ground crew station locations
We analyzed the placement of ranger guard stations to determine their response to poaching threats in high risk areas. We based this on analysis of the time it took to travel over the conservation area terrain given the roads, land cover type, and the presence of water features. Roads were assigned the least amount of associated travel cost, while densely forested land covers and river water features were assigned the highest.
Finally, we created a single travel time cost surface layer by combining the individual travel time costs using the ArcGIS Reclassify and Raster Calculator tools [16]. This travel time cost layer was used with the ArcGIS Cost Distance tool [16] to test various ground crew station locations, and to analyze response times to different high risk areas.
Results
Point pattern analysis on elephant poaching incidents
First order effects – kernel density analysis
Density analysis using the large 25 km bandwidth (Fig. 2) reveal that there is an increased density of poaching incidents in the northwestern and southwestern areas of the Tsavo East National Park. In addition, the analysis also helped us identify an area of high density in the eastern portion of the Ngai Ndethya National Park. A comparison of high density poaching regions with data describing the physical environment showed that approximately 60 % of all poaching incidents occurred in areas where the vegetation is largely open low shrubs (65 % - 40 % crown cover). This comparison also showed that higher density areas occur within close proximity to both water and road features. For example, over 62 % of poaching incidents occurred within 4,000 m of a water source, while approximately 89 % of poaching incidents occurred within 4,000 m of a road. Land cover type, proximity to roads, and proximity to water, appear to be strong factors influencing the spatial distribution of poaching incidents and indicate the presence of first order effects.
Large bandwidth (large area), set to 25 km, Kernel Density Estimation results that display poaching incident density at a regional scale. Density results are shown as the number of poaching incidents per km2. Each of the protected areas that make up the Tsavo ecosystem is outlined in red.
Second order effects – G and K Function distance based measurements
The G function analysis indicates that approximately 50 % of points have a nearest neighbor distance less than 3,500 m (Appendix 1-1). A steep rise, indicative of clustering, was seen in the plot of the G function for poaching incidents as the G value increased from 0 to 0.5 between the short distance of 500 m and 3,000 m. The results are nearly identical for the three separate edge-corrected outputs including the Kaplan-Meier, border, and Hanisch corrected estimates which are displayed as black, red, and green lines respectively. Clustering is also indicated by the great difference between the random Poisson process, represented as a blue line, and the three edge-corrected poaching incident outputs. The K function analysis also reveals clustering in the point pattern as indicated by the difference between the random Poisson process (blue line) and the edge-corrected poaching incident pattern output at nearly all radial distances (Appendix 1-2). The K function edge corrected outputs include the Ripley isotropic (black line), translation (red line), and border corrected (green line) estimates.
Deterministic versus stochastic process – The Monte Carlo Procedure and IRP/CSR
The Monte Carlo procedure revealed that point clustering observed in the G function analysis (Gobs) occurred outside the IRP/CSR envelope for most distances (Fig. 3), which indicates that the elephant poaching incident pattern is a deterministic process that should be conducive to predictive modeling. For the K function analysis, the Monte Carlo procedure also shows that observed point clustering (Kobs) occurred outside the IRP/CSR envelope for all radial distances, corroborating the finding that the elephant poaching point pattern is part of a deterministic process. Appendix 1-3 displays the results of the K function Monte Carlo procedure.
Results of Monte Carlo procedure for the G function analysis. The x-axis represents the nearest neighbor measurement distance (r) in meters while the y-axis represents the fraction of points that have a nearest neighbor distance less than the given (r) distance. The black line represents the observed G function results. The combined IRP/CSR simulation values envelope is represented as the grey shaded area, bounded by the high and low simulation values. The red line represents the theoretical value of a random Poisson process for the study area. The observed G function results are outside the random IRP/CSR value envelope for all distances. This is a strong indication that the elephant poaching incident pattern is derived from a deterministic process.
Spatial auto correlation – Moran's Index and local indicators of spatial association
The Moran's Index analysis revealed that the majority of scatter plot points representing grid cell poaching incident counts have positive x and y-axis values, and fall within quadrant one (Appendix 1-4). This indicates spatial autocorrelation. An additional indication of spatial autocorrelation is a Moran's I value of 0.351, which is greater than the 0.3 autocorrelation threshold [15].
The LISA procedure used in this study produced a LISA cluster map and a LISA significance map to identify specific areas of positive spatial autocorrelation (Fig. 4). The cluster map's red “High-High” lattice grid cells, indicating areas contributing to positive spatial autocorrelation, were found primarily in the high density poaching incident areas in the northwest, west central, and southeastern areas of Tsavo East, the eastern region of Ngai Ndethya, and southern region of Lumo.
Local Indicators of Spatial Association (LISA) cluster map and significance map. The cluster map on the left reveals individual sample cells that contribute to positive spatial autocorrelation (clustering). The significance map on the right shows individual cells that are statistically significant after performing the Monte Carlo procedure.
The significance map shows grid cells and their statistical significance after 100 Monte Carlo simulations (Fig. 4). The most significant grid cells (p = 0.001) that also contributed positively to spatial autocorrelation were mainly found in the dense poaching incident areas of northwest and west central Tsavo East Park.
Identify poaching high risk areas based on physical environment and incident intensity
Land cover type
The number of poaching incidents in each land type were correlated with the total area of each land type. For example, the most common land cover type, open low shrubs (65 % - 40% crown cover), occupies approximately 45 % of the Tsavo protected ecosystem area and also has the highest proportion of poaching incidents (60 %). The number of incidents decreased as total area decreased (Appendix 1-5).
High risk values were assigned to the land cover types with the largest total area, including open low shrubs (65 % - 40% crown cover) and trees and shrubs savannah. These two land cover types comprised 74 % of the protected ecosystem area and contained 81 % of the poaching incidents. A complete list of each land type, their incident count, and associated risk value are displayed in Table 1.
Land cover types and their associated poaching incidents, risk value, and percentage of the total Tsavo ecosystem area. The number of poaching incidents strongly correlates to the total area of each land type. The two most abundant land cover types comprise 74% of the protected ecosystem area and contain 81% of the poaching incidents.
Proximity to roads and water features
There was a strong correlation between poaching incidents and proximity to water and road features. Approximately 69 % of poaching incidents occurred < 2,000 m from a road, and 62 % of incidents occurred < 4,000 m from water features. Tabular results showing the number of poaching incidents and associated risk values for each water and road buffer zone are presented in Appendices 2-1 and 2-2 respectively. Risk analysis mapping for water and road features are displayed in Appendices 1-6 and 1-7 respectively.
Small bandwidth density analysis
The small 7 km bandwidth density analysis (Appendix 1-8) reveals specific high poaching intensity hotspots throughout the study area. The northwestern area of the Tsavo East National Park, and the eastern portion of Ngai Ndethya National Park, contained the highest poaching incident densities and were assigned the largest risk values. An isolated area of moderately high density was located in the southeastern area of Tsavo East (Appendix 1-8).
Calculating total risk
The total poaching risk map (Fig. 5) displays the cumulative poaching risk for all areas within the Tsavo ecosystem. Areas with the highest cumulative risk were located on the “open low shrubs land cover type (65 % - 40 % crown cover)” land cover type < 2,000 m from a water source and road, and in hotspot areas of high density incident occurrence. Some of the areas with high total risk values did not have poaching incidents. This is because a high risk value range can result from physical environmental factors alone, enabling the total risk analysis to predict high risk poaching areas without relying entirely on the presence of existing incidents.
Total poaching risk map based on summed risk values. The results of this map will increase the efficiency of anti-poaching efforts by focusing resources on high risk areas.
Surveillance drone flight paths
Flight path modeling indicated that two flights traveling a total of 152.4 km would be required to cover the selected 23 km2 area. Using a cruising speed of 59.5 km/h, surveillance for the selected area would take a total of 2 hours and 34 minutes. Figure 6A shows the grid cells that simulate the horizontal field of view of the UAV and the center points of each grid as created by the ArcGIS Create Fishnet tool. Figure 6B displays how the grid center points are used as guides for the creation of the drone flight path line features.
Fishnet grid and grid center points covering high risk area.
Drone flight lines following grid path. Two flights traveling a total of 152.4 km would be required to cover the selected 23 km2 area. Using a cruising speed of 59.5 km/h, surveillance for the selected area would take a total of 2 hours and 34 minutes.
Ground crew station locations
The location of seven ranger stations and their associated cost distance response times are displayed in Appendix 1-9. These guard stations were located near high risk areas and have an approximate maximum response time of 30 minutes to the surrounding high risk areas. A table of assumed travel speeds is presented in Appendix 2-3.
Flight Path Modeling
Discussion
Accurate prediction of spatial point phenomenon (poaching incidents) requires that the point pattern being analyzed is driven by an underlying deterministic process. Attempting to predict specific areas of future incidents when the context is an IRP/CSR process is not feasible, as any point has an equal probability of occurring at any location in a study area, uninfluenced by the existence of other points [15].
Our initial analyses identified the characteristics of the analyzed poaching incident point pattern through spatial statistical analysis. First order effects using a large bandwidth kernel density analysis revealed that higher density poaching incident areas were influenced by the underlying environmental factors of water features, roads, and land cover type, indicating a nonrandom distribution of points throughout the study area. The observed strong second order effects for poaching incidents as seen through the G and K function distance-based measurements indicated clustering within short distances, suggesting that multiple poaching incidents occurred in localized areas. This type of pattern recognition facilitates a fuller understanding of poaching behavior, and greater accuracy in predicting where future incidents will occur. The Monte Carlo procedure suggested second order effects occurred outside the IRP/CSR envelope for all of the distance-based measurements, evidence that the process behind the point distribution is deterministic, and that extrapolation of future point distributions is a valid and useful endeavor.
The Moran's Index and LISA analyses for spatial autocorrelation advanced our understanding of poaching incident clustering for the entire study area by highlighting the specific sampling unit grid cells that contribute most to positive spatial autocorrelation. The Moran's Index confirmed that autocorrelation of poaching incidents existed within the study areas (MI = 0.351 > 0.3). In addition, both the LISA cluster and significance maps revealed the specific grid cells exhibiting the most positive spatial autocorrelation and statistically significant clustering when compared with Monte Carlo simulations.
A high correlation was found with poaching incidents and each of the physical environmental features in the study area including land cover type, proximity to water, and proximity to roads. Approximately 60 % (93) of all the poaching incidents occurred on the most abundant land cover type open low shrubs (65 % - 40 % crown cover). These results indicate poachers were focusing their activity on land cover types with the largest total area that are known to be preferred by elephants [21], and which dominate the majority of the study area. Poaching incidents also occurred within close proximity to water, suggesting that poachers took advantage of the known elephant migration behavior of staying close to vital water resources [23]. Thirty-four percent (53) of poaching incidents occurred < 2,000 meters from water sources, and 62 % (96) of incidents occurred < 4,000 meters from them. Roads were the physical feature with the greatest correlation to poaching incidents: approximately 69 % (107) of poaching incidents occurred within 2,000 meters of a road. These results concur with previous findings that poachers rely on roads for access to wildlife and for expedited escape after poaching has occurred [22].
The occurrence of poaching incidents in relation to physical environmental features was used to construct the total risk model. These risks were based on the percentage of total poaching incidents occurring within specific land cover types, and the percentage of incidents occurring within multiple buffer zones radiating from water sources and roads at 2,000 m increments. An additional risk layer was created using a small bandwidth (search radius) kernel density analysis that revealed areas of high poaching intensity. The resulting total risk map allows conservation planners to identify specific areas for increased anti-poaching efforts, an essential need for the efficient use of resources for large regional areas.
Drone flight path modeling offers a novel, GIS-based approach to plan for optimal UAV surveillance. The flight model can determine the total distance that a UAV could fly on one fueling (battery or gas tank) and the surveillance area observed by a selected camera. This type of flight modeling can assist conservation planners at the initial stage of deciding on UAV and camera hardware based on the size of a given protected area. The flight model will also help in the deployment of a UAV by providing specific coordinates to use as waypoints with autopilot software such as ArduPilot [25]. If long range drones are selected for use, this flight path modeling would also assist in determining where additional radio control towers would be needed to maintain flight control over large areas. The advantages of UAV use are clear when surveillance flight times are compared to the time a foot patrol would take to cover the same area.
The response time recorded in the guard post analysis demonstrates how GIS methods help enable effective and complete coverage for protected wildlife areas. This type of analysis could help conservation planners determine ideal ground crew station locations and thus minimize response times to high risk areas and increase the ability of ground personnel to intercept poachers. It can also be used to reveal weaknesses in current protection coverage based on existing guard station location response times.
This study would benefit from additional research to generate more detailed data, particularly regarding poaching incidents with related temporal data, such as the date of the observed poaching incident and the estimated carcass age. The ability to differentiate seasonal poaching habits would increase the accuracy and effectiveness of risk modeling and possibly reveal new trends associated with seasonal weather data. Other beneficial data would include the conditions of roads to provide a more accurate estimate of response times from guard stations. The difference in condition between a paved road and an uneven dirt road has a significant effect on response time estimates.
Implications for Conservation
Our results indicate that poaching incidents follow a deterministic, non-random process and that spatial autocorrelation (clustering) is occurring. There is a strong correlation between certain physical features of the landscape and poaching incidents. The majority of incidents occurred in areas where vegetation cover was of open low shrubs with 65 % - 40 % crown cover, and within 2,000 m of roads and 4,000 m of water features.
These strong correlations could contribute to the prediction of future poaching behavior and the identification of high risk areas. As we have shown, high risk areas can be identified for increased drone surveillance and ranger patrols. Haines et al. [22] noted that identifying these high risk areas will allow for a more efficient expenditure of resources to protect threatened wildlife. This type of analysis and risk identification has the potential to enable small conservation groups with limited budgets to protect wildlife within their borders more efficiently.
The use of GIS was integral to this study as it allowed for the visualization and analysis of spatial data. GIS offers unique tools for analyzing spatial data and discovering trends and relationships among various data sets. Many of the methods we used are based on ArcGIS tools and commands that could be reproduced by wildlife conservation groups.
The drone flight modeling methods we used will be of use for the conservation community as drone popularity continues to increase. Aerial surveillance through the use of drones is significantly less expensive compared to using human pilots and full size aircraft. Modeling drone capabilities and understanding their limitations should assist conservation managers in choosing suitable equipment and in developing realistic aerial surveillance programs.
Wildlife poaching is a problem that will not be solved easily. The illegal ivory trade is too profitable to resist for many poachers living in extreme poverty. A recent study identified poaching as most likely to occur in areas of high poverty, where infant mortality rates are high, and standards of living are low [27]. Recent economic growth notwithstanding, Kenya is still an underdeveloped country with minimal resources for combating poaching. Working towards social equality and towards improving the standard of living for people affected by poverty is an important endeavor, as combating poverty is part of combatting poaching.
Footnotes
Acknowledgements
The authors would like to acknowledge the Penn State University Masters of Geographic Information Systems (MGIS) program and instructors for their excellence in graduate level education. The authors would also like to thank the many conservation groups and rangers throughout Africa that are working tirelessly to protect vital wildlife.