A cellular automata model of land cover change to integrate urban growth with open space conservation

Landscape and Urban Planning 99 (2011) 141–153

Contents lists available at ScienceDirect

Landscape and Urban Planning journal homepage: www.elsevier.com/locate/landurbplan

A cellular automata model of land cover change to integrate urban growth with open space conservation Diana Mitsova a,∗ , William Shuster b , Xinhao Wang c a

Florida Atlantic University, School of Urban and Regional Planning, Fort Lauderdale, FL 33301, United States U.S. Environmental Protection Agency, National Risk Management Research Laboratory, Sustainable Technology Division/Sustainable Environments Branch, Cincinnati, OH 45268, United States c University of Cincinnati, School of Planning, 6210 DAAP, Cincinnati, OH 45221, United States b

a r t i c l e

i n f o

Article history: Received 13 May 2010 Received in revised form 1 October 2010 Accepted 18 October 2010

Keywords: Cellular automata Urban growth Green infrastructure Markov transition probabilities Multi-criteria evaluation Landscape metrics

a b s t r a c t The preservation of riparian zones and other environmentally sensitive areas has long been recognized as one of the most cost-effective methods of managing stormwater and providing a broad range of ecosystem services. In this research, a cellular automata (CA)—Markov chain model of land cover change was developed to integrate protection of environmentally sensitive areas into urban growth projections at a regional scale. The baseline scenario is a continuation of the current trends and involves only limited constraints on development. The green infrastructure (GI) conservation scenario incorporates an open space conservation network based on the functional boundaries of environmentally sensitive areas. It includes variable buffer widths for impaired streams (as identified on the USEPA 303d list for stream impairment), 100-year floodplain, wetlands, urban open space and steep slopes. Comparative analysis of each scenario with landscape metrics indicated that under the GI conservation scenario, the number of urban patches decreased while the extent of interspersion of urban land with green infrastructure patches increased leading to improved connectivity among open space features. The analysis provides a quantitative illustration of how our process contributes towards achieving urban planning objectives while incorporating green infrastructure. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Landscape alterations due to urban development have a profound effect on ecosystem functions, and can regulate or modulate the benefits that humans derive from physical, biological or chemical properties and/ or processes occurring in natural systems (Costanza et al., 1997). These benefits, often termed ecosystem services, include infiltration and evapotranspiration of stormwater runoff, groundwater recharge, protection of water resources against sedimentation, reduced risk of nutrient enrichment and contamination, and preservation of landscapes that provide aesthetic and recreational value (Osborne and Kovacic, 1993; Dosskey, 2001; Randolph, 2004). A wider vegetated riparian buffer zone can lead to a longer travel time and distance for runoff, increasing the number of infiltration opportunities, and thereby promoting deposition of eroded soil material, facilitating nutrient removal while also cooling stormwater before it reaches the streambed (Osborne

∗ Corresponding author at: Florida Atlantic University, School of Urban and Regional Planning, 111 East Las Olas Blvd. HEC 1009J, Fort Lauderdale, FL 33301, United States. Tel.: +1 954 762 5674; fax: +1 954 762 5673. E-mail addresses: [email protected] (D. Mitsova), [email protected] (W. Shuster), [email protected] (X. Wang). 0169-2046/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.landurbplan.2010.10.001

and Kovacic, 1993; Dosskey, 2001). Moreover, terrain characteristics such as slope, vegetative cover, infiltration rates, soil moisture storage capacity, and slope can impact the effectiveness of a riparian buffer zone for water quantity and quality objectives (Polyakov et al., 2005). However, the importance of preserving and enhancing natural areas in an urbanizing landscape has often been understood after development has removed or altered vegetated areas. A key principle of smart growth which attempts to minimize the negative side effects of urbanization through a set of planning and policy options, is preservation of open space, prime agricultural land and environmentally sensitive areas (Smart Growth Network, 1996). The Smart Growth Manual specifically emphasizes the role of regional planning in prioritizing areas for development beginning with the creation of a green footprint, i.e., mapping a region’s “rural reserve” and other environmentally important areas (Duany et al., 2009). The concept of green infrastructure describes the interdependence of land conservation and land development (Benedict and McMahon, 2006) and refers to a contiguous, interconnected green network consisting of riparian areas, floodplains, aquifer recharge zones, wetlands, and forested areas. To incorporate green infrastructure into land management states of Florida, Georgia and Maryland have adopted programs aimed at enhancing connectivity of protected natural areas as part of a state-scale initiative to establish and maintain green infrastructure networks (Randolph, 2004;

142

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

Benedict and McMahon, 2006). The Trust for Public Land works in collaboration with local governments to implement a Greenprinting program which envisions a symbiosis of future growth and conservation principles. A major component of the program is acquiring conservation lands on which communities depend for drinking water supply, stormwater management, recreation, and agricultural production (TPL, 2005). There are planning opportunities to preserve, integrate, and otherwise connect green infrastructure as a part of a comprehensive regional conservation plan for landscapes anticipating development pressure. A model for expediting this type of analysis and forecasting involves the use of cellular automata, which are wellsuited to represent geographic processes due to the similarities between a two-dimensional lattice and a raster grid (Longley and Batty, 1996; Clarke and Gaydos, 1998; Torrens, 2003; Batty and Xie, 2005). Torrens (2003) suggests that cellular automata can be useful in simulating urban systems because land use, concentrations of employment location and population change can all be modeled as automata; cells can effectively aggregate economic, demographic and transportation data; neighborhoods as part of the cityscape can be successfully simulated by equivalent neighborhoods of cells on the cellular lattice; and even urban theories based on spatial interaction models can be tested and represented. Claggett et al. (2004) applied the results from SLEUTH (Clarke et al., 1997) to the assessment of future development pressure on forested and agricultural lands for the purposes of vulnerability assessment of productive lands at a regional scale. Syphard et al. (2005) used Clarke’s Urban Growth Model (UGM) to estimate possible impact of urbanization around Santa Monica Mountain National Recreation Area through metrics such as total core habitat area, mean core patch area and number of patches, and concluded that important migration corridors will be affected if development of urban clusters does not take into consideration the spatial requirements of wildlife habitat (Syphard et al., 2005). Eppink et al. (2004) coupled a dynamic model of population growth and projected that draining wetlands can result in significant biodiversity loss and a cutback on such activities may counter to a certain extent the impact of rapid and sustained urbanization. In this research, we develop a cellular automata model of urban growth to incorporate environmentally sensitive areas into a regional planning analysis of urban development. We established transition rules for a Markov-chain based CA model with multicriteria evaluation (MCE), then apply this model to the task of better retaining green infrastructure in an urban build-out plan. Two scenarios are examined and compared. The baseline scenario projects urban growth without environmental constraints while the GI conservation scenario incorporates an open-space conservation network intended to incorporate environmentally significant areas. The GI conservation scenario includes establishment of variable buffer widths on impaired streams as identified on the USEPA 303d list for stream impairment, protection of floodplains, and preservation of wetlands, steep slopes, and urban open space. We then quantitatively compare landscape configuration under the two different scenarios with structural metrics at the landscape level and discuss how this approach may guide not only retention of preserved land, but also increase the provision and quality of the ecosystem services offered by this land mass. In this research, we made an effort to simulate the changes of several land cover classes simultaneously. 2. Materials and methods 2.1. Study area The study area included twenty-seven counties in Ohio, Indiana, and Kentucky (Fig. 1). Fifteen of these are part of the CincinnatiMiddletown metropolitan statistical area. Twelve adjacent counties

lie on the fringe of the MSA and were incorporated in order to account for the edge effect (Zheng and Chen, 2000), which defines dynamics at the boundary where spatial characteristics of the landscape structure change. The geography of urban development in the Cincinnati-OH-KYIN Metropolitan Statistical Area (MSA) at the end of the last century was markedly characterized by two prevailing features: rapid expansion of low-density greenfield development and decline of the urban core and older suburbs. As Fig. 2 indicates, between 1990 and 2000 the population and employment densities in the majority of census tracts in and around the City of Cincinnati have declined. Some urban core tracts have lost as many as 14,465 persons/sq km while the urban fringe continued to expand outwards. A recent study on measuring sprawl found that the Cincinnati-OH-KY-IN MSA ranked 23rd most sprawling among 83 metropolitan regions in the United States in terms of residential density, 25th in terms of accessibility to road networks, and 32nd most sprawling in terms of job, housing and services mix (Ewing et al., 2003). The northern part of the region is characterized by agricultural land of high productivity, and this is also where development has been concentrated over the past decades as large proportion of the productive land had been converted to urban uses. Not surprisingly, most of the streams in this area tend to appear on the USEPA 303d listing of impaired streams. Given that the southern and western parts of the study region are particularly vulnerable to development – due to erosion risk, poor drainage, and seasonally shallow water tables – we focused on these areas as they might potentially profit from planning efforts informed by the results of our methods. 2.2. Data and data processing Land-cover land-use (LULC) data were obtained from the MultiResolution Land Characteristics Consortium (Vogelmann et al., 2001; Homer et al., 2007) which provided 1992 and 2001 National Land Cover Dataset (NLCD) in ArcGIS grid format. The classifications schemes of the two datasets were regrouped to reduce the number of classes and prepare the data for the simulation. Five classes were derived from the initial 1992 and 2001 NLCD datasets: water, woodland/open space, cropland/agricultural, wetlands, and urban/built-up land. For brevity, these five categories are further referred to as water, woodland, cropland, wetland and urban land. The water layer was based on the open water features included in the 1992 and 2001 NLCD datasets and the perennial streams network provided by the USGS National Hydrography Dataset (NHD) (USGS, 2004). The surface waters on the USEPA 303d list of impaired streams obtained through the BASINS data download tool (USEPA, 2004) were included as part of the GI conservation scenario. Woodland included deciduous, evergreen and mixed forests as identified by the 1992 and 2001 NLCD datasets. For the purposes of simulating the open space GI conservation scenario, this layer was combined with the open space data layer, provided by the Cincinnati Area Geographic Information Systems (CAGIS) as the latter included conservation easements, state parks, wildlife areas, forests, scenic parks, camp areas, bike trails, nature centers, conservancy districts, preserves, and other designated urban green space. Cropland and wetlands were derived from the 1992 and 2001 NLCD datasets to include cultivated crops and pastures, and woody and herbaceous wetlands, respectively. Finally, urban land was derived from commercial/industrial/transportation and low and high density residential categories as identified by the 1992 NLCD dataset, and low, medium and high intensity developed land as identified by the 2001 NLCD dataset. The road network for the study area was obtained from the USGS Seamless Server—National Atlas Roads database (USGS, 1999), which included major roads and ferry crossings (USGS, 1999). Population and employment data for the year 2000 and boundary files

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

143

Fig. 1. Location of the study area.

by census tracts were obtained for the U.S. Census Bureau (2000). Demographic and employment data for the year 1990 were downloaded from the National Historical Geographic Information System (NHGIS) database maintained by the Minnesota Population Center at the University of Minnesota (MPC, 2004). The NED elevation layer at 10 m resolution obtained from USGS was used to derive a topographic map for the study area, and areas susceptible to erosion hazard were determined based on slope gradient. The SSURGO database provided land productivity and habitat potential data in adjectival categories (e.g., good, fair, poor and very poor), and helped identify areas with seasonally high water table and water table depth annual minimum as well as areas with shallow depth to bedrock. SSURGO-derived data layers were derived to delineate floodplains based on flooding frequency observed as both dominant and maximum flood stage conditions. 2.3. Modeling approach Cellular automata (CA) models are based on the interaction of several components: the grid space, also known as “the lattice”, is 1, 2, or multi-dimensional space; the “cell” or the “automaton” is a discrete variable that represents the structural units of the lattice. A “cell state” is a description of the cell characteristics subject to change. The change occurs according to specific transition rules. Transition rules are mathematical expressions that govern changes

in the cell state. A self-reproductive cell is a point on the lattice that is able to assume a finite number of different states based on the unambiguous transition rule which may be interactive with the states of the immediate neighbors on the same lattice. The urban growth simulation model presented in this research incorporates cellular automata, Markov chain analysis, and multi-criteria evaluation techniques. We used a cellular automata—Markov chains (CA MARKOV) model in IDRISI, Andes, v15.0 (Eastman, 2006) to simulate the urban growth patterns in the study region. The module allows the simultaneous simulation of a group of land use/ land cover categories. The CA MARKOV module requires a land cover dataset to represent the initial states, a Markov transition matrix, a group of suitability images (one for each land cover class), a number of iterations and a contiguity filter. The transition rules are set up using a multi-criteria evaluation (MCE) and fuzzy membership functions to develop suitability maps for each simulated land cover class (Eastman, 2006). Suitability analysis ranks available land in a systematic procedure according to which the combined effects of various factors assumed to determine locational preferences are derived through evaluation, weighting, and overlay (Wang and vom Hofe, 2007). Five land cover classes were simulated simultaneously. During the iteration process, each land cover class became consecutively a host category competing against the “claimant” categories (Eastman, 2006). Fig. 3 displays a flow-chart of the simulation process:

144

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

Fig. 2. Change in population and employment densities per census tract (persons/sq km) in 1990–2000 Data source: U.S. Census Bureau.

Step 1: A Markov transition matrix establishes the frequency with which cells within a specific land cover class will transition to any other land cover class during the simulation period. The initial Markov transition probabilities matrix is based on the observed frequencies of transition between the initial and the latter land cover datasets using the MARKOV module in IDRISI (Eastman, 2006). A specific Markov chain matrix was derived for each ten-year period of land cover projections. The end year of projection is 2030. Step 2: Since the Markov process does not establish where land use transitions would occur, we used multi-criteria evaluation (MCE) involving fuzzy membership functions to determine

Fig. 3. A flow-chart of the modeling process.

the suitability and location of transitioning cells (Wu, 1998; Wu and Webster, 1998; Yeh and Li, 2001; Liu and Phinn, 2003). MCE employs constraints and factors. Constraints, expressed as Boolean images, in this context are used to constrain or limit the extent of development. Factors used in MCE account for suitability, accessibility, and neighborhood effects (Jiao and Boerboom, 2006). Local field characteristics such as soil properties, slope and elevation are used to determine and assign physical suitability scores. Accessibility is measured in terms of proximity to urban growth centers and major roads. We have taken into consideration residential and employment mobility between 1990 and 2000 by assigning higher weights to the areas that have experienced higher growth rates in terms of increasing residential and employment density. The neighborhood effect is accounted for as proximity to environmentally sensitive areas. We used monotonically increasing and decreasing S-shaped, J-shaped and linear fuzzy-logic functions (available within IDRISI) to develop suitability layers based on proximity for each land cover class included in the simulation. For example, proximity to water bodies for the purposes of suitability analysis was represented as a monotonically decreasing S-shaped function to account for a non-linear decrease in desirability with distance. Step 3: Assigning the number of iterations. The time span between the two land cover images for which Markov transition probabilities were estimated determined how many iterations would be used in the simulation. A contiguity filter assigned lower weights to suitability scores of isolated cells that did not belong to any of the nearby clusters of similar land cover cells (Eastman, 2006). The down-weighing did not surpass 90 percent which provided each floating pixel a certain chance to be assigned a category

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

145

Table 1 Structure and composition landscape metrics used in the analysis. Patch Analyst 4.0 (Rempel, 2008) Structure and composition metrics

Unit

Description

Number of patches (NP) Mean patch size (MPS)

n/a ha

Extent Patch size coefficient of variation (PSCoV) Patch size standard deviation (PSSD) Shannon’s Diversity Index (SHDI)

Percent Dimensionless ha Dimensionless

Shannon’s Evenness Index (SHEI) Edge density (ED)

Dimensionless m/ha

Total patch edge (TE)

in

Mean patch edge (MPE)

m/patch

Area-Weighted Mean Shape Index (AWMSI)

Dimensionless

Patch abundance and richness. The ratio of the total area to the number of patches. As it decreases, the fragmentation of the landscape increases. The proportional abundance of each class. Variability relative to the typical patch size. A measure of dispersion of the patch size distribution. SHDI describes a relationship between the number of classes, the total number of patches, and the relative abundance of patches in each class. It has a value of 0 when no diversity is present and increases as the landscape becomes more fragmented. SHEI reaches its maximum value of 1 when all classes have an equal number of patches. ED seeks to describe the total edge per unit area. It increases as the total patch perimeter increases. TE is defined as total perimeter of all patches. It increases as the number of patches and irregularity of the patch shape increases. The index is computed as the ratio of TE to NP. It increases as the landscape fragmentation increases. The index describes the complexity of the patch shape. It equals 1 when all patches have a circular shape and increases as the irregularity of the shape increases.

even if it was in close proximity to an existing cluster of cells of the same land cover class. Step 4: Scenario Development. Two scenarios were developed to compare the impact of our land use modeling approach. Both scenarios used fuzzy membership functions to evaluate a number of commonly applied factors such as distance to roads, employment growth centers, residential growth locations, streams and wetlands. The baseline scenario (the baseline scenario) involved no environmental constraints. By contrast, the green infrastructure conservation scenario incorporated an extended network of environmentally sensitive areas and a calculation of variable buffer widths for impaired streams included in the USEPA 303d list. Under this scenario, areas excluded from development included those with exceedingly shallow depth to seasonally high water table (less than 60 cm), shallow depth to bedrock (less than 80 cm), poorly drained soils (hydrologic soil group C and D), steep slopes (above 15 percent) and floodplains. In addition, factors based on open space features and variable buffers were also included in the GI conservation scenario. A method developed by Phillips (1989) and Xiang and Stratton (1996) evaluates the effectiveness of riparian buffers for non-point source pollution control based on soil characteristics, soil moisture, land cover and slope gradient. The method employs Phillips’ (1989) Riparian Buffer Delineation Equation (RBDE) which computes the buffer effectiveness ratio (dimensionless) in relation to the recommended buffer width (L−1 ), Manning roughness coefficient (dimensionless), slope (in percent), saturated hydraulic conductivity (LT−1 ), and soil storage moisture capacity (L). The data for the RBDE calculation were derived from the SSUSRGO database and the land cover dataset. We used this method to calculate variable buffer widths for the impaired streams identified under USEPA 303d list for stream impairment within the study area. Step 5: Model validation. We used the VALIDATE module in IDRISI to assess the level of agreement between the observed 2001 and simulated 2001 land cover datasets based on the Kappa spatial correlation statistic, which ranges between zero (random) and unity (perfect agreement) and is based on the pixel-by-pixel comparison of the observed data and the modeled result. 2.4. Comparison of the two scenarios using spatial statistics and landscape metrics Landscape analysis was performed using Patch Analyst 4.0 (Rempel, 2008), which integrates FRAGSTATS 2.0 (McGarigal and Marks, 1995) as an extension of ESRI’s ArcGISTM , and Conefor

Sensinode 2.2 (Saura and Torne, 2009). Tables 1 and 2 provide a description of the landscape metrics used in the analysis. The most common measures of spatial composition include the number of patches (NP), mean patch size (MPS), and extent which reflects the proportional abundance of each. The patch size coefficient of variation (PSCoV) indicates the degree of variability (low or high) relative to the typical patch size and is expressed as a percentage of MPS. The patch size standard deviation (PSSD) is a measure of dispersion of the patch size distribution. Shannon’s Diversity Index (SHDI) seeks to describe a relationship between the number of classes and the area covered by the patches within each class (Rempel, 2008). The index value is maximized when there is an even, proportional distribution of the area among all classes. Subsequently, the SHDI increases as the number of patches and the proportional distribution of the area covered by an individual patch type increases. Shannon’s Evenness Index (SHEI) emphasizes the equitable distribution of the areas covered by each class of patches (Rempel, 2008). A common measure of patch structure and distribution is the size and density of patch edges. Edge density (ED) is a measure of the total edge (as a perimeter measure) per unit area (Rempel, 2008). An increase in total patch edge (TE), mean patch edge (MPE), or edge density indicates proliferating fragmentation of the landscape. Patch shape complexity is measured by the area-weighted mean shape index (AWMSI) which is calculated as the summation of the individual patches perimeter-area ratios divided by the square root of patch area adjusted by the number of patches (Rempel, 2008). Table 1 describes the Patch Analyst 4.0 landscape metrics used in the analysis. Land Use Mix Index is a measure of dispersion developed by Frank and Pivo (1994) to assess the distribution of land use/land cover classes across the landscape. The index is calculated as follows:

 n

Land Use Mix = −

1

(pi )ln(pi )



ln(n)

where pi is the percentage of land use of type i in the area, and n is the number of land cover classes. Values close to 0 indicate a dominance of a particular land cover class, while values close to 1 reflect a blend of various similarly sized land cover classes across the landscape. Connectivity metrics were derived using Conefor Sensinode 2.2 (CS22) (Saura and Pascual-Hortal, 2007; Saura and Torne, 2009). CS22 calculates the number of links (NL) and number of components (NC), the Harary Index (H), the class coincidence probability

146

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

Table 2 Landscape connectivity metrics used in the analysis. Conefor Sensinode 2.2 (Saura and Torne, 2009) Connectivity metrics

Unit

Description

Number of links (NL)

n/a

Number of components (NC)

n/a

Harary Index (H)

Dimensionless

Class Coincidence Probability (CCP)

Dimensionless 0–1

Integral Index of Connectivity (IIC)

Dimensionless 0–1

Probability of Connectivity (PC)

Dimensionless 0–1

A link is the path between a pair of landscape patches. NL increases as the landscape connectivity increases. A component consists of interlinked patches. NC decreases as the landscape becomes more interconnected. H represents the number of minimum distances between patches that belong to the same component. Higher values of H indicate increased landscape connectivity. CCP is defined as the likelihood that two randomly selected patches belong to the same component. It yields higher values as connectivity improves. IIC associates patch area with the shortest paths between each pair of landscape patches. It represents both ease of access and habitat availability and increases as connectivity strengthens. PC is based on a stochastic landscape graph defined as the maximum probability paths between each pair of landscape patches. PC yields higher values as connectivity increases.

(CCP), the integral index of connectivity (IIC) and the probability of connectivity (PC). NL indicates the number of inter-patch connections within a specified threshold dispersal distance. In order to select the most appropriate threshold dispersal distance, we carried out a sensitivity analysis at 5000 m, 1000 m and 500 m. With the increase of inter-patch dispersal distance the open space distribution tend to coalesce into a single graph extending over the entire landscape, thus precluding valid conclusions about connectivity at more localized scales (Keitt et al., 1997). As a result of sensitivity analysis, we selected a threshold distance of 500 m, also used in previous studies (Cogan, 2001). NC provides the number of interconnected areas, i.e., the areas in which paths connect every pair of patches. H is based on the number of patches and the number of shortest topological distances between them (Saura and Torne,

2009). CCP indicates the probability that two randomly selected patches would fall within the same component. IIC is an improved habitat assessment measure as it links the area of each individual patch to the number of minimum distances between patches. PC applies the same logic but instead of shortest paths it uses the maximum product probability of all possible paths (Saura and Torne, 2009). All indices increase as the connectivity improves with the exception of NC which shows fewer components as the landscape becomes more interconnected. Conefor Sensinode 2.2 connectivity metrics used in the analysis are described in Table 2. We extracted CS22 indices for wetlands and five comparison areas to examine changes in connectivity under the two scenarios. The comparison areas were selected using the following criteria: areas with high potential for measuring change, areas where envi-

Fig. 4. Impaired streams on the 303d list and calculated variable buffer widths.

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

ronmentally sensitive areas are present, areas that would allow to predict conditions over a larger region, and areas that can show changes in response to management actions (i.e., variable buffers). Based on these criteria, we selected four areas located at the fringe of urbanization where the fastest land conversion occurs (near Taylor, Four Mile, Elk and Shaker creeks) and one in the urban core (Mill Creek) (Figs. 4 and 8). Four of them, with the exception of Shaker Creek, drain into streams that are on the 303d list of stream impairment. All five areas have known environmentally sensitive attributes and provide examples of changes that would occur over larger areas if no environmental constraints are considered. 3. Results and discussion 3.1. Land cover simulation results under two scenarios In this study, we used cellular automata—Markov chains model to simulate change in five land cover classes as a result of urban development (Fig. 5). Table 3 provides a summary of the Markov

147

probability matrix used for simulating the transitions between the initial cell state and the remaining four cell states represented as land cover categories. Table 3 indicates that during 1992–2001 there was 12.5 percent chance that cropland pixels would transition to urban land and 17.2 percent chance that they would transition to woodland. During the same period, there was 8.2 percent chance that woodland pixels would transition to cropland and 15.3 percent chance that they would transition to urban land. The contiguity filter which assigns lower weights to pixels that have high suitability scores but are not in close proximity to the land cover class which they belong to, helps reduce the number of unassigned pixels. As a result, in subsequent simulations the probability of transition from cropland to urban stabilizes at approximately 0.09, and from woodland to urban at approximately 0.04 (Table 3), thus reducing the variance in these elements of the transition matrix. The transition rules allowed the number and location of water pixels to remain unchanged, the wetland pixels to be preserved under the GI conservation scenario, the number of woodland and cropland pixels to decrease over time

Fig. 5. 2030 land cover projections for the baseline (a) and GI conservation (b) scenarios.

148

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153 Table 4 Summary statistics of projected land cover change in the study area by the year 2030.

Table 3 Markov transition probability matrix. Given

Probability of changing to Water 1 Woodland 2 Cropland 3 Wetland 4 Urban 5

1992–2001 Class 1 Water Class 2 Woodland Class 3 Cropland Class 4 Wetland Class 5 Urban 2001–2010 Class 1 Water Class 2 Woodland Class 3 Cropland Class 4 Wetland Class 5 Urban 2010–2020 Class 1 Water Class 2 Woodland Class 3 Cropland Class 4 Wetland Class 5 Urban 2020–2030 Class 1 Water Class 2 Woodland Class 3 Cropland Class 4 Wetland Class 5 Urban

0.7871 0.0046 0.0033 0.053 0.0045

0.1166 0.7586 0.1722 0.6393 0.1182

0.047 0.1515 0.6987 0.1758 0.0858

0.0106 0.0027 0.0007 0.021 0.0018

0.0387 0.0826 0.1251 0.111 0.7897

0.9979 0.0004 0.0013 0.0034 0

0.0011 0.8751 0.1134 0.3038 0.0019

0.0009 0.0805 0.7909 0.0514 0.0013

0.0001 0.0013 0.0004 0.6344 0.0001

0 0.0425 0.094 0.007 0.9966

0.9989 0.0005 0.0013 0.0033 0

0.001 0.9452 0.0428 0.0429 0

0.0001 0.0128 0.8624 0.0056 0

0 0.0001 0.0002 0.9411 0

0 0.0414 0.0934 0.0071 1

0.9998 0.0006 0.0012 0.0032 0

0.0001 0.9472 0.0417 0.0385 0

0.0001 0.012 0.8712 0.0044 0

0 0 0.0001 0.9474 0

0 0.0402 0.0858 0.0065 1

while the number of urban pixels to increase as urban development expands. The MARKOV module in IDRISI calculates Markov transition probabilities and outputs a text file with a transition probability matrix, a text file with the number of transitioning cells from one land cover class to another, and raster grids indicating Markov transition areas. Markov transition probability matrix is computed from cross-tabulation of earlier and subsequent land cover images. Markov transition areas are derived by multiplying each column representing a land cover category in Markov probability matrix (Table 3) by the number of cells of the same land cover class in the subsequent image (Eastman, 2006). As a CA-Markov model is set in motion, a 5 × 5 contiguity CA filter re-weights the suitability maps during each iteration increasing the suitability of pixels in close proximity to contiguous areas of the same category. The reweighted maps undergo a multi-objective land allocation process that resolves land allocation conflicts using the highest suitability score (Eastman et al., 1998). Each category consecutively becomes a “host” category while the remaining four complete through the multi-objective land allocation procedure (MOLA) (Eastman et al., 1998; Eastman, 2006). During each iteration, pixels with the highest transition probability and highest suitability score for a particular class transition to a new class while pixels with lower probabilities and lower suitability scores remain unchanged. If the input consists of 10 iterations, the model allocates 1/10 of all cells expected to transition to another land cover class during each iteration (Eastman, 2006). Fig. 6 illustrates Markov transition areas and suitability scores for four land cover classes near Four Mile Creek, Ohio, and the contiguity filter used in the study. Fig. 7 illustrates Markov transition areas from cropland to urban and woodland to urban per time span for the same area. Our projections indicate that approximately 30 percent of the study area will be urbanized by the year 2030 (Table 4), though the extent of urban development under the GI conservation scenario is slightly less than that under the baseline scenario While the total amount of urban development in the GI conservation scenario is similar to that of the baseline scenario, the proportion of urban patches in environmentally sensitive areas is reduced from

Variables

Baseline scenario

GI conservation scenario

Total urban area (ha) Total sensitive areas (ha) Percent urban land Percent sensitive areas Urban development in the sensitive areas (ha) Percent urban development in sensitive areas

438,874 226,723 0.33 0.16 108,989

421,319 226,723 0.29 0.16 48,715

0.48

0.22

48 to 22 percent. Herold et al. (2005) reported similar results when comparing five alternative urban growth scenarios including MSQ (maintaining the status quo) and EEP (maximum protection of environmentally sensitive areas). The study found general loss of open space and natural corridors under the MSQ scenario, and increase in the number of urban patches under the EEP scenario due to the spatial regulation of environmentally sensitive areas (Herold et al., 2005). The validity of the model results have been evaluated by comparing the projected land cover image for 2001 with the existing 2001 land cover map. The overall Kappa statistics for the five projected land cover classes was 0.7 which showed a very good agreement between projected and observed land cover layers. A separate Kappa statistics was calculated for the agreement between observed and projected built-up areas. The overall Kappa statistics for the urban class was 0.75 which indicated a very good agreement between observed and projected urbanized areas. The results for urban land show that the model predicts accurately the location of 88 percent of the urban pixels. The possible sources of the unexplained variation are potentially due to the inherent uncertainty of the land cover data as well as the uncertainty associated with the modeling process. The accuracy of NLCD maps, processing errors as pixel assignment, and the edge effect may have contributed to variance in these elements of the transition matrix. 3.2. Comparing the two scenarios using spatial statistics and landscape metrics To further describe differences between the two scenarios, we used a set of landscape-level metrics to quantify the distribution of land uses and their spatial relationships. The number of patches and the patch density in the projected land cover under the GI conservation scenario increased, compared to the baseline scenario (Table 6), indicating a slightly more fragmented landscape (Furberg and Ban, 2008), which we attribute to an increase in the number of green infrastructure patches under the set of constraints prescribed to foster conservation and linkages among natural areas. For example, the number of wetland patches increased from 2769 to 6763 as a result of their inclusion in the green infrastructure network, woodland patches increased from 18 to 24 percent under the GI conservation scenario, and the percent of cropland patches from 22 to 41 percent. Importantly, an increase in all of these more natural land uses apparently led to the urban patches becoming interspersed with these “greener” land cover categories. However, the cropland category is somewhat ambiguous, as it implies the implementation of any of several agricultural management practices including row-crops, horticulture, silviculture, among other more specialized practices. Given the size of the relatively small size and continuity of these patches, however, it is possible that management of this land would be limited to practices that are economical at smaller scales, which may include small-scale diversified farms, conservation set-aside, greenways. These practices

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

149

Fig. 6. An illustrative example of a contiguity filter, Markov transition probabilities and suitability scores resulting from the MCE procedure near Four Mile Creek, Ohio.

would be expected to have less impact on the environment than that of large-scale, mechanized row-crop agriculture. The value of the Shannon’s Diversity Index (SHDI) for the baseline scenario and the GI conservation scenario does not change because under both scenarios the number of classes remains the same, and there are no large differences in the proportional distribution of the area covered by each land cover class. However, the increase in interspersion among land classes is quantified by the values of the Shannon’s Evenness Index (SHEI), whereas a slight increase in SHEI under the GI conservation scenario suggests that the landscape composition has become more even indicating a

transition from a dominant land cover class towards a more diversified landscape that is heterogeneously distributed. The increased interspersion of land uses is also reflected in the area-weighted mean shape index (AWMSI) which increased from 13.34 under the baseline scenario to 25.49 under the GI conservation scenario. The increased value of AWMSI can be attributed to increased morphological complexity of the landscape. As morphological complexity increases, the smoothness, regularity and compactness of the landscape decrease. As a result, the more subdivided and interspersed the landscape is, the more morphologically complex it becomes. Consequently, the AWMSI value for the baseline scenario

150

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

Fig. 7. An illustrative example of transitioning areas per time span near Four Mile Creek, Ohio: cropland to urban and woodland to urban.

is considerably lower than that for the GI conservation scenario (Table 5). The edge density (ED) (m/ha) which represents the patch perimeter length per unit area increased substantially from 61.21 m/ha under the baseline scenario to 92.79 m/ha under the GI

conservation scenario (Table 5). This result suggests that the relatively homogeneous urban landscape simulated under the baseline scenario has become more complex and heterogeneous as it is interspersed with elements of the open space conservation network. As Table 5 indicates, the mean patch edge (m/patch) does

Table 5 Landscape metrics for various land cover classes. Landscape metrics

SHDI

Baseline scenario Landscape 1.23 Water Woodland Cropland Wetlands Urban G1 conservation scenario Landscape 1.22 Water Woodland Cropland Wetlands Urban

SHEI

AWMSI

ED (m/ha)

MPE (m/patch)

MPS (ha)

PSCoV

PSSD (ha)

% Area

% NumP

0.69

13.34 9.91 13.38 5.30 1.67 21.45

2.32 18.20 19.69 0.67 20.15

716.79 1792.48 1594.60 349.71 689.29

5.78 36.93 23.13 0.93 10.45

3518.33 2942.97 1198.74 324.35 11314.91

203.27 1086.84 277.32 3.01 1182.87

0.02 0.36 0.29 0.002 0.33

0.06 0.18 0.22 0.03 0.51

25.49 11.68 20.25 3.78 1.56 50.27

2.89 31.82 28.87 1.19 27.56

799.31 1721.24 902.93 252.19 1444.26

5.10 31.07 17.88 0.39 17.04

3327.33 3348.78 1167.46 349.36 9775.14

169.63 705.73 91.96 1.37 1665.34

0.02 0.41 0.28 0.002 0.29

0.05 0.24 0.41 0.06 0.25

0.78

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

151

Table 6 Reported metrics for open space connectivity at five comparison areas. Connectivity metrics Baseline scenario NL NC H CCP IIC PC GI conservation scenario NL NC H CCP IIC PC

Wetlands

Great Miami River near Taylor Creek

Mill Creek

Four Mile Creek

Great Miami River b/n Four Mile and Elk Creek

Little Miami River near Shaker Creek

565 47 841 0.198 0.146 0.158

19 15 24 0.071 0.122 0.102

1467 31 3360 0.184 0.290 0.193

120 10 392 0.570 0.651 0.4348

131 13 243 0.524 0.380 0.567

42 13 51 0.330 0.472 0.357

2778 65 4468 0.342 0.262 0.284

210 1 1085 1.000 0.922 0.877

5194 20 18325 0.652 0.585 0.714

1340 4 10072 0.988 0.675 0.738

1611 17 6978 0.935 0.972 0.901

739 6 3274 0.475 0.523 0.685

not change for woodland but drops significantly for cropland which confirms the observation that development in the study area occurs primarily on prime agricultural land. The mean patch edge of urban areas increases twofold and since there is no change in the amount of the total urbanized area, we assume that the increase in mean patch edge is due to greater variability in the arrangement of urban patches under the GI conservation scenario. The mean patch size (MPS) slightly drops under the GI scenario for woodland, wetlands and agricultural land, but increases for urban land which suggests the presence of fewer but larger urban patches. This observation is confirmed by the number of patches which for the baseline scenario (without environmental constraints) is substantially higher (41,981) compared to that for the GI conservation scenario (27,281).

Although the amount of urbanized land under the two scenarios is only slightly different, the proportional abundance of each class varies considerably. The proportion of urban patches drops from 0.51 (the baseline scenario) to 0.25 (the GI conservation scenario) as a result of greater fragmentation. Under the GI conservation scenario, the patch size standard deviation decreases for all land cover types except build-out areas. The coefficient of variation shows slight variability for woodland, cropland and wetland but declines considerably for urban land. The discrepancy in these two results can be explained by the way the two metrics are derived. As PSCoV indicates variation as proportion of the mean, it is preferred when the standard deviation is influenced by the number of patches. Therefore, we can assert that the variability in patch

Fig. 8. Comparison of three selected areas with and without environmental constraints.

152

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

size distribution of urban patches decreases for the GI conservation scenario. The patterns revealed by these metrics suggest that under the baseline scenario there is a consolidated urbanized area in which the existing intra-urban open space is gradually buildout, thus creating a smaller number of larger urban patches and significant variability in patch size. The number of urban patches increases for the GI conservation scenario as the build-out area is subdivided by open space, natural corridors and variable riparian buffer zones. Both scenarios reveal that development occurs primarily outward of the existing urban core following an approximately concentric pattern. Smaller or satellite urban cores also formed around existing urban centers as the expansion of the urban land is projected into the future. Urban cores under both scenarios are characterized by larger urban patches at the center and a mosaic of more fragmented land cover classes at the rural/urban periphery. Fig. 5 indicates that the fragmentation of the urban core under the GI conservation scenario reflects the presence of more natural areas and connected vegetated corridors along the major watercourses in the region. The land cover mix index increased from 0.68 under the baseline scenario to 0.73 under the GI scenario, indicating that under the GI scenario, the built environment has retained more of its green features thus creating a more diverse and balanced mix of landscape features. Table 6 provides a summary of the connectivity metrics for wetlands and open space under the two scenarios. The results indicate that the number of links in wetland areas grow from 565 to 2778 and the Harary Index increases from 841 to 4468. The class coincidence probability, the integral index of connectivity, and the probability of connectivity also improve, reflecting an improvement in interconnectedness under the GI conservation scenario. The number of links between the green patches in all five comparison areas increases in order of magnitude. The number of components drops for four of the five areas, thus indicating that the paths between patches have increased creating larger interconnected regions under the GI conservation scenario. Under the GI conservation scenario, the Harary Index increases exponentially for all selected comparison areas. The most significant change occurs near Four Mile Creek, between Four Mile and Elk Creek, and near Shaker Creek (Fig. 8). These areas are located in the northern part of the study region where the rate of land conversion to urban uses is among the highest. CCP indicates a substantial increase in the same portion of the study region as well as at Taylor Creek. IIC and PC have higher value for all five comparison areas, thus reflecting a stable trend towards improved connectivity in the study area after the inclusion of a green infrastructure network in the simulation model. The importance of the variable buffer widths for the improved connectivity is difficult to assess directly. 4. Conclusion A focal objective of urban growth modeling is evaluating possible future paths of development (Herold et al., 2005). In this study, we have examined the spatial structure of continued urban growth under two scenarios: with and without green infrastructure conservation. The approach based on cellular automata, Markov probabilities and multi-criteria evaluation allowed the simulation of five land cover classes simultaneously and thus provided a basis for meaningful examination of the landscape dynamics under different scenarios. Setting transition rules in the form of suitability maps based on MCE allows for consideration of various factors that are commonly used in land use planning and decision-making such as changes in population and employment density, proximity to roads, proximity to waterbodies as well as protection of open space, water resources, riparian corridors, and other environmentally sensitive areas. This

type of systematic evaluation in a GIS environment combined with foresight scenarios and projections equips environmental and land use planners with a useful tool that can integrate activities consistent with the smart growth implementation options and strategies. Two limitations of the proposed approach require further investigation. First, the model is scale-dependent. Its results hinge on the scale at which grids are aggregated assuming that all the processes in the landscape are occurring at that scale. However, both urban development and ecological effects are multi-scalar processes that involve dynamic responses at various levels. In order to account for these responses, in this research, we integrated three levels of spatial resolution: a pixel, comparison area, and region. Other approaches are worth further investigating. Alberti and Waddell (2000), for example, suggested treating land use at a parcel level, and land cover—at a patch level in order to respond to the underlying dynamics. Second, preserving open space and critical environmental areas depends largely on local government regulation and local-level land use planning. Each regulation must incorporate acceptable methodology for critical area identification based on current scientific research. Therefore, the results of this research are descriptive rather than prescriptive. The results show that under the GI conservation scenario the overall extent of urban development decreased, the mean patch size for urban land increased, and the percent of urban patches decreased as they become interspersed with a variety of other patches. The strength of the GI conservation scenario is also related to its ability to incorporate variable buffer widths for the impaired streams on the USEPA 303d list. As the impaired streams oftentimes cross jurisdictional boundaries, their incorporation in an urban growth simulation suggests a more holistic approach towards the management of such impaired waters in the projected future. For municipalities that require comprehensive planning account for present and future environmental regulations, this type of modeling and prediction will play a significant role in illustrating the probable outcomes of different land use choices upon expanding the extent of green infrastructure and the numerous ecosystemlevel services that are provided by such area. Local governments have at their disposal several planning tools that can help construct, operate and maintain the green infrastructure network. Acquisition through easements, special area zoning, transferable development rights, specific water quality standards for designated areas, wetland restoration, wide-ranging watershed planning, and outreach programs can become part of directing regional (USEPA, 2005). In addition, Arendt (2004) recognized that preserving open space at a site scale does not prevent the cumulative development impacts on natural systems and suggested that a community-wide approach be used for preserving the connectivity of green networks, and that local land-use regulations that would identify and protect potential conservation areas would play a role in this process. Combined with the capabilities of urban growth models, including cellular automata these planning instruments can be used to generate and assess different scenarios of urban growth. A public discussion of these scenarios involving stakeholder participation can bring awareness of the expected change and valuable inputs to the planning process. Given transparency of the process and ability to visually and quantitatively compare development scenarios accounting for the non-linearities that are not considered in linear GIS map algebra, this may offer an improvement over the more qualitative and subjective suitability overlays. Urban growth projections under various conservation scenarios can be used to resolve conflicts of interests in land development and build consensus among stakeholders about the future patterns and directions of urban development, if visions differ. As part of future research, the model developed in this study can incorporate various growth management strategies, zoning ordinances or other regulatory tools

D. Mitsova et al. / Landscape and Urban Planning 99 (2011) 141–153

that may influence and redirect the patterns and pace of urban development. References Alberti, M., Waddell, P., 2000. An integrated urban development and ecological simulation model. Integr. Assess. 1, 215–227. Arendt, R., 2004. Linked landscapes: creating greenway corridors through conservation subdivision design strategies in the northeastern and central United States. Landscape Urban Plan. 68, 241–269. Batty, M., Xie, Y., 2005. Urban growth using cellular automata models. In: Mcguire, D., Batty, M., Goodchild, M.F. (Eds.), GIS, Spatial Analysis and Modeling. ESRI Press, Redlands, CA, pp. 151–172. Benedict, M., McMahon, E.T., 2006. Green Infrastructure: Linking Landscape and Communities. Island Press, Washington, DC. Claggett, P., Jantz, C.A., Goetz, S.J., Bisland, C., 2004. Assessing development pressure in the Chesapeake Bay watershed: an evaluation of two land-use change models. Environ. Monitor. Assess. 94, 129–146. Clarke, K.C., Hoppen, S., Gaydos, L., 1997. A self-modifying cellular automaton model of historical urbanization in the San Francisco Bay area. Environ. Plann. B 24, 247–261. Clarke, K.C., Gaydos, L., 1998. Loose-coupling a cellular automaton model and GIS: long-term urban growth prediction for San Francisco and WashingtonBaltimore. Int. J. Geogr. Inf. Sci. 12, 699–714. Cogan, C., 2001. Biodiversity predictions: integrating urban growth models with land cover data and species habitat information. USGS GAP Anal. Bull. 10 (4), 33–36. Costanza, R., d’Arge, R., de Groot, R., Farber, S., Grasso, M., Hannon, B., Limburg, K., Naeem, S., O’Neill, R.V., Paruelo, J., Raskin, R.G., Sutton, P., van den Belt, M., 1997. The value of the world’s ecosystem services and natural capital. Nature 387, 253–260. Dosskey, M., 2001. Toward quantifying water pollution abatement in response to installing buffers on crop land. J. Environ. Manage. 28 (5), 577–598. Duany, A., Speck, J., Lydon, M., 2009. The Smart Growth Manual. McGraw-Hill, New York, NY. Eastman, J.R., Jiang, H., Toledano, J., 1998. Multi-criteria and multi-objective decision making for land allocation using GIS. In: Beinat, E., Nijkamp, P. (Eds.), Multi-criteria Analysis for Land-Use Management. Kluwer Academic Publishers, Dordrecht, pp. 227–251. Eastman, J.R., 2006. IDRISI Andes. Clark University, Worcester, MA. Eppink, F., van den Bergh, J., Rietveld, P., 2004. Modeling biodiversity and land use: urban growth, agriculture and nature in a wetland area. Ecol. Econ. 51, 201–216. Ewing, R., Pendall, R., Chen, D., 2003. Measuring Sprawl and Its Impacts. Smart Growth America, Washington, DC, Available at www.smartgrowthamerica.org (accessed August 14, 2010). Frank, L.D., Pivo, G., 1994. Impacts of mixed use and density on utilization of three modes of travel: single-occupant vehicle, transit, and walking. Transport. Res. Rec. 1466, 44–52. Furberg, D., Ban, Y., 2008. Satellite monitoring of urban sprawl and assessing the impact of land cover changes in the Greater Toronto Area. Int. Arch. Photogram. Rem. S., Spatial Inf. Sci. 37 (Part B8), 131–137. Herold, M., Couclelis, H., Clarke, K.C., 2005. The role of spatial metrics in the analysis and modeling of urban land use change. Comput. Environ. Urban 29, 369–399. Homer, C., Dewitz, J., Fry, J., Coan, M., Hossain, N., Larson, C., Herold, N., McKerrow, A., Van Driel, J.N., Wickham, J.D., 2007. Completion of the 2001 National Land Cover Database for the conterminous United States. Photogramm. Eng. Rem. S. 73, 337–341. Jiao, J., Boerboom, L., 2006. Transition rules elicitation methods for urban cellular automata models. In: Van Leeuwen, J.P., Timmermans, H.J.P. (Eds.), Innovations in Design & Decision Support Systems in Architecture and Urban Planning. Springer, Dordrecht, pp. 53–68. Keitt, T.H., Urban, D.L., Milne, B.T., 1997. Detecting critical scales in fragmented landscapes. Conserv. Ecol. 1 (1), 4. Liu, Y., Phinn, S.R., 2003. Modeling urban development with cellular automata incorporating fuzzy-set approaches. Comput. Environ. Urban 27 (6), 637–658. Longley, P., Batty, M., 1996. Analysis, modeling, forecasting, and GIS technology. In: Longley, P., Batty, M. (Eds.), Spatial Analysis: Modeling in a GIS Environment. Geo-Information International, Cambridge, UK, pp. 1–15.

153

McGarigal, K., Marks, B.J., 1995. FRAGSTATS: spatial pattern analysis program for quantifying landscape structure. USDA For. Serv. Gen. Tech. Rep. PNW-351. Minnesota Population Center (MPC), 2004. National Historical Geographic Information System: Pre-release Version 0.1. University of Minnesota, Minneapolis, MN, http://www.nhgis.org. Osborne, L., Kovacic, D.A., 1993. Riparian vegetated buffer strips in water quality restoration and stream management. Freshwater Biol. 29, 243–258. Phillips, J.D., 1989. An evaluation of the factors determining the effectiveness of water quality buffer zones. J. Hydrol. 107, 133–145. Polyakov, V., Fares, A., Ryder, M.H., 2005. Precision riparian buffers for the control of nonpoint source pollutant loading into surface water: a review. Environ. Res. 13, 129–144. Randolph, J., 2004. Environmental Land Use Planning and Management. Island Press, Washington, DC. Rempel, R., 2008. Patch Analyst v4.0, Available at http://flash.lakeheadu.ca/∼rrempel/patch/. Saura, S., Pascual-Hortal, L., 2007. A new habitat availability to integrate connectivity in landscape conservation planning: comparison with existing indices and application to a case study. Landscape Urban Plan. 83 (2–3), 91–103. Saura, S., Torne, J., 2009. Conefor Sensinode 2.2: a software package for quantifying the importance of habitat patches for landscape connectivity. Environ. Model. Software 24, 135–139. Smart Growth Network, 1996. Principles of Smart Growth, Available at http://www.smartgrowth.org/about/principles/ (accessed August 26, 2010). Syphard, A., Clarke, K.C., Franklin, J., 2005. Using a cellular automaton model to forecast the effects of urban growth on habitat pattern in southern California. Ecol. Compl. 2, 185–203. Torrens, P., 2003. Automata-based models of urban systems. In: Longley, P.A., Batty, M. (Eds.), Advanced Spatial Analysis: the CASA Book of GIS. ESRI Press, Redlands, CA, pp. 61–81. Trust for Public Land (TPL), 2005. GreenPrint for King County, WA. WLRD Visual Communications and GIS Unit Archives, Seattle, WA, Available at http://your.kingcounty.gov/dnrp/library/2005/KCR1856/0505 Greenprint.pdf (accessed August 15, 2010). U.S. Census Bureau, 2000. American FactFinder for Selected Ohio Counties. U.S. Environmental Protection Agency (USEPA), 2004. Better Assessment Science Integrating Point and Nonpoint Sources (BASINS), Version 4.0. Office of Water, Washington, DC. U.S. Environmental Protection Agency (USEPA), 2005. Riparian Buffer Width, Vegetative Cover, and Nitrogen Removal Effectiveness: A Review of Current Science and Regulations. EPA/600/R-05/118. Office of Research and Development, National Risk Management Research Laboratory, Ada, OK. U.S. Geological Survey (USGS), 1999. Major Roads of the United States. U.S. Geological Survey, Reston, VA, Available at http://nationalatlas.gov/atlasftp.html (accessed June 2008). U.S. Geological Survey (USGS), 2004. National Hydrography Dataset (NHD). USGS, Lakewood, CO, Available at http://nhdgeo.usgs.gov/viewer.htm (accessed June 2008). Vogelmann, J.E., Howard, S.M., Yang, L., Larson, C.R., Wylie, B.K., Van Driel, N., 2001. Completion of the 1990s National Land Cover Data Set for the Conterminous United States from Landsat Thematic Mapper data and ancillary data sources. Photogramm. Eng. Rem. S. 67, 650–662. Wang, X., vom Hofe, R., 2007. Research Methods in Urban and Regional Planning. Springer and Tsinghua University Press, Beijing. Wu, F., 1998. Simulating urban encroachment on rural land with fuzzy-logiccontrolled cellular automata in a geographical information system. J. Environ. Manage. 53 (4), 293–308. Wu, F., Webster, C.J., 1998. Simulation of land development through the integration of cellular automata and multi-criteria evaluation. Environ. Plann. B 25, 103–126. Xiang, W., Stratton, W., 1996. The b-function and variable stream buffer mapping: a note on A GIS method for riparian water quality buffer generation. Int. J. Geogr. Inf. S. 10 (4), 499–510. Yeh, A.G, Li, X., 2001. A constrained CA model for simulation and planning of sustainable urban forms by using GIS. Environ. Plann. B 28, 733–753. Zheng, D., Chen, J., 2000. Edge effects in fragmented landscapes: a generic model for delineating area of edge influences (D-AEI). Ecol. Model. 132, 175–190.