Linking GIS with Models
This paper was presented at the Third International Conference on Integrating GIS and Environmental Modeling, at Santa Fe, New Mexico, on January 23, 1996.
Linking GIS with Models of Ecological Risk Assessment for Endangered Species
Applied Biomathematics
100 North Country Road
Setauket, NY 11733, USA
Abstract
A model that links GIS to models for viability analysis and risk assessment is applied to endangered species, including the Spotted Owl in the northwestern US, the Red-cockaded Woodpecker in Louisiana, and the California Gnatcatcher in Orange County, California. The model integrates landscape data on habitat requirements with demographic data to analyze risks of extinction, evaluate management options, and assess human impact on wildlife populations. Other applications of the model involve design of nature reserves, wildlife management, and population viability analysis. The model analyzes habitat data exported from a GIS, and identifies the patches of habitat that can support a population. The structure of these patches, including their locations, sizes and distances from each other, define the spatial structure of the metapopulation. The spatial structure is combined with demographic data and other information on the ecology of the species to complete a metapopulation model, which incorporates age or stage structure and density dependence for each population, spatial correlation and dispersal among populations, environmental and demographic stochasticity and catastrophes. The model performs a risk analysis, and runs multiple simulations, automatically changing parameters to analyze the sensitivity of risks to input data.
Introduction
One of the main tasks of conservation biologists is to evaluate the viability of endangered and threatened species under different natural conditions, and under alternative options for wildlife management, reserve design, and habitat protection plans. These evaluations usually ask questions about the predicted future abundance, risk of extinction, or chance of recovery of the species; and are addressed by population viability analysis (PVA), which is a systematic examination of interacting factors that place a population or species at risk of extinction (Gilpin and Soul‚ 1986, Shaffer 1990). The factors that a PVA examines may be both natural and anthropogenic in origin, and their analysis often involves mathematical or computer models that predict the future changes in the abundance and distribution of the species in question, given information about its ecology and demography (Burgman et al. 1993).
Habitat loss and fragmentation are among the most common threats facing endangered species, making GIS-based evaluations an essential component of population viability analyses. Often habitat loss and fragmentation, combined with the natural heterogeneity of landscapes, forces species to exist in multiple populations inhabiting relatively isolated habitat patches. Such a collection of populations of the same species is called a metapopulation. The existence of multiple populations usually introduce complexities that make it impossible to evaluate the viability of the species based on PVAs performed on separate populations.
This paper describes a computer program for building metapopulation models and performing GIS-based PVAs, and discusses its application to cases involving endangered species.
The Model
The PVA program RAMAS GIS is designed to link GIS-generated landscape data with a detailed metapopulation model for extinction risk assessment, viability analysis, reserve design and wildlife management. Descriptions of the first version of the program (Akçakaya 1994) can be found in Akçakaya (1995) and Akçakaya et al. (1995a); it was also reviewed by Kingston (1995). This section describes the second version of the program, which is being tested with applications to various endangered species (see the section on "Applications").
The program operates in four steps. First, landscape data is analyzed and the patch structure is exported to a metapopulation model. Second, temporal changes in habitat characteristics are modeled. Third, a metapopulation model is built by combining spatial and demographic information. Fourth, simulations are run to estimate risks of extinction or decline, and to predict the abundance and distribution of individuals in the metapopulation. Below the essential aspects of these four model components are summarized.
Landscape data
The function of the "Landscape data" component is to analyze GIS data to determine the spatial structure of the metapopulation as well as several population-specific demographic parameters that may depend on habitat quality. This component works in 6 steps described below.
1. Imports landscape data from a GIS
RAMAS GIS operates under MS-DOS [note: current version requires Windows], and can import landscape data in the form of raster maps exported in ASCII format from ARC/INFO, GRASS, IDRISI (and any GIS that supports bare grid format), and can also import in binary format from IDRISI. The requirements are: (i) maps must be in raster format, with the "cells" (i.e., locations given a value by the map) must be arranged in a square (not hexagonal) format; (ii) map values must be numerical; (iii) all maps must describe the landscape with the same precision (i.e., number of cells in both north-south and east-west directions must be the same); (iv) coverage of all maps must be identical (the corners of the rectangular area described by all maps must be the same); (v) maximum map size depends on the available memory, but cannot exceed 5000 rows and 5000 columns in version 2.0.
2. Creates habitat map
Information in different maps are then combined to make a habitat suitability (HS) map, with a user-defined habitat suitability function. This function must be estimated outside the program. Most methods of estimating the HS function involve statistical procedures, using species occurrence or abundance at each location as the dependent variable and the habitat characteristics as the set of independent variables. The statistical procedures most commonly used are stepwise multiple regression methods (including logistic regression) and stepwise discriminant function analysis. The function is then entered in RAMAS/GIS, which creates a habitat suitability map by calculating the value of the function for each cell of the raster map.
3. Finds habitat patches
The program then employs a patch-recognition algorithm to find clusters or groups of nearby cells that have HS values higher than or equal to a threshold habitat value, and labels them as patches. For this procedure, the program uses two parameters that determine how the species perceives (or reacts to) the patchiness of the habitat. These are threshold HS value and neighborhood distance. Threshold HS is the minimum HS value (as defined by the HS function) below which the habitat is not suitable for reproduction and/or survival (although individuals may disperse or migrate through habitat that has a lower HS than this threshold). Neighborhood distance is used to identify nearby cells that belong to the same patch. Suitable cells (as defined by the threshold parameter) that are separated by a distance less than or equal to the neighborhood distance are regarded to be in the same patch. For an animal species, the neighborhood distance parameter may represent the foraging distance.
4. Calculates demographic parameters
Next, the program calculates five demographic parameters for each patch. These are carrying capacity, initial abundance, maximum growth rate, relative fecundity and relative survival. The meaning of these parameters will be described below (see section "Metapopulation model"). The program allows these five parameters to be calculated as any arbitrary function of patch characteristics, such as total habitat suitability in the patch, average habitat suitability, area as the number of cells in the patch, length of the perimeter (edge) of the patch, average values of input maps for the patch. In addition to these landscape variables, the functions for the five demographic parameters can also include any of the standard mathematical functions (such as min, max, log, exp, etc.)
5. Calculates the spatial structure
In addition to the population-specific demographic parameters described above, the program calculates several parameters related to the spatial structure of the metapopulation. These include the location of patches, the distances among patches and migration (dispersal) rates, and spatial correlations based on these distances. These two metapopulation-level parameters, dispersal and correlation, will be discussed below. The rate of dispersal between two population can depend on one of three distance measures: edge to edge (minimum distance between the two patches), center to edge (distance from the center of source patch to the closest edge of the target patch), or center to center (distance between the centers of the two patches).
6. Exports the results to a metapopulation model
Finally, the program exports the patch coordinates, demographic parameters (carrying capacities, growth rates, etc.), the dispersal (migration) rates and the correlation structure to the metapopulation model.
Habitat dynamics
The second component of the program is designed to be used in modeling temporal changes in habitat. It allows the calculation of a time-series of carrying capacities for each population. It reads data files saved from the Landscape data subprogram, and outputs a set of data files for the Metapopulation model subprogram. One of these data files contains the main metapopulation model, and others contain the temporal changes in carrying capacities.
The Metapopulation model subprogram inputs these files of carrying capacities and vital rates, and uses them to model habitat dynamics.
Metapopulation model
The main component of the program is where the spatial and demographic parameters calculated by the Landscape data component are combined with other ecological and demographic information about the species to develop a metapopulation model. The model may incorporate the following factors and parameters discussed below.
Age structure or Stage structure within populations is modeled by a matrix model (Caswell 1989) that incorporates age- or stage-specific vital rates (survival rates and fecundities). Each population in the model can have a different stage matrix, and a different initial age or stage structure (initial number of individuals in each age or stage). The initial structures can be specified as the stable age or stage distributions. The population-specific stage matrix can be specified to change through time, by reading two files (one each for fecundities and survival rates) that contain the temporal change in the relative values of the vital rates.
Density dependence in population dynamics is modeled by modifying the mean values of survival rates and fecundities as a function of the population size (N). Density-dependent population growth may involve a simple ceiling model, logistic-like functions that describe contest- or scramble-type intraspecific competition (including Ricker and Beverton-Holt functions), Allee effects (i.e., density dependence at low population sizes), or Allee effects combined with density dependence at high population sizes. All density dependence functions are parameterized with the same set of parameters that include maximal growth rate (Rmax) and carrying capacity or ceiling (K), random variation in K, and temporal trend in K. Each population can have a different set of parameters.
Habitat change, e.g., habitat loss (as a result of human impact) or habitat increase (as a result of vegetation growth) can be modeled by specifying how the carrying capacity of the population changes through time. This can be done either with a constant rate, or as a time-series of carrying capacities saved in a disk file.
Environmental stochasticity is modeled by random fluctuations in vital rates and in carrying capacities. The random fluctuations can be normal- or lognormal-distributed, and can be correlated among populations. They are assumed to be perfectly correlated among age classes or stages within each population.
Demographic stochasticity is modeled by sampling the number of survivors from a binomial and the number of offspring from a Poisson distribution (Akcakaya 1991).
Catastrophes can affect abundances (a proportion of all individuals die), vital rates (survival rates and fecundities are reduced after a catastrophe), or carrying capacities (which are reduced after a catastrophe). The spatial extent of catastrophes may be local or regional. The impact of catastrophes can be population-specific (some populations may be more prone to, or more affected by catastrophes than other populations), or they may be stage-specific (some stages, or even certain vital rates may be more affected by catastrophes than others).
Density dependence in migration is modeled by making the total rate of emigration or dispersal a function of the size of the population. The rate of dispersal/emigration can be specified either to increase or decrease as N increases. Migration rates can also be stage-specific.
Geographic configuration is specified by the coordinates of each population. The distance of populations from each other and their relative positions are utilized in modeling the effects of the two spatial factors discussed below.
Dispersal (migration) describes the movement of individuals among populations. In RAMAS GIS, migration is modeled by specifying the proportion of individuals that move from each population to each other at every time step. These rates are input in the form of a migration matrix. In most cases, the rate of dispersal may be a function of the distance between source and target populations. RAMAS GIS allows users to specify a function that describes the dependence of dispersal rates on distance. The matrix can filled according this to function, and can be edited to account for habitat corridors (by increasing the rate between specific pairs of populations) and for obstacles or geographic barriers to migration (by decreasing the rate). The migration rates may also be specified to be dependent on population size to allow for density-dependent migration (see above), or on the age or stage of the individuals to allow for age- or stage-specific dispersal tendencies.
Correlations among populations describe the similarity of environmental patterns experienced by each population. This factor is important in the "rescue effect" in metapopulations: when fluctuations are spread over a number of separate populations, the overall risk faced by the metapopulation is reduced. If the fluctuations in the environment are at least partially independent (uncorrelated), it will be less likely that all populations go extinct at the same time than if the fluctuations are dependent (i.e., synchronous or correlated). In uncorrelated environments, extinct populations will have a chance to be recolonized. Conversely, if the fluctuations are partially synchronous (correlated), then models based on an assumption of independent population dynamics among patches will underestimate extinction risk. Thus, correlation of environmental fluctuations has important effects on metapopulation persistence and viability (Gilpin 1988; Akçakaya and Ginzburg 1991; Burgman et al. 1993).
RAMAS GIS models correlations by sampling the vital rates of each population from a normal or lognormal distribution which is correlated with the vital rates of other populations according to a correlation matrix specified by the user. This matrix can be specified as a function of the distance between populations (as closer populations are more likely to experience similar environmental patterns).
Risk assessment
In this component, metapopulation models built in RAMAS are used to simulate the dynamics and to predict the future of the metapopulation. RAMAS GIS summarizes the results of a simulation with 16 types of output, some of which are superimposed. Most of these results are related to risk assessment and report risk analytical measures such as risk of extinction and time to extinction (Akçakaya 1992). Types of output include
Risk of extinction (or decline to any level) any time during, or at the end of, the simulated time period
Probability of exceeding a range of population sizes any time during, or at the end of, the simulated time period
Distribution of times to extinction (or decline to a specific level), and cumulative time to extinction or decline
Distribution of times that the metapopulation will exceed a specified threshold, and cumulative time to increase
Abundance (of the metapopulation and each of its populations) through time
Metapopulation occupancy (number of extant populations through time)
Local occupancy (for each population, the number of time steps that it was extant)
Local extinction duration (for each population, the maximum number of consecutive time steps during which a population remained unoccupied (below local threshold defined by user) during a typical replication
Final stage structure (histogram of the number of individuals in each stage of a specific population at the end of the simulated time period)
Metapopulation structure at a specific time step (histogram of the number of individuals in each population at a specific time).
The model can be run several times, to analyze the sensitivity of results to input parameters by varying them automatically, to compare management options, or to assess anthropogenic impact by comparing outputs from simulations with parameters for impacted and non-impacted situations. The sensitivity analysis facility also allows the user to superimpose graphs from different simulations to make the comparisons easier.
Applications of RAMAS GIS, the second version of the program, is currently being tested with applications to four endangered birds (California gnatcatcher, cactus wren, red-cockaded woodpecker, and northern spotted owl). These applications are briefly described below.
California gnatcatcher and cactus wren
The program is applied to California gnatcatcher Polioptila californica and cactus wren Campylorhynchus brunneicapillusmetapopulations in central and coastal regions of Orange County, California. The California gnatcatcher (a federally-listed threatened species) and cactus wren are two of the three "target species" in the State of California's Natural Community Conservation Planning (NCCP) program. The target species are considered as surrogates for the conservation of coastal sage scrub, a plant community that has decreased substantially from its historic coverage as a result of urban and agricultural development, and has become fragmented.
The application started with a statistical analyses of observation locations, and GIS maps for elevation, slope, and vegetation (exported from ARC/Info) to characterize the habitat of the two birds. The habitat maps were then analyzed to determine the spatial structure of the metapopulation model. These data were combined with demographic parameters estimated from field work conducted by Atwood et al. (1995) to produce spatially-explicit, stage-structured, stochastic metapopulation models for the two species(Akçakaya and Atwood, in prep.). Simulations and sensitivity analyses with these models pointed out to the relative importance of temporal variation in vital rates (including the frequency and the magnitude of catastrophic declines in vital rates as was observed in the winter of 1994-95).
Red-cockaded woodpecker
The goal of this application of the program was to evaluate the impact of timber prescriptions (forest management practices) on the viability of a red-cockaded woodpecker Picoides borealis metapopulation in Louisiana. The habitat of the red-cockaded woodpecker was characterized based on the stand type (i.e, the dominant tree species), stand condition, and basal area of pine species. These variables were input into the program in the form of GIS maps exported from GRASS (Akçakaya et al. 1995b). Currently the habitat characterization is being extended to include other habitat variables such as distance to streams, and age of dominant trees in the stand. The habitat descriptions are used to evaluate the impact of forest management practices, such as a planned cut, on the carrying capacities of habitat patches and other parameters of the red-cockaded woodpecker populations inhabiting these patches. The model will then be run with different management options to compare their impacts in terms of the predicted risk of extinction, or rate of recovery of the metapopulation.
Northern spotted owl
This application focused on factors affecting the viability of the Northern Spotted Owl Strix occidentalis caurina throughout its range in the U.S. A second goal in this analysis was to incorporate two sources of variability in determining the threat the species faces. The study incorporated natural variation (resulting from temporal fluctuations in environmental factors) in the form of randomly distributed vital rates (survivals and fecundities). In addition, demographic stochasticity was modeled to describe chance variations in reproduction, survival and dispersal. These types of natural variation (environmental and demographic) were used to express the model results in probabilistic terms such as the viability of the species (for example in terms of the chance of survival or risk of extinction). Uncertainties that result from a lack of knowledge were incorporated in the form of parameter ranges, and were used to estimate upper and lower bounds on the estimated viability of the species.
Based on the habitat maps provided by the Forest Service, the program found 18 habitat patches. The size distribution of the patches was very skewed, with the 4 largest patches making up about 95% of the total area of all patches, and the seven largest making up about 97%. Because of the large differences in sizes of neighboring populations, the model results (risk of decline) were not very sensitive to the rate of inter-patch dispersal of juvenile spotted owls. The model predicted a large difference between lower and upper bounds on the viability of the northern spotted owl, based on the best-case and worst-case scenarios which were parameter combinations that resulted in best and worst chance for survival. According to sensitivity analyses, the viability of the metapopulation was most sensitive to the set of vital rates used (the dependence of fecundities and survival rates on habitat), and also sensitive to the degree of spatial correlation among vital rates of the populations, and to the carrying capacities of the populations. In addition, metapopulation occupancy was sensitive to dispersal and Allee effects (Akçakaya, unpublished).
Other applications
The methodology used in the models, and the relative effects of factors discussed above on metapopulation dynamics have been discussed in several previous publications (Akçakaya 1991, 1992; Akçakaya and Ginzburg 1991; Burgman et al. 1993). Previous versions and precursors of the program have been applied to a wide variety of organisms. RAMAS GIS (version 1.0) was used in assessing the effectiveness of translocating individuals as a management option for the endangered helmeted honeyeater (Lichenostomus melanops cassidix) in Australia (Akçakaya et al. 1995a).
The metapopulation model component was used to build a model for a land snailArianta arbustorum metapopulation in northeastern Switzerland. The model was used to investigate the effect of population subdivision on the persistence of a land snail metapopulation and to analyze the interaction between spatial factors, population subdivision, and catastrophes (Akçakaya and Baur, in press).
A precursor to the metapopulation model component (RAMAS Space; Akçakaya and Ferson 1990) has been used to model metapopulation dynamics the California spotted owl (LaHaye et al. 1994). California spotted owl (Strix o. occidentalis) is a subspecies of the endangered Northern Spotted Owl (S. o. caurina) found in the Pacific Northwest. Gibbs (1993) used RAMAS Space to model dynamics of spotted salamanders, bullfrogs, snapping turtles, swamp sparrows and water shrews in a mosaic of wetlands.
References
Akçakaya HR. 1991. A method for simulating demographic stochasticity. Ecological Modelling 54:133-136.
Akçakaya HR, McCullough DR, Barrett RH.1992. Population viability analysis and risk assessment. Wildlife 2001: Populations:148-157. Elsevier Publishers, London.
Akçakaya HR. 1994. RAMAS/GIS: Linking Landscape Data With Population Viability Analysis (version 1.0). Applied Biomathematics, Setauket, New York.
Akçakaya HR. 1994b. GIS enhances endangered species conservation efforts. GIS WORLD Vol.7, November 1994:36-40.
Akçakaya HR, Baur B. Effects of population subdivision and catastrophes on the persistence of a land snail metapopulation. Oecologia.
Akçakaya HR, Ginzburg LR. 1991. Ecological risk analysis for single and multiple populations. Species Conservation: A Population-Biological Approach: 73-87. Birkhaeuser Verlag, Basel.
Akçakaya HR, Ferson S. 1990. RAMAS/space User Manual: Spatially Structured Population Models for Conservation Biology. Applied Biomathematics, New York.
Akçakaya HR, McCarthy MA, Pearce J. 1995. Linking Landscape Data with Population Viability Analysis: Management Options for the Helmeted Honeyeater. Biological Conservation 73:169-176.
Akçakaya HR, Jackson J, Ginzburg LR. 1995b. Modeling the Red-cockaded Woodpecker Populations at Fort Polk. Report to U.S. Army Corps of Engineers.
Atwood JL, Fugagli MR, Reynolds CH, Luttrell JC, Tsai SH. 1995. Distribution, population dynamics, and dispersal of California Gnatcatchers on the Palos Verdes Peninsula, 1993-1995. Proceedings of CalGnat'95 Symposium.
Burgman MA, Ferson S, Akçakaya HR. 1993. Risk Assessment in Conservation Biology. Chapman and Hall, Population and Community Biology Series.
Caswell H. 1989. Matrix Population Models: Construction, Analysis, and Interpretation. Sinauer Associates, Sunderland, Massachusetts.
Gibbs JP. 1993. Importance of small wetlands for the persistence of local populations of wetland associated animals. Wetlands 13:25-31.
Gilpin ME. 1988. A comment on Quinn and Hastings: extinction in subdivided habitats. Conservation Biology 2:290-292.
Gilpin ME. and M.E. Soule. 1986. Minimum viable populations: processes of species extinction. Conservation Biology: the science of scarcity and diversity:19-34. Sunderland, Massachusets.
Kingston T. 1995. Valuable modeling tool: RAMAS/GIS. Conservation Biology 9:966-968.
LaHaye WS, R.J. Gutierrez, and H.R. Akcakaya. 1994. Spotted owl metapopulation dynamics in southern California. Journal of Animal Ecology 63:775-785.
Shaffer ML. 1990. Population viability analysis. Conservation Biology 4:39-40.