Black-tailed prairie dog (BTPD) (Cynomys ludovicianus) populations in Kansas have declined significantly due to both natural and human-induced threats. To minimize the risk of future population declines, it is necessary to identify existing BTPD habitats in the state as well as areas suitable for BTPD habitat. This paper presents a method for modeling BTPD habitats in Kansas using geographic information systems (GIS), remote sensing, and ecological niche modeling with the Genetic Algorithm for Rule-Set Prediction (GARP). Environmental variables incorporated into the ecological niche modeling process include composite biweekly Normalized Difference Vegetation Index (NDVI) layers derived from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite imagery, slope, soil depth, and soil texture. Species occurrence training and validation data were selected from an aerial survey of BTPD colonies by the Kansas Department of Wildlife and Parks (KDWP). Accuracy assessment methods, including Receiver Operating Characteristic (ROC) analysis, omission calculation, and validation with an independent BTPD colony dataset collected for the Cimarron National Grassland in Morton County, indicate a high degree of accuracy for the GARP models. A map of BTPD habitat suitability produced by the ecological niche modeling has the potential to aid state agencies and organizations in their efforts to prevent further population declines in the species.
Over the past century, black-tailed prairie dog (BTPD) (Cynomys ludovicianus) populations in Kansas have decreased significantly due to both natural threats such as bubonic plague as well as human-induced threats including unregulated hunting, poisoning, and the conversion of native grasslands to other types of land cover. In an effort to curtail further declines in BTPD populations and in response to findings by the U. S. Fish and Wildlife Service that the BTPD may warrant protection by the Endangered Species Act, the Kansas Black-Tailed Prairie Dog Working Group was created and recently formulated the “Kansas Black-Tailed Prairie Dog Working Group was created and recently formulated the “Kansas Black-Tailed Prairie Dog Conservation & Management Plan” (Kansas Black-Tailed Prairie Dog Working Group 2002). The primary goal of the plan is “to maintain biologically viable populations of black-tailed prairie dogs at selected sites across the historical range in Kansas” (Kansas Black-Tailed Prairie Dog Working Group 2002: 8). Key objectives of the plan include to “determine and monitor species distribution and status” and to “identify, maintain, and promote existing and additional suitable prairie dog habitats” (Kansas Black-Tailed Prairie Dog Working Group 2002: 8).
In order to achieve key objectives of the “Kansas Black-Tailed Prairie Dog Conservation & Management Plan,” and to prevent future population declines, it is necessary to identify existing BTPD habitats as well as areas in the state with environmental parameters that correspond to high suitability for the BTPD. The primary objective of this paper is to propose a method for modeling habitat suitability for BTPDs in Kansas using geographic information systems (GIS), remote sensing, and ecological niche modeling with the Genetic Algorithm for Rule-Set Prediction (GARP) (Stockwell and Peters 1999). The modeling technique presented is a cost-effective method that may be used to assess, monitor, and manage BTPD habitats.
Black-tailed prairie dogs in Kansas
The historic extent of the BTPD in Kansas includes the portion of the state west of the Flint Hills and the tallgrass prairie ecosystem. Lantz (1903) estimated that the BTPD once occupied over two million acres in sixty-eight of the 105 Kansas counties. However, BTPD populations have declined dramatically over the past century. An aerial survey by the Kansas Department of Wildlife and Parks (KDWP) in 2001 estimated that the BTPD only occupied approximately 130,521 acres in the state (Pontius 2002). Based on a sampling of BTPD towns from the same aerial survey, Pontius (2002) estimated fewer than 5,000 total BTPD towns in western Kansas.
BTPD population declines have had adverse effects on other species, since the BTPD is a keystone species in the shortgrass prairie ecosystem in Kansas and the broader Great Plains region (Agnew, Uresk and Hansen 1986; Desmond, Savidge and Eskridge 2000; Kotliar et al. 1999; Miller, Ceballos and Reading 1994; Wuerthner 1997). Several endangered, threatened, or species of concern that may be found in parts of the Great Plains depend on the BTPD for survival. For example, BTPDs might provide critical food and habitat resources for the black-footed ferret (Mustela nigripes), ferruginous hawk (Buteo regalis), and swift fox (Vulpes velox). Monitoring and assessing BTPD populations may provide key information regarding the status of these other species.
Genetic Algorithm for Rule-set Prediction
Ecological niche modeling is based on the premise that a given species distribution is limited, in part, by the presence of suitable habitat on the landscape. This suitable habitat may be understood as the joint description of the environmental conditions that allow a species to satisfy its minimum requirements (Chase and Leibold 2003), and is defined by various ecological parameters. A benefit of such modeling is the ability to extend a limited number of known habitat locations for a species (or species occurrence data) to “predict” the geographic distribution of the species for a larger area. The Genetic Algorithm for Rule-Set Prediction (GARP) (Stockwell and Peters 1999) provides one such approach to ecological niche modeling. The GARP algorithm has been implemented in the Desktop GARP software of Scachetti-Pereira (2001). GARP is a genetic algorithm that is intelligent, iterative, and based on machine-learning. Ecological niche modeling with GARP requires the selection of two types of GIS data layers for input: 1) environmental variables that collectively characterize a species niche and 2) species occurrence data in the form of x and y geographic coordinates (latitude/longitude, UTM, etc.). A particular strength of GARP is the ability to generate accurate models with a limited number of species occurrence data points (Stockwell and Peterson 2002).
GARP describes relationships between species occurrence data and environmental variables through a number of rule types (atomic, range, negated range, and logistic regression [logit]), and, following internal testing and refinement, develops a set of rules to create a model (Stockwell and Peters 1999). GARP refines the rule-set until either the maximum number of iterations or model convergence is achieved, and outputs a user-defined number of models for each species. Each model is a raster map layer that predicts either absence or presence of the species for each pixel in the layer. An automatic “best subsets” option in GARP may be implemented to select the best model set from all models (see Anderson, Lew and Peterson (2003) for discussion of the best subsets selection procedures). Typically, a 10-best model set is selected from all GARP models (Anderson, Lew and Peterson 2003).
GARP has yielded successful predictive results in a variety of ecology, biodiversity, and related applications. For example, Peterson, Ball and Cohoon (2002) modeled ecological niches for twenty-five bird species in Mexico with GARP. Peterson et al. (2001) projected ecological niche models developed in GARP to predict the potential effects of global climate change on the distribution of the Cracidae bird species in Mexico. In addition to terrestrial species, GARP has been extended to aquatic environments for ecological modeling of fish species in both freshwater (McNyset 2005) and marine (Wiley et al. 2003) environments.
Generally, BTPDs prefer sites with flat or gentle slopes, deep soils with a low sand content, and vegetation characteristic of the shortgrass prairie. Environmental variables selected for ecological niche modeling were Normalized Difference Vegetation Index (NDVI) data layers for the spring, summer, and fall growing seasons (Fig. 1); slope (Fig. 2A); soil depth (Fig. 2B); and soil texture (Fig. 2C). Environmental variables were selected based on known characteristics of BTPD habitats, in addition to consultation with Kansas BTPD experts and a literature review of related studies for other regions of the Great Plains. Proctor, Beltz and Haskins (1998) used similar variables (vegetation, slope, soil depth, and soil texture) to model BTPD habitats in the northern Great Plains with GIS, and through a classification tree analysis and logistic regression found correlations between three of these variables (vegetation, slope, and soil texture) with BTPD habitat presence.
NDVI, a measure of the amount and condition of biomass and vegetation, was selected to serve as a surrogate for land cover in the ecological niche modeling. Although the 250 m spatial resolution of the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite sensor is coarser than the 30 m resolution of the Landsat Thematic Mapper (TM) sensor, a benefit of MODIS is improved temporal resolution. Other land cover mapping studies (e.g., Wardlow 2006; Wessels et al. 2004) have noted the advantages of the MODIS temporal resolution for discriminating among land cover types throughout the growing season. NDVI data layers were derived from biweekly MODIS satellite imagery composites for the year 2001 from March through October (fifteen dates in total) to encompass typical vegetation of the growing season in Kansas. NDVI was calculated for each of the fifteen MODIS images using the following equation: [(Near-Infrared Band − Red Band) / (Near-Infrared Band + Red Band)].
Soil depth and soil texture layers were derived from the Soil Survey Geographic (SSURGO) database, a detailed digital soil survey (1:24,000 scale) compiled on a county basis by the Natural Resources Conservation Service of the U.S. Department of Agriculture. Attribute data from the Soils Lite database for Kansas, a subset of commonly used SSURGO soil attributes, were joined to the SSURGO map unit polygons. The Soils Lite database includes measures of the average soil depth (in) to bedrock and an ordinal classification (1–21) of soil surface texture ranging from clay (1) to coarse sand (21) for each SSURGO map unit polygon. The SSURGO and Soils Lite databases were available for all Kansas counties in the study area.
Slope (%) was calculated from 30 m digital elevation models (DEMs) available from the U.S. Geological Survey's (USGS) National Elevation Dataset (NED). Striping, rows of cells with identical elevation values, and other data artifacts were present in the NED dataset in portions of Dickinson, Morton, Saline, and Sedgwick Counties. Corrected 10 m DEMs, available for each of the USGS topographic quadrangles with striping problems, were resampled to 30 m and replaced the quadrangles with striping in the NED dataset.
Species Occurrence Data
Species occurrence point data were selected from a database of known BTPD colonies provided by the KDWP. Latitude and longitude coordinates for each colony were visually estimated by observers in airplanes flying along north-south oriented transects across western Kansas. Although an advantage of the data collection method is a large sample size with coverage for much of the western half of the state, one limitation of the aerial survey data is the positional accuracy of each estimated coordinate. For this reason, the Kansas GAP Analysis Project land cover map (Egbert et al. 2001) was used to eliminate any species occurrence points that corresponded to land cover classes unsuitable for BTPD habitation (cropland, Conservation Reserve Program (CRP), forest/woodlands, marshland, urban, and water) as a result of positional error introduced during the coordinate estimation procedures. The remaining species occurrence points (455 total points) were randomly subset into training (80%, 364 points) and validation (20%, 91 points) datasets (Fig. 3). The training dataset was input into GARP for training and internal testing during the modeling process, whereas the validation dataset was withheld entirely for an independent validation of model results through Receiver Operating Characteristic (ROC) analysis. An independent validation data subset is preferable to data re-substitution for assessing model quality (Fielding and Bell 1997). Given the number of occurrence points available, an 80/20 training/testing subset was used as this satisfies the minimum required number of training points necessary to achieve maximal GARP model accuracy (Stockwell and Peterson 2002), while maintaining a testing subset of sufficient size to achieve statistical power in the ROC analysis (Hanley and McNeil 1982).
Desktop GARP requires that all environmental layers be in raster data format at the same spatial resolution or pixel size. Therefore, the slope layer was resampled from 30 m to 250 m and the soil texture and soil depth to bedrock layers were converted from vector to raster format and sampled to 250 m to match the native spatial resolution of the NDVI layers derived from the MODIS satellite imagery.
A mask layer was created to exclude all geographic areas in Kansas east of the Flint Hills (i.e., beyond the historic range of the BTPD) in the resulting raster layers from the GARP modeling. In addition, cropland, CRP, forest/woodlands, marshland, urban, and water land cover classes were extracted from the Kansas GAP land cover map and added to the mask layer in order to exclude all land cover types unsuitable for BTPD habitation. CRP was included in the mask since CRP land is generally considered unsuitable for BTPD habitation (USDA Forest Service 2005).
In order to achieve more regional variation in the resulting ecological niche models, the environmental variables and species occurrence datasets for the state were divided into two zones, one north and one south of 38° 42′N latitude. Preliminary results by Stockwell and Peterson (2002) have demonstrated that such regional subdivision may improve GARP model results. A north-south division of the state was preferable to an east-west split in order to ensure an equivalent number of species occurrence data points in each zone.
Jack-knifing was performed in Desktop GARP separately for the north and south zones to determine if any environmental layers should be excluded from the ecological niche modeling due to high correlations with omission, i.e., an error corresponding to an absence prediction where species occurrence data actually occur, or a “false negative”. Environmental layers with a correlation (r) to omission of greater than 0.10 were removed from the ecological niche modeling. One NDVI layer (October 16) for the south zone (r = 0.11), and one NDVI layer (September 14) for the north zone (r = 0.13) were removed due to correlations with omission that exceeded 0.10.
GARP was run separately for the north and south zones with a convergence level of 0.001 and 10000 maximum iterations for both runs. The “best subsets” option in Desktop GARP was enabled (omission threshold = 5%, commission threshold = 50% of distribution, 20 total models under hard omission threshold) to select automatically the 10-best model set for each of the two zones. The 10-best model sets for both zones were merged together to create the final BTPD habitat suitability map.
Results and Discussion
Values on the resulting BTPD habitat suitability map (Fig. 4) range from 0 (0 of 10 models predicted presence or the area was masked during GARP modeling) to 10 (10 of 10 models predicted presence). Areas with a high density of BTPD colonies, such as Logan County in western Kansas and portions of the Cimarron National Grassland in Morton County in the southwest corner of the state, are characterized by the intersection of all models predicting presence (values of 10). Portions of central Kansas west of the Flint Hills, and other areas beyond the range of the BTPD are generally characterized by values of 1–9, indicating variation between the model predictions, or values of 0 for the intersection of all models predicting absence.
Model Validation with Receiver Operating Characteristics (ROC) Analysis and Average Omission Error
Two methods of model validation commonly used in ecological niche modeling, Receiver Operating Characteristics (ROC) analysis (e.g., McNyset 2005; Wiley et al. 2003) and average omission error calculation (e.g., McNyset 2005), were used to assess the overall accuracy of the 10-best model set. ROC analysis compares model results with the species occurrence validation data withheld from the GARP modeling. The area under the curve (AUC) statistic is a measure of the overall accuracy of the model set; an AUC statistic of 1.0 indicates total agreement between the model set and the species occurrence test data whereas a model set developed randomly yields an AUC of 0.5. Average omission error is a measure of the average omission for the 10-best models. The BTPD model set for Kansas yielded an AUC = 0.701 and an average omission = 7%, indicating good overall model set predictions.
Model Validation with an Additional Independent Dataset: An additional dataset of thirty-seven active BTPD colonies in the Cimarron National Grassland in Morton County was available for an independent validation of the 10-best model set (Fig. 5). These colonies were on publicly held portions of the Cimarron National Grassland, and were identified by visually and audibly locating BTPDs in the area. Perimeters of these BTPD colonies were mapped as polygon features with a differential Global Positioning System (GPS) receiver and then post-processed with data from a base station in Elkhart, KS. The use of differential GPS for mapping these colonies resulted in locations that were considerably more accurate than the statewide BTPD species occurrence dataset.
The centroid of each BTPD colony polygon was computed and used to calculate the AUC and average omission for Morton County (AUC = 0.628, average omission = 6%). In addition, the 10-best model set values for each pixel in each of the polygons were summed and converted to percentages of total area for each BTPD colony (Table 1). The results indicate a high degree of accuracy in the model predictions. For example, each of the 10-best models predicted presence (i.e., values of 10) for at least 75% of the total area of the colony for thirty of the thirty-seven (81.1%) colonies. Likewise, each of the 10-best models predicted presence for 100% of the total colony area for seventeen of the thirty-seven (45.9%) colonies.
Calculated area (in square meters) of masked and 10-best model values for each BTPD colony in Morton County, Kansas as well as percentage of total colony area.
Interestingly, seventeen of the thirty-seven (45.9%) colonies fell partially within the masked land use/land cover types with percentages ranging from 2.4% to 99.7% of the total colony area. To evaluate this artifact, colonies were overlaid on top of the GAP land cover map used to generate the mask. Perhaps since the landscape is highly fragmented in western Kansas, several BTPD colony boundaries slightly extend beyond prairie land cover types into cropland or other land use/land cover types that are deemed unsuitable for BTPD habitat. These particular colonies contain relatively small percentages of the masked area.
Colonies with higher percentages of the masked area typically fell on land classified as CRP by the GAP land cover map. While land enrolled in CRP is planted to grassland species, it is generally considered to be unsuitable for BTPD habitat (USDA Forest Service 2005). However, it is possible that BTPDs have successfully established colonies on CRP in western Kansas and would suggest that CRP should not have been masked a priori to the GARP modeling. Perhaps a more likely explanation is that native prairie in these areas were misclassified as CRP in the GAP land cover map used for the mask. While general land cover classes (herbaceous, scrubland, forest, cropland, urban, and water) were mapped at an overall accuracy of 88%, accuracy levels for specific grassland alliance-level communities were substantially lower (Egbert et al. 2001). Field verification of the GAP land cover types that correspond to BTPD colonies in Morton County is necessary to verify these possible explanations.
Conclusions and Future Work
Maps of BTPD habitat suitability resulting from ecological niche modeling in GARP have the potential to help state agencies and organizations, including the Kansas Black-Tailed Prairie Dog Working Group and the KDWP, in their efforts to protect the species in the state through implementation of the “Kansas Black-Tailed Prairie Dog Conservation & Management Plan.” This paper has presented one approach for modeling BTPD habitat suitability through GIS, remote sensing, and ecological niche modeling with GARP, a powerful modeling tool that produces robust results with relative ease of use. This project also indicates the potential benefits of integrating high temporal resolution satellite imagery from MODIS into ecological niche modeling for deciphering land cover types. An avenue of future research for this project is a direct comparison of GARP models derived from high temporal resolution satellite imagery such as MODIS to those derived from the higher spatial resolution but lower temporal resolution of Landsat satellite imagery to determine the effect of spatial and temporal resolution on model accuracy. In addition, future research will also examine the effect of higher resolution GIS datasets (e.g., 10 m DEMs) and the improved accuracy of species occurrence data collected with GPS on resulting GARP models. Field verification of model results is an additional area of future research that will be beneficial for determining the effect of spatial and temporal resolution of datasets on the resulting ecological niche models.
Funding for this project was provided by the AmericaView Program (grant AV04-KS01) through the U.S. Geological Survey. We thank the following agencies and individuals for project assistance: the Kansas Department of Wildlife and Parks; Mike Houts and Bill Busby of the Kansas Biological Survey; Brian Wardlow of the Kansas Applied Remote Sensing Program; Rich Sleezer of Emporia State University; A. Townsend Peterson of the Natural History Museum & Biodiversity Research Center and Department of Ecology and Evolutionary Biology at the University of Kansas; Ricardo Scachetti-Pereira of the University of Kansas; and John Dunham of the Kansas Data Access and Support Center. We thank the three anonymous reviewers for their comments and suggestions that improved the final manuscript.