Detailed Methods and Results

Manage For Ecological Integrity

Comer, P. J. & J. Hak.2012. Landscape Condition in the Conterminous United States. Spatial Model Summary. NatureServe, Boulder, CO.

Conceptual Basis: Integrated, quantitative expressions of anthropogenic stress over large geographic regions can be valuable tools in environmental management. When they take the form of a map, they characterize ecological conditions on the ground; from highly disturbed to apparently unaltered conditions. They can be particularly helpful for identifying relatively intact landscape blocks or for screening ecological reference sites;i.e., a set of sites where anthropogenic stressors range from low to high. Ecological condition of reference sites are often further characterized in the field to determine how ecological attributes are responding to apparent stressors. This knowledge may then apply in other similar sites for all forms of environmental decision making.

This Landscape Condition Model integrates readily available spatial data in order to express common ecological stressors. The intent ofthe model is to enable spatial expression of common knowledge and experience regarding the relative effects of land uses on natural ecosystems and habitats. Expert knowledge forms the basis of stressor selection, and relative weighting in the model. This model has been calibrated westwide, and continues to be evaluated with field samples.

Technical Description: Table 1 includes a brief summary of data sets, settings, and assumptions included in this model. We selected a limited set of stress-inducing land use classes for which we have regionally consistent coverage. Our aim here is to characterize the primary local scale stressors. We have not attempted to factor in regional stressors, such as air pollutants or climate change. Stressors are organized into thematic groupings of Transportation, Urban and Industrial Development,and Managed & Modified Land Cover. Transportation features, derived from ESRI StreetMap data circa 2006, depict roads of five distinct sizes and expected traffic volume. These data provide a practical measure of human population centers and primary transportation networks that link those centers. Ecological stress induced by built infrastructure (through habitat loss, fragmentation, altered ecological processes, etc.) are well known.

As a compliment to Transportation infrastructure,Urban and Industrial Development includes industrial (e.g., mines, energy development) and built infrastructure across a range of densities, from high density urban and industrial zones, to suburban residential development, to exurban residential and urban open spaces (golf courses, for outdoor recreation. These data were derived mostly from national land cover data through combined efforts of USGS (National Land Cover and Gap Analysis Programs) and the LANDFIRE effort. Other data sets included oil/gas well and transmission line right-of-ways.

The third category, Managed and Modified Land Cover, includes the gradient of land cover types that reflect land use stressors at varying intensities. Again, national data from USGS and LANDFIRE provide a consistent depiction of these varying land cover classes,from intensive (cultivated and/or irrigated) agriculture, vineyards and timber tree plantations, various forms of introduced non-native vegetation in upland and wetland environments, and finally, areas where native vegetation predominates, but modifications have clearly taken place. These modifications include recently logged areas, or areas that have seen historic conversion, but have recovered some combination of mainly native vegetation (old fields, etc.).

Model Development: Each data layer is given a site impact value, scaled from 0.0 to 1.0 reflecting expert assumptions of the generalized ecological impact where the land use occurs, with values closer to 1.0 expressing relatively little ecological impact from the land use. A second "distance decay"function calculates and applies a decreasing value for each input layer with distance away from its location. Therefore, a given land use,such as a road of a given size and presumed traffic volume will be given two values, one for its relative impact where it occurs, and a second for the rate of decay of its presumed impact with distance. The result for each input layer is a map surface indicating relative scores between 0.0 and 1.0. Distance decay settings may vary from 0.0 - effectively no presumed ecological impact within one pixel of its location - out to a maximum of 2,000 meters, where presumed effects of a given land use would finally reach zero.

Individual spatial models for each input layer are then combined and normalized to a 0.0 to 1.0 scale. Where the lowest individual layer score is lower than the resulting normalized score, that lower score overrides the normalized score.The combination of per-pixel scores results in a continuous map surface.

Table 1. Date inputs and values integrated together for the NatureServe Landscape Condition Model.

Theme	Impact Score	Presumed Relative Stress	Decay Score	Impact Decays to Zero
Transportation
Dirt roads, 4-wheel drive	0.7	Low	0.5	200m
Local, neighborhood and connecting roads	0.5	Medium	0.5	200m
Secondary and connecting roads	0.2	High	0.2	500m
Primary Highways with limited access	0.05	Very High	0.1	1000m
Primary Highways without limited access	0.05	Very High	0.05	2000m
Urban and Industrial Development
Low Density Development	0.6	Medium	0.5	200m
Medium Density Development	0.5	Medium	0.5	200m
Powerline/Transmission lines	0.5	Medium	0.9	100m
Oil /gas Wells	0.5	Medium	0.2	500m
High Density Development	0.05	Very High	0.05	2000m
Mines	0.05	Very High	0.2	500m
Managed and Modified Land Cover
Ruderal Forest & Upland	0.9	Very Low	1	0m
Native Veg. with introduced Species	0.9	Very Low	1	0m
Recently Logged	0.9	Very Low	0.5	200m
Managed Tree Plantations	0.8	Low	0.5	200m
Introduced Tree & Shrub	0.5	Medium	0.5	200m
Introduced Upland grass & forb	0.5	Medium	0.5	200m
Introduced Wetland	0.3	High	0.8	125m
Cultivated Agriculture	0.3	High	0.5	200m

Invasive Annual Grass Models

Wildfire Regime Condition Class

Wildfire is a key natural process for many terrestrial CEs within each ecoregion but land use patterns commonly result in significant departure from expected fire frequency and intensity. In a limited way, we will develop spatial models of wildfire risk based on lightning strike and landscape information, as was completed in the Northern Basin and Range ecoregion. However,most aspects of these CAs are best addressed within the context of major coarse-filter CEs since existing knowledge and modeling centers around their characteristic fire regimes. This knowledge forms the basis for conceptual tabular and spatial models of fire regime departure and enables us to summarize these effects by appropriate landscape units (e.g., watersheds by 5^thlevel hydrologic unit codes or HUC10). Fire regime models also provide one key mechanism for translating measured and predicted trends in climate regimes as they affect these critical ecological dynamics.

Manage for Climate Change Refugia

2060 Climate Refugia for Terrestrial Ecosystems

The 4 kilometer California Academy Terrestrial Ecosystem Climate Envelope change summary data, which predicts where vegetation ranges will contract, expand or remain the same, were used to create a 2060 terrestrial ecosystem climate refugia map. Fifteen terrestrial ecosystem species were evaluated (these are all the terrestrial ecosystems that were evaluated in the California Academy Climate Envelope project): Great Basin Pinyon Juniper, Great Basin Xeric Mixed Sagebrush Shrubland, Inter-Mountain Basins Aspen Mixed Conifer Forest and Woodland, Inter-Mountain Basins Big Sagebrush Shrubland, Inter-Mountain Basins Big Sagebrush Steppe, Inter-Mountain Basins Curl-leaf Mountain Mahogany Woodland and Shrubland, Inter-Mountain Basins Montane Sagebrush Steppe, Inter-Mountain Basins Mixed Salt Desert Scrub,Inter-Mountain Basins Semi-Desert Shrub Steppe, Inter-Mountain Basins Subalpine Limber Bristlecone Pine Woodland, Mojave Mid-Elevation Mixed Desert Scrub,Rocky Mountain Aspen Forest Woodland, Sonora Mojave Creosotebush White Bursage Desert Scrub, Sonora Mojave Semi Desert Chaparral, and Sonoran Mojave MixedSalt Desert Scrub.

Each terrestrial ecosystem climate envelope change summary map was reclassified to identify the areas that are predicted to remain the same and assigned a value of 1, while all other pixels (areas of expansion or contraction or no occurrence) were assigned a value of 0.

In the Raster Calculator, each of the fifteen reclassified terrestrial ecosystem grids were added together to produce a terrestrial ecosystem climate refugia grid. Pixel values ranged from 0 to 8. The numeric value of a pixel represents how many different terrestrial ecosystems overlap;the higher the number the more terrestrial ecosystems that are predicted to remain in that location.

2060 Climate Refugia for Landscape Species

The 4 kilometer California Academy Landscape Species Climate Envelope change summary data, which predicts where species range will contract,expand or remain the same, was used to create a 2060 landscape species climate refugia map. Eight landscape species are evaluated in the Yale study. bighorn sheep, desert tortoise Mojave, greater sage grouse, mule deer summer, mule deer winter, mule deer yearlong and pygmy rabbit.

Each landscape species climate envelope change summary map was reclassified to identify the areas that are predicted to remain the same and assigned a value of 1, while all other pixels (areas of expansion or contraction or no occurrence) were assigned a value of 0.

In the Raster Calculator, each of the eight reclassified landscape species grids were added together to produce a landscape species climate refugia grid. Pixel values ranged from 0 to 6. The numeric value of a pixel represents how many different landscape species overlap; the higher the number the more landscape species that are predicted to remain in that location.

GEOPHYSICAL/BIOPHYSICAL HETEROGENITY DENSITY MAPPING

Abstract Summary:

Natureserve modeled the relative complexity of the physical landscape across the Yale study area for use in identifying potential climate refugia. The assumption is that the more complex the underlying physical landscape, the more likely it will be able to sustain species in the future if when species ranges shift due to climate change (Anderson and Ferree 2010). A quot;Geophysical Heterogenity Index" (GHI) map was derived from a combination of landform, solar radiation and flow accumulation maps. The relative density of the resultant GHI classes was then summarized in hexagon maps of varying scales (e.g. 1/2 km, 1 km, 4 km, 8 km, 16km area hexagons) to identify the geophysical heterogeneity density (GHD) across the Yale landscape.  In addition, a Biophysical Heterogeneity Density (BHD) map was produced from LANDFIRE Biophysical Setting vegetation data, to evaluate if the resultant maps are comparable to the Geophysical Heterogeneity Density(GHD) map, to evaluate if this readily available source of data on vegetation could be effectively utilized to identify climate refugia.

Methods and Results:

The landforms and solar radiation maps were produced using models (with slight modifications) from Jeff Jenness's Topographic Toolbox 9.3 (i.e. Landform Classification (Jenness).tbx and Solar Radiation (McCune 2002).tbx). The Flow accumulation map was modeled using standard ArcGIS hydrology modeling tools. These three maps were then combined to create a Geophysical Heterogeneity Index (GHI) map. The Geophysical Heterogeneity Density (GHD) hexagon map was created from a simple count of the total number of GHI classes within each hexagon. Detailed descriptions of the technical methodologies for modeling the landforms, solar radiation and flow accumulation maps, and producing the geophysical heterogeneity index map, and the final geophysical and biophysical heterogeneity density maps are provided in the appendix.

The landforms map effectively models the physical complexity of the Central Basin and Range / Mojave Basin and Range landscape into ten landform classes. The nuances of the landscape are fairly well defined (see Map 1). However, the landform map tends to over-estimate the extent of the marco-scale landform classes (i.e. u-shaped valleys, plains, open slopes and upper slopes/mesas), as well as the extent of canyons/deeply incised streams, and mountain tops/high ridges, but tends to under-estimate the extent of micro-scale landforms (i.e. midslope drainages/shallow valleys, upland drainages/headwaters, local ridges/hills in valleys, and midslope ridges/small hills in plains). The model parameters were adjusted to provide the best compromise between reducing the macro-scale landforms and increasing the micro-scale landforms. On alluvial plains, upperslopes/mesas extend down too far into the open slopes, and should more correctly be identified as open slopes. This issue was unresolvable, and simply reflects that alluvial fan elevations are similar in structure to upperslopes/mesas, and therefore will tend to be incorrectly classified in this type of landform model using digital elevation model data.

The solar radiation map models the relative level of solar energy across the Yale study area (see Map2). The solar radiation map is a continuous surface of predictive/modeled solar energy, unlike an aspect map which is a classified map of the eight cardinal directions of the compass (north, north/east, east, south/east, etc.),traditionally used in modeling as a proxy for solar radiation. The solar radiation map values match example values for similar latitudes/aspects from McCune and Keon's study (2002). However, this map should be reviewed because the Yale study area is relatively drier/hotter than other locations at similar latitudes, and therefore it could be expected that the relative solar radiation would likely be higher throughout the Yale study area. The McCune and Keon model is based solely on elevation and latitude and does not consider relative precipitation, temperature or biogeographic location.

The Geophysical Heterogeneity Index (GHI) map is a combination of three physical characteristics of the landscape: landforms (10 classes), solar radiation(reclassified to 3 classes) and flow accumulation (reclassified to 2 classes) grid maps (see Map 3). For example, a site could be classified as acombination of "Upper Slope/Mesa landform + South Facing Solar Radiation + noflow accumulation". The GHI map produced 46 unique classes (see the appendix for additional details about the values).

The Geophysical Heterogeneity Density (GHD) hexagon map aims to identify sites of relatively higher densities of geophysical heterogeneity (i.e. site of multiple physical characteristics). It was displayed by standard deviation to try and parse out hexagons with relatively higher densities of GHI classes (see map 4).

The Biophysical Heterogeneity Density (BHD) hexagon map aims to identify sites of relatively higher densities of ecosystem heterogeneity (i.e. sites of multiple vegetation types). It also was displayed by standard deviation to try and parse out hexagons with relatively higher densities of Landfire Biophysical Settings(BPS) vegetation classes (see map 5).

Map1. Landforms, 1:100,000 (a small part of the Yale study area draped over a hillshade)

Map2. Solar Radiation, 1:100,000 (a small part of the Yale study area, draped over a hillshade)

Map 3. Geophysical Heterogeneity Index (GHI) Map, 1:100,000 (a small part of the Yale study area NOT draped over hillshade)

Map 4. Geophysical Heterogeneity Density (GHD) map, displayed by standard deviation

Map 5. Biophysical Heterogeneity Density (BHD) map, displayed by standard deviation

APPENDIX:

Source Data Inputs:

• Yale Study Area boundary
• USGS NED 10 meter Digital Elevation Model

Landforms

A modified version of Jeff Jenness's Landform Classification Model (Tagil and Jenness, 2008) from Topography Tools 9.3 was used to model landforms. This tool models 10 landform classes using a digital elevation model (DEM) as source input. The 10 landform classes and their grid values (in parenthesis) include:

• canyons, deeply incised streams (1)
• midslope drainages, shallow valleys (2)
• upland drainages, headwaters (3)
• u-shapedvalleys (4)
• plains,less than 2% slope (5)
• openslopes, over 2% slope (6)
• upperslopes, mesa (7)
• localridges, hills in valleys (8)
• midslope ridges, small hills in plains (9)
• mountaintops, high ridges (10)

The USGS 10 meter NED digital elevation model was used as source data for modeling landforms in the Yale study area. This data has significant artifacts — banding/steps - that affect any derived topographic datasets (it creates significant salt and pepper in the derived grids rather than smooth continuous surfaces). It was necessary to first smooth the NED10 data using a filter to try and remove these artifacts(using a focal majority function where a moving window moves across the grid evaluating the majority value within a specified window). Smoothing will remove very high peaks and very low sinks. A 3x3 circular neighbourhood analysis window (NAW) was moved over the NED10m to remove the artifacts, and this step was then repeated second time using the first smoothed DEM as input.

Jenness Landform Classification model uses a 5% slope threshold to distinguish between Plainsand Open Slopes, but this value was modified to 2% based on expert opinion of how best to differentiate between Plains and Open Slopes landforms in the Yale study area.

Jenness landform classification tool classifies landform based on Topographic Position Index(TPI). The TPI were calculated using a moving window and is the difference between a cell elevation value and the average elevation of the neighborhood around that cell. Positive cell values meant the cell was higher than surrounding cells, while negative cell values meant it was lower.

Landform category can be determined by classifying the landscape using 2 TPI grids at different scales.The combination of TPI values from different scales suggest various landformtypes.For example, a high TPI value in a small neighborhood, combined with alow TPI value in a large neighborhood, would be classified as a local ridge or hill in a larger valley, while a low small neighborhood TPI plus a high large-neighborhood TPI would be classified as an upland drainage or depression.

Various landform classification trials were conducted using different small and large NAW for calculating the 2 Topographic Position Index (TPI) maps, to try and identify optimal small and large NAWs for developing landforms in the Yale landscape.

The large NAW generally identifies the extensive/macro scale landform features: plains, open slopes,upper slopes/mesas, and u-shaped valleys. Using the twice smoothed DEM, as the size of the small NAW increased, it increased the extent of midslope drainages/shallow valleys, upland drainages/headwaters, local ridges/hills in valleys, and midslope ridges/small hills in plains. These classes are under-represented in the landform map and therefore it was desirable to see their extents expand. But when the small NAW increased, it also significantly increased the extent of canyons/deeply incised streams, and mountain tops/high ridges, as well as upper slopes/mesas and u-shaped valleys which was not desired. The best compromise between small NAW and large NAW — where it balance pulling out the former classes, without unduly expanding the latter classes — was a small NAW of 55 and a large NAW 400. The final landform dataset was reclassified with values in the 1000s to represent each landform(i.e. 1000 = canyons, deeply incised streams. 2000 = midslope drainages,shallow valleys, etc.)

Diagram 1. Landform Classification (Jenness) ArcGIS toolbox model from Jeff Jenness's Topography Tools 9.3.

Solar Radiation

A modified version of the McCune and Keon's Solar Radiation model (McCune and Keon 2002) from Jeff Jenness's Topography Tools 9.3 was used to model solar radiation. This tool presents a GIS version of McCune's model, using a digital elevation model and a grid of latitude (decimal degrees) as source data to derive slope, aspect (folded) and latitude to model potential annual direct solar radiation (MJ/cm2/year).

It is based on McCune and Keon's (2002) equation 1 for modeling solar radiation:

Exp(-1.467 + 1.582 *COS([Latitude Radians])*COS([Slope Radians]) - 1.5*COS([Aspect Radians]) *SIN([Slope Radians]) * SIN([Latitude Radians]) - 0.262 * SIN([Latitude Radians]) * SIN(Slope Radians]) + 0.607 * SIN([Aspect Radians]) * SIN([Slope Radians]))

Where

Latitude Radians = Latitude Raster * (pi/180)
Slope Radians = Slope Degrees * (pi/180)
Aspect Radians = Folded Aspect (180 - (Aspect-180)) * pi/180

But in their paper McCune and Keon present three equations for modeling solar radiation:

Equation 1 can be used anywhere on the planet, irrespective of slope and latitude.
Equation 2, however, can be implemented at any latitude, but only on dems with slopes from 0-60.
Equation 3 is the most restrictive, and can be utilized at latitudes from 30-60 and slopes from 0-60.

The study suggests that their equation 3 produces the most robust results and that it is generally a better option. Another researcher suggested if you are in the 30-60 latitudes it is best to reclassify any slope values over 60, to 60 degrees, and therefore be able to use equation 3. There were few slopes over 60 degrees in the Yale study area, therefore I followed the suggestion and reclassified these slopes to 60 degrees, and used equation 3.

McCune and Keon's equation 3 for modeling solar radiation:

0.339 +0.808*COS([Latitude Radians])*COS([Slope Radians])-0.196*SIN([Latidude Radians])*SIN([Slope Radians])-0.482*COS([Aspect Radians])*SIN([Slope Radians])

Diagram 2. Solar Radiation(McCune 2002) ArcGIS toolbox model from Jeff Jenness's Topography Tools 9.3

The USGS NED 10 meter digital elevation model, and a 10 meter grid of latitude (decimal degrees),were used as the source data.  

Within the model any slopeover 60 was reclassified as 60 degrees (using the following CON statement(CON(Slope Degrees > 60), 60, Slope Degrees)) and replaced equation 1 with equation 3.

The resultant map is a 10 meter raster grid of potential annual direct solar radiation (MJ/cm2/year).Solar radiation is measured as a value of energy per unit area. In this case it is megajoules per cm2, per year. In the Yale study area the values ranged from 0.131 to 1.086. The solar radiation map is a continuum from least potential annual direct solar radiation (0.131 - North facing slopes) to most potential solar radiation (1.086 - South facing slopes).

North	Flat	South
0.131--------------	0.958--------------	1.086

The result showed very good correlation with McCune and Keon's example values for equation 3:

From McCune and Keon(2002), page 605:

Latitude/Slope/Aspect	McCune and Keon Example Values	Average Yale Example Values
40N, 30 slope, N aspect	0.571	0.672
40N, 30 slope, S aspect	1.053	1.055
40N, Flat	0.958	0.967

The geometrical interval classification method in ArcGIS 9.3 was used to classify the solar radiation output into three classes. The geometrical interval classification is appropriate for continuous data that is not distributed normally (for more information see http://blogs.esri.com/esri/arcgis/2007/10/18/about-the-geometrical-interval-classification-method/)

A visual review of the classified map draped over a hillshade shows that the geometric class breaks appear to be suitable thresholds for defining north aspect versus flat versus south aspect in the solar radiation grid. The solar radiation map was reclassified to three classes based on these threshold values (see map 6).

Map6. Classified Solar Radiation Map (a small part of the Yale study area, draped over hillshade)

Flow Accumulation

The USGS NED 10 meter digital elevation model was used in the ArcGIS Flow Direction/Flow Accumulation tools to model accumulated flow within the Yale study area.. Accumulated flow is the accumulated weight of all grid cells flowing into each downslope cell in the final raster.

Of note, this tool models all hypothetical stream networks based on the digital elevation model,regardless of whether a stream/river actually exists on the landscape. In this particular study area, many of the streams identified in this modeled flow accumulation map will not exist or be ephemeral.

A threshold for identifying a stream network (versus the surrounding upland) was identified in the resultant flow accumulation map based on a visual review of the data/expert opinion and the flow accumulation map was reclassified into two classes based on this class break threshold value (17).

Geophysical Heterogeneity Index (GHI) Map

The landform, classified solar radiation and classified flow accumulation maps were then combined to create a geophysical heterogeneity index map.. The landforms are represented by the 1000s values (1000 = canyon, deeply incised streams; 2000 = midslope drainages, shallow valleys; 3000 = upland drainages, headwaters; 4000 =u-shaped valleys; 5000 = plains; 6000 = open slopes; 7000 = upper slopes, mesas; 8000 = local ridges, hills in valleys; 9000 = midslope ridges, small hills in plains; 10000 = mountain tops, high ridges), the reclassified solar radiation are represented by the 100s values (100 = north facing, 200 = flat,and 300 = south facing) and the reclassified flow accumulation classes are represented by the 1s (1 = stream, 0 = non-stream)

Geophysical Heterogeneity Density (GHD) Map

A series of hexagon maps were created at various scales — 16km, 8km, 4 km, 1 km and ½ km — and the total number of GHI classes that occurred within each hexagon was calculated.

Map 7. Hexagon Map — 4 km Area, displayed over the GHI map

The distribution of the resultant Geophysical Heterogeneity Density (GHD) hexagon map was then displayed by standard deviation to identify areas of relatively higher geophysical complexity/potential climate refugia. The goal of mapping the densities at different scales was to try and find the scale at which the density per hexagon best reflected the heterogeneity on the ground.A very small hexagon would likely have occurrences of only 1 or 2 different GHI classes, whereas a very large hexagon could have occurrences of all 46 GHI classes. Both of these results are meaningless. The challenge is to identify the scale (sweet spot) somewhere in the middle of these two extremes that will model the geophysical heterogeneity of the landscape in a realistic manner.

Biophysical Heterogeneity Density (BHD) Map

A biophysical heterogeneity density map was created for comparison with the geophysical heterogeneity map, to see if this readily available vegetation map could be used as a proxy for identifying climate refugia. The Landfire Biophysical Settings map, which represents the vegetation that may have been dominant on the landscape prior to Euro-American settlement, was used as the source data for creating a biophysical heterogeneity density map. The same methodology used to create the Geophysical Heterogeneity density map was applied to produce the biophysical heterogeneity density map. The distribution of the resultant Biophysical Heterogeneity Density (BHD) hexagon map was then displayed by standard deviation to identify areas of relatively higher biophysical complexity/potential climate refugia.

Source References:

Topography Tools

Jeff Jenness.TPI_Documentation_online.pdf.

McCune, B. & Keon, D.2002. Equations for potential annual direct incident radiation and heat load. J.Veg. Sci. 13: 603-606.

Downloaded Testrad.xl. -McCune and Keon's spreadsheet with constants and equations for potential annual direct incident radiation

Sean Parks 2004. Solar Radiation ArcView Script (based on McCune and Keon, 2002) (asr_readme.txt ;)

Sermin TAGIL and JeffJenness, 2008. GIS-Based Automated Landform Classification and Topographic,Landcover and Geologic Attributes of Landforms Around the Yazoren Polje,Turkey. Journal of Applied Sciences, 8: 910-921.