**[NASA Office of Earth Science -- Applications Program]{.smallcaps}**

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

####### National Invasive Species Forecasting System {#national-invasive-species-forecasting-system .unnumbered}

##### Evaluation Report {#evaluation-report .unnumbered}

*Version 1.0*

########  {#section .unnumbered}

###### ![](media/image1.png){width="0.875in" height="0.7319444444444444in"}Science for Society {#science-for-society .unnumbered}

> ![](media/image2.png){width="1.125in"
> height="0.4173611111111111in"}*Accelerating the realization of
> economic*
>
> *and societal benefits from Earth science,*
>
> *information, and technology*

#### July 31, 2004 {#july-31-2004 .unnumbered}

# Table of Contents {#table-of-contents .unnumbered}

Executive Summary 3

1.  Introduction

    1.  Overview of Application Domain xx

    2.  NASA and Application Mission Traceability xx

    3.  Introduction to Decision Support Tool (DST) xx

    4.  Initial Partnering Activities xx

    5.  System Engineering Approach xx

2.  Overview of DST

    1.  Background xx

    2.  Detailed Description xx

3.  Consideration of NASA Inputs

    1.  Overview of NASA Inputs Matched Against Technical Requirements
        xx

    2.  Detailed Description and Analysis of NASA Inputs to Improve DST
        xx

    3.  Description of Past and Current NASA Contributions xx

4.  Gaps in Meeting DST Needs xx

5.  Conclusions and Recommendations

    1.  Findings xx

    2.  Recommendations xx

    3.  Next Steps xx

Supporting Documents

# Executive Summary {#executive-summary .unnumbered}

# 1 Introduction {#introduction .unnumbered}

The NASA Office of Earth Science and the US Geological Survey (USGS) are
working together to develop a National Invasive Species Forecasting
System (ISFS) for the management and control of invasive species on all
Department of Interior and adjacent lands. This decision support tool
(DST) will use early detection and monitoring protocols and predictive
models developed at the US Geological Survey Fort Collins Science Center
to process NASA and commercial data and create on-demand, regional-scale
assessments of invasive species patterns and vulnerable habitats useful
for a wide range of science, management, and policy decision-making
activities. The community of ecologists and land managers who will use
this tool are being involved in all stages of ISFS development. The
system will be made broadly available to the natural resource management
community through the Internet-based USGS National Biological
Information Infrastructure program.

This project builds upon NASA's unique remote sensing capabilities and
current joint USGS/NASA work on invasive species. The end product will
be a dynamic and flexible decision support system to generate electronic
and paper maps of "hot spots" for potential exotic species invasions.
These maps will be used by private and public national, state, and local
management agencies for remediation, management, and control. Over the
next 10 years, we envision that the National Invasive Species
Forecasting System will provide revolutionary capabilities of
fundamental importance to understanding and managing the effects of
invasive species on human health, the economy, native biodiversity, and
ecosystem processes. The techniques developed here will have immediate
use in other current and planned NASA Earth Science Enterprise research
programs, including mapping climate change, temporal and spatial ocean
change, disease outbreaks, and national and global carbon storage, and
will create important new opportunities for commercial development.

## Overview of Application Domain

During the past century, non-indigenous plants, animals, and pathogens
have been introduced at increasing rates into all US ecosystems. A
growing number of these species are becoming invasive, and contribute to
declines in native species, changes in ecosystem function, and
cumulative economic impacts currently estimated at more than \$137
billion annually.

An "invasive species" is defined as a non-native species whose
introduction causes or is likely to cause harm to the economy,
environment, or human health. The cost of infestations of leafy spurge
alone to agricultural producers and taxpayers is \$144 million/year in
the Dakotas, Montana, and Wyoming. Aggressive invasive fishes in the
Great Lakes threaten a commercial fishery valued at \$4.5 billion which
supports 81,000 jobs. Invasive Norway rats cause up to \$19 billion/year
in environmental and economic damage. Non-native livestock diseases cost
\$9 billion/year. In the coming decades, increasing human travel and
trade and changing types and patterns of environmental disturbance are
expected to exacerbate these impacts. Because of its high diversity of
environmental conditions and habitats, the US is particularly vulnerable
to invasions.

The Federal Government has begun to mount an organized effort to address
the invasive species threat, coalescing around Executive Order 13112
(1999). There is now a National Invasive Species Council, which has
issued a draft Management Plan, and has assembled several technical
working groups. The National Biological Information Infrastructure has
several regional programs developing invasive species information
systems as their highest priority initiative, and plans to establish a
national node for invasive species in 2002. These efforts are
coordinated with international initiatives under the United Nations (the
Global Invasives Species Programme, GISP), NAFTA (the North American
Biodiversity Information Network, NABIN), the Summit of the Americas the
Inter-American Biodiversity Information Network, IABIN), and a number of
bilateral agreements, to develop international exchange on invasive
species information. Globalization has greatly increased the
international movements of harmful species through travel and
agricultural, horticultural, and pet industries, and has become a
principal impediment to international trade agreements.

## NASA and Application Mission Traceability

All of these efforts recognize the central role of NASA's space-based
sensors and advanced computational, modeling, and information
technologies in addressing invasive species science and policy on a
national and global scale. Both the potential for movements of invasive
species, and the susceptibility of sensitive habitats to new invaders
are known to be strongly influenced by climate warming, changes in
rainfall, soil moisture, and runoff, and are increasingly driven by
extreme events. As a result, a wide range of existing and emerging NASA
Earth Science Enterprise technologies will be central to a better
understanding of invasive species risks. At the same time, many invasive
species greatly alter the water relations, carbon storage, fire cycle,
and reflectance properties of landscapes, and may be an important
feedback link to climate.

## Introduction to Decision Support Tool (DST)

High resolution mapping of biological resources is central to
confronting the invasive species threat. To meet the needs of the
invasive species science and policy communities, we must be able to
identify dominant plants and vegetation structures with a reasonable
ability to distinguish between species. This is becoming possible with
hyperspectral instruments at meter-scale resolutions, particularly when
combined with LIDAR and other active sensors that can detect meter-scale
vegetation structure, landforms, and soil surface properties. USGS uses
advanced geostatistical methods, many derived from NASA Earth Science
Enterprise support, to integrate multiple types and scales of data,
including satellite images, aerial photography, and ground data of
various resolutions in order to map biological resources.

## Initial Partnering Activities

The government partner in this activity is the USGS Geological Survey's
National Institute for Invasive Species Science located at the USGS Fort
Collins Science Center. Dr. Thomas J. Stohlgren is the ecologist lead
and is working closely with the senior management of USGS located at
USGS Headquarters in Reston, VA. Numerous partners are involved in the
USGS NIISS program.

This effort is proceeding through three major phases. The first phase
involves developing a prototype Invasive Species Forecasting System at
selected test sites, such as Rocky Mountain National Park, CO, the Cerro
Grande Wildfire Site, NM, and Grand Staircase-Escalante National
Monument, UT. Sites have chosen in consultation with the USGS
ecologists, land managers in bureaus of the Departments of Agriculture
and the Interior and private stakeholder groups. At each test site, the
project is compiling existing field survey data, ground measurements,
and airborne and satellite data. Local and regional models are being
refined and tested.

The second phase will expand the system to include MODIS time-series
data and high-resolution hyperspectral and other measurements in the
modeling protocols. Doing periodic, on-demand, national assessments of
management-scale risks will require unprecedented data integration and
computing capacity, the infrastructure for which will be developed
during this second phase.

In phase three, the system will be deployed beyond the test sites and
expanded into a fully operational National Invasive Species Forecasting
System for use throughout the Department of Interior. The system will be
made widely available to the scientific and management communities
through the USGS National Biological Information Infrastructure program,
which provides US interagency coordination for online biological
databases and information services.

## Systems Engineering Approach

{Most of this introductory material has been taken from the [EOM
article]{.underline} on the NASA/USGS Project. Becoming dated, but still
a good source for a general overview.}

# Overview of DST

{Here, in the case of ISFS, we need to emphasize that what USGS had when
we first became involved with them was a geostatistically-based
"modeling process" -- rather than a mature DSS. In order for this
process to be in a position to use more NASA data, three bottlenecks had
to be overcome: (1) improved access to the software and hardware
required to implement the process, (2) improved computational
performance for the geostatisical code, and (3) better access to
combined USGS and NASA data and a more sophisticated use of existing and
new data products. The transition from a DSProcess to a DSTool thus
requires in the case of this application three mutually dependent and in
a sense equally important NASA contributions: (1) Engineering -- to
create a web-based workflow management environment more widely
accessible to DOI/USGS partners and the extended community, (2) High-End
Computing -- to allow analysis algorithms to run fast on commodity
cluster computers (ie consume more NASA data and model larger extents),
and (3) Data Services -- to facilitate the integration of new remote
sensing measurements and data products into the USGS process.}

## Background

The ISFS system will have users from government agencies, universities,
industry, and the general public. To the users, the ISFS is deployed as
a web browser based system that will present options for applying a
series of models to available datasets yielding predictive result sets.
The system can ingest data from different sources and in different
formats and existing models can be either run, or new ones created. The
system outputs maps and additional information depicting the applied
model and the predicted species distributions.

![](media/image3.png){width="3.8868055555555556in"
height="3.467361111111111in"}

## Detailed Description

{Most of this comes from our [Baseline Software Engineering
Document]{.underline}. The language needs to be updated to reflect
current thinking. Some of the equations and symbols haven't copied
correctly but can be repaired. I'm using .jpg place holders for the
figures, but we have PostScript versions that print much better. Also
need to add figure/table captions.}

## *USGS Modeling Process*

In this section, we provide a general description of the core analysis
approach that is the focus of our development efforts. Additional
information about the geostatistical methods described here may be found
in Isaaks and Srivastava's *Applied Geostatistics* (1989). The USGS
modeling approach is currently being used to identify areas at risk for
exotic plant species invasions at the Cerro Grande Fire Site near Los
Alamos, New Mexico (Kalkhan *et al.*, 2002), in Rocky Mountain National
Park, Colorado (Chong *et al*., 2000; Kalkhan *et al*., 2000a, Kalkhan
*et al*., 2000b), and Grand Staircase-Escalante National Monument, Utah
(Kalkhan *et al*., 2000c).

## [Modeling Large-Scale Spatial Variability]{.underline} *--* As shown in Figure 2, the process begins with stepwise regression and trend surface analysis for geographical variables and measures of focal taxa to evaluate large-scale spatial variability in a study area. The functional form of this model is defined as: {#modeling-large-scale-spatial-variability-as-shown-in-figure-2-the-process-begins-with-stepwise-regression-and-trend-surface-analysis-for-geographical-variables-and-measures-of-focal-taxa-to-evaluate-large-scale-spatial-variability-in-a-study-area.-the-functional-form-of-this-model-is-defined-as .unnumbered}

![](media/image4.pct)

where, *β*~ij~ are the regression coefficients associated with the trend
surface component of the model, *γ*~k~ are the regression coefficients
associated with the *q* auxiliary variables, *y*~k0~, available as a
coverage in the GIS data base, and *η*~0~ is the error term which may or
may not be spatially correlated with its neighbors (Kallas, 1997;
Metzger, 1997).

Stepwise multiple regression analysis is used first to identify the best
linear combination of independent variables. It also allows us to
explore the variation in predicting total, exotic, and native plant
species richness as a function of the TM bands, derived vegetation
indices, tasseled cap transformation indices, slope, aspect, and
elevation. The selected independent variables are used in an Ordinary
Least Square (*OLS*) procedure to describe large-scale variability
estimates.

*OLS* estimators are used to fit the model if the variable of interest
has a linear relationship with the geographical coordinates of the
sample plots, the digital number (DN) value of any of the Landsat TM
bands, and the topographic data. In addition, the least squares method
fits a continuous, univariate response as a linear function of the
predicted variable. This trend surface model represents continuous first
order spatial variation. Akaike's Information Criteria "AIC", (Brockwell
and Davis 1991, Akaike 1997) is used as a guide in selecting the number
of model parameters to include in the regression model where:

AIC = -2 (max log likelihood) + 2 (number of parameters) (2)

When using maximum likelihood as a criterion for selecting between
models of different orders, there is the possibility of finding another
model with equal or greater likelihood by increasing the number of
parameters (Metzger 1997). Therefore, the AIC allows for a penalty for
each increase in the number of parameters. Using this criterion, a model
with a smaller AIC is considered to have a better fit. While, the model
is kept as simplistic as possible, a more complex model could be used if
the situation warrants it.

![](media/image5.png){width="5.315972222222222in"
height="3.9868055555555557in"}

[Modeling Small-Scale Spatial Variability]{.underline} -- In the next
stage of the model building process the residuals from the trend surface
models are analyzed for spatial dependencies. This is accomplished using
spatial auto-correlation and cross-correlation statistics. If the
residuals are cross-correlated with other variables, we can use
co-kriging to interpolate the residuals. However, if the residuals are
not cross-correlated, we use ordinary kriging. Finally, the weights
associated with the kriging and co-kriging models are estimated as a
function of the spatial continuity of the data (Isaaks and Srivastiva
1989). This estimation can be accomplished using a sample variogram to
describe spatial continuity. With spatial data, the variation of the
samples generally changes with distance. In other words, the variogram
is a measure of how the variance changes with distance. The variogram
and cross-variogram models used in this analysis were considered "basic"
models, meaning they are simple and isotropic (Reich et al. 1999). They
include, Gaussian, spherical, and exponential models. Since the primary
focus of our work in this project involves improving the performance of
kriging, we now describe this technique in some detail.

Kriging is a method of interpolation named after a South African mining
engineer named D. G. Krige who developed the technique in an attempt to
more accurately predict ore reserves. Over the past several decades
kriging has become a fundamental tool in the field of geostatistics.
Kriging is based on the assumption that the parameter being interpolated
can be treated as a regionalized variable. A regionalized variable is
intermediate between a truly random variable and a completely
deterministic variable in that it varies in a continuous manner from one
location to the next and therefore points that are near each other have
a certain degree of spatial correlation, but points that are widely
separated are statistically independent (Davis, 1986). Kriging is a set
of linear regression routines which minimize estimation variance from a
predefined covariance model.

The first step in ordinary kriging is to construct a variogram from the
scatter point set to be interpolated. A variogram consists of two parts:
an experimental variogram and a model variogram. Suppose that the value
to be interpolated is referred to as f. The experimental variogram is
found by calculating the variance (g) of each point in the set with
respect to each of the other points and plotting the variances versus
distance (h) between the points. Several formulas can be used to compute
the variance, but it is typically computed as one half the difference in
f squared.

![](media/image6.jpeg){width="4.777777777777778in" height="3.0in"}

Once the experimental variogram is computed, the next step is to define
a model variogram. A model variogram is a simple mathematical function
that models the trend in the experimental variogram. As can be seen in
the above figure, the shape of the variogram indicates that at small
separation distances, the variance in f is small. In other words, points
that are close together have similar f values. After a certain level of
separation, the variance in the f values becomes somewhat random and the
model variogram flattens out to a value corresponding to the average
variance.

Once the model variogram is constructed, it is used to compute the
weights used in kriging. The basic equation used in ordinary kriging is
as follows:

![](media/image7.png){width="0.93125in" height="0.3486111111111111in"}

where n is the number of scatter points in the set, f~i~ are the values
of the scatter points, and w~i~ are weights assigned to each scatter
point. The weights used in kriging are based on the model variogram. For
example, to interpolate at a point P based on the surrounding points
P~1~, P~2~, and P~3~, the weights w~1~, w~2~, and w~3~ must be found.
The weights are found through the solution of the simultaneous
equations:

![](media/image8.png){width="2.43125in" height="0.19791666666666666in"}

![](media/image9.png){width="2.43125in" height="0.19583333333333333in"}

![](media/image10.png){width="2.43125in" height="0.20694444444444443in"}

where S(d~ij~) is the model variogram evaluated at a distance equal to
the distance between points i and j. For example, S(d~1p~) is the model
variogram evaluated at a distance equal to the separation of points P~1~
and P. Since it is necessary that the weights sum to unity, a fourth
equation:

![](media/image11.png){width="1.30625in" height="0.15763888888888888in"}

is added. Since there are now four equations and three unknowns, a slack
variable, l, is added to the equation set. The final set of equations is
as follows:

![](media/image12.png){width="2.68125in" height="0.19722222222222222in"}

![](media/image13.png){width="2.68125in" height="0.1986111111111111in"}

![](media/image14.png){width="1.625in" height="0.16458333333333333in"}

The equations are then solved for the weights w~1~, w~2~, and w~3~. The
f value of the interpolation point is then calculated as:

![](media/image15.png){width="1.63125in" height="0.17777777777777778in"}

By using the variogram in this fashion to compute the weights, the
expected estimation error is minimized in a least squares sense. For
this reason, kriging is sometimes said to produce the best linear
unbiased estimate (BLUE). However, minimizing the expected error in a
least squared sense is not always the most important criteria and in
some cases, other interpolation schemes give more appropriate results
(Philip & Watson, 1986).

An important feature of kriging is that the variogram can be used to
calculate the expected error of estimation at each interpolation point
since the estimation error is a function of the distance to surrounding
scatter points. The estimation variance can be calculated as:

![](media/image16.png){width="2.50625in" height="0.2111111111111111in"}

When interpolating to an object using the kriging method, an estimation
variance data set is always produced along with the interpolated data
set. As a result, a contour or iso-surface plot of estimation variance
can be generated on the target mesh or grid.

Since kriging is a rather complex interpolation technique and includes
numerous options, a complete description of kriging is beyond the scope
of this report. The reader is strongly encouraged to consult the UNCERT
User Guide (Wingle, et.al, 1995) and the GSLIB textbook (Deutsch and
Journel, 1992) for more information. Other good references on kriging
include Royle et. al. (1981), Davis (1986), Lam (1983), Heine (1986),
Olea (1974), Journel & Huijbregts (1978). Isaaks and Srivastiva's (1989)
chapter on "Ordinary Kriging" is particularly helpful.

## *Modeling Workflow*

The practical implementation of the modeling process, as developed and
used by colleagues at USGS, consists of a pre-processing step, a
modeling step, and a post-processing step.

[Pre-Processing]{.underline} -- Pre-processing activities merge ingested
datasets to create a data product that can be analyzed in the subsequent
modeling step. In the baseline scenario, the field data are merged with
the Landsat and DEM information at the same UTM x,y coordinates.
Resampling may be performed at this time if the input data are not at
the same resolution, and the Landsat data may be processed to higher
level products, e.g. tassel cap coefficients, principal components,
atmospherically corrected reflectance values, etc. The merged data
product is written to backing store in a common analysis format.

For the most part, this pre-processing step is a straightforward
application of common techniques, and is thus not a major focus of our
current work. A possible exception to be explored in the future is the
pre-processing that might be needed to accommodate new data sets in the
model, e.g., atmospherically corrected hyperspectral data which presents
different types of computational challenges.

[Modeling]{.underline} -- The primary modeling pipeline uses the merged,
flat file resulting from the pre-processing step. The file is logically
arranged with one row of data for each field survey point. The data
sampled at each point are arranged in columns. The file contains a
simple internal ASCII header that contains the number of rows and
columns along with a one-word column descriptor. The columns include a
subset of the following data: location, plant, soil, digital elevation
model (DEM), and remote sensing information such as Landsat DN and
derived quantities (e.g. tassel cap coefficients or NDVI values). A
series of statistical analyses are then performed in S-plus as follows:

1)  Read the input field data to create an object within S-plus.

2)  Compute the distance matrix, which is the Euclidean distance between
    each sample point.

3)  Perform a stepwise multiple regression with total plants as the
    dependent variable and the DEM, remote sensing, soils, etc. data as
    the independent variables.

4)  Perform a weighted ordinary least squares (OLS) fit to the total
    plants for the independent variables that are found to be
    significant predictors.

5)  Compute Moran's I coefficient to determine if there is spatial
    structure in the residuals of this OLS fit.

6)  If there is no spatial structure then skip to step 10.

7)  Compute the variogram of the residuals to determine the spatial
    structure.

8)  Determine whether a gaussian, exponential, and spherical model best
    fits the variogram.

9)  Perform kriging to estimate the residual surface across the entire
    study area. The kriging can be performed either in S-plus or using a
    Fortran program.

    a.  If S-plus kriging is performed:

        i.  The kriged residual surface is created directly as an S-plus
            object. Error estimates are also calculated.

        ii. The kriged residuals and estimated errors are then written
            to separate ASCII files. These files contain headers listing
            the number of rows and columns in the kriged surface along
            with georeferencing information.

    b.  If Fortran kriging is performed:

        i.  The residuals at each field sample point are written to an
            intermediate ASCII file, along with a header containing
            parameters that control the kriging. These parameters
            include the number of rows and columns to be kriged, the
            spatial resolution and georeferencing information for the
            kriged surface, the number of field sample points, the
            number of nearest neighbors to include in the estimation,
            the parameters describing the variogram model (e.g. range,
            nugget, sill, gaussian), and finally a flag indicating
            whether to calculate the error estimates for the kriged
            residuals.

        ii. The compiled and linked Fortran executable is invoked. It
            reads the above file of residuals, performs the kriging, and
            writes the kriged residuals to an output ASCII file. The
            format of this file is identical to that created if S-plus
            kriging is performed. The error estimates for the kriged
            residuals are written to a separate file using the same
            format.

    c.  The ASCII files with the results of the kriging are converted to
        a binary raster format using a simple filter. Separate header
        files are created so that the kriging results can be easily
        viewed the using version 3.5 of ENVI, the Environment for
        Visualizing Images from Research Systems, Inc. (see
        [http://www.rsinc.com](http://www.rsinc.com/)).

10) Apply the results of the OLS fit to the maps of the significant
    independent variables to create an estimate of the spatial
    distribution of the total plants.

11) If the residuals were kriged in step 9, then add the kriged residual
    surface to the estimate of the total plants.

[Post-Processing]{.underline} -- The post-processing step applies the
results of the above modeling activities to generate products such as
the map shown in Figure 5. In the canonical case, we:

12) Create a JPEG rendering of the total plants map. Steps 10 -- 12 are
    performed using version 5.5 of IDL, the Interactive Data Language
    from Research Systems, Inc.

13) Create separate header files so that the total plants and error
    estimates can be easily read using ENVI.

![](media/image17.jpeg){width="4.754861111111111in" height="3.675in"}

## *Baseline Performance Evaluation*

In this project, we are working with three "canonical" study sites: the
Cerro Grande Fire Site in Los Alamos, NM (CGFS), Rocky Mountain National
Park, CO (RMNP), and Grand Staircase Escalante National Monument, UT
(GSENM). The three sites provide contrasting ecological settings and
analysis challenges and vary in the types and scales of data used, areas
covered, and maturity of the investigation.

Initially, the most comprehensive modeling efforts involved the Cerro
Grande site which covers approximately 50,000 acres. Cerro Grande
modeling activities integrate a range of environmental attributes and
data sets including, as described below, data from 79 field sample
plots. The Rocky Mountain study site, at nearly 25,000 acres, is a
smaller study area but includes data from nearly 200 field sample sites.
RMNP has been studied the longest of the three sites. The newest study
area is Grand Staircase Escalante National Monument, which at 1.9
million acres is by far the largest investigation to be undertaken with
these modeling techniques. The gathering of field data started three
years ago and continues with over 350 plots being studied to date.

Each study site is being used to examine various aspects of the analysis
approach being refined here. Cerro Grande provides the best opportunity
to analyze baseline performance characteristics of the modeling system.
RMNP is a rich and mature dataset, but appears to show only weak spatial
anisotropy for the dependent variables currently being considered,
making it a poor candidate for kriging. GSENM will be a particularly
abundant dataset, but field work for GSENM will not be completed until
the end of this summer. At later stages of the project, both RMNP and
GSENM will be important for evaluating new and larger applications of
the modeling system.

[Cerro Grande Fire Site]{.underline} -- Investigating spatial
relationships among fuels, wildfire severity, and post-fire invasion by
exotic plant species through linkage of multi-phase sampling design and
multi-scale nested sampling field plots, pre- and post-fire, has been
accomplished on the CGFS. The technique provides useful information and
tools for describing ecological and environmental characteristics
including landscape-scale fire regimes, invasive plants, and hot spots
of diversity (native and non-native plants) for the site.

The Cerro Grande Fire Site is located near Los Alamos, New Mexico with
elevation range from 1932 to 3200 meters. The Cerro Grande fire began as
a prescribed fuel treatment by Bandelier National Monument, Los Alamos,
NM on 4 May 2000. The fire escaped control and was declared a wildfire
on 5 May 2000. The fire was contained on 24 May after burning about
19,300 ha of lands managed by seven different agencies, including the
town of Los Alamos, NM. However, 60% of the fire area burned 10-11 May,
2000, and 60% of the fire was on the Española District of the Santa Fe
National Forest (Burned Area Emergency Rehabilitation \[BAER\] Team
2000). Initial remotely sensed estimates of burn severity were
classified as high (35%), moderate (9%), and low (56%).

To predict the distribution, presence, and patterns of native and exotic
species, we used 79 data points (based on Modified-Whittaker nested
plots of 1000 m^2^) to represent different variables that were extracted
from Landsat TM data (eight bands, six vegetation indices, and six bands
of tasseled cap transformation indices), topographic data (elevation,
slope, and absolute aspect), and vegetation characteristics. Spatial
statistics were used to integrate these data to model large- and
small-scale variability.

In the canonical case, we use trend surface models that describe the
large-scale spatial variability using stepwise multiple regressions
based on the Ordinary Least Squares (*OLS*) method. Models with small
variance were selected. In addition, the residuals from the trend
surface model based on the *OLS* estimates were modeled using ordinary
kriging for modeling small-scale variability based on a Gaussian
semi-variogram. The final surfaces were obtained by combining two models
(the trend surface based on the *OLS* and the kriging surface of
residuals). All models were selected based on lowest values of standard
errors, AICC statistics, and high *R^2^*. For large-scale spatial
variability models using the *OLS* procedure, *R2* values ranged from
10.04% to 58.6% and all variables were significant at α \< 0.05 level.
When adding the kriging model with the OLS model, *R2* values ranged
from 60% to 84%.

[Baseline Hardware and Software Environment]{.underline} -- The baseline
processing system uses a combination of COTS and public domain software
to generate maps of estimated biodiversity or ecosystem parameters. The
major software components include the following:

-   S-plus version 6.0.1 for Linux has been used for the baseline test
    cases. S-plus is a commercially available statistical package the
    common use in many scientific communities. The functionality of
    S-plus has been enhanced by a large collection of spatial
    statistical functions developed and maintained by Drs. Robin Reich
    and Richard Davis of the Colorado State University (see
    <http://www.cnr.colostate.edu/~robin/>).

-   Version 3.5 of ENVI, the Environment for Visualizing Images from
    Research Systems, Inc.

> [http://www.rsinc.com](http://www.rsinc.com/)). ENVI is a commercially
> available image analysis application in common use in the remote
> sensing community.

-   Version 5.5 of IDL, the Interactive Data Language from Research
    Systems, Inc. IDL is a commercially available image analysis
    application in common use in the remote sensing community.

The baseline processing hardware is a singe-processor AMD machine
running RedHat Linux 7 (kernal 2.4.9-31) at 1.2 GHz with 1.5G of RAM and
60G of disk space.

[Computational Performance Characteristics]{.underline} --Three factors
influence the performance of the model: the size of the output surface
area over which kriging occurs (area), the total number of number of
sample points in the data set (pts), and the number of "nearest
neighbor" (nn) sample points from the total data set actually used to
compute a kriged value for any given point in the output area. When we
first began work with colleagues at USGS, a scalar, single-processor run
of this model using S-plus took approximately two weeks. The major
computational bottleneck in the model is the kriging routine. Solving
for the weights in the equations that form the ordinary kriging system
(Eq. 6) uses LU decomposition with backsubstitution to do matrix
inversions. The overall computational complexity of ordinary kriging is
thus O(*n*^3^), and the time required to compute a result is strongly
influenced by the number of sampled data points used to estimate the
residual surface across the entire study area.

In order to achieve the two-week result described above, initial USGS
model runs limited kriging (and thus the size of the computed covariance
matrices) to only 18 out of the total 79 sample points for Cerro Grande.
The kriging process iterates over the rows and columns of the output
surface. For each (i,j) point of the output surface area, 18 nearest
neighbor sample points were found and ultimately transformed into the
appropriate weighted average to estimate the kriged value at point
(i,j). This sub-sampling of 18 nearest neighbors is significant because
it represents an accommodation that may be appropriate and exploitable
in some circumstances while other types of applications may require the
use of significantly more sampled points or the entire set of sampled
points. The implications of these options will be explained in more
detail below.

The output surface representing the entire CGFS area consists of 652 x
715 points, and the canonical CGFS dataset consists of 79 data points.
Appendix B and Figures 6a and 6b show the results of baseline runs of
the CGFS case and confirm that processing time scales linearly with both
the size of the output surface (area) and the number of total data
points (pts). In contrast, processing time increases order *n*^3^ with
respect to the number of nearest neighbors (nn) being used in the
kriging routine (Figure 7). In order to do our baseline run, a FORTRAN
kriging routine was developed that ran approximately two orders of
magnitude faster than the original S-plus routine (see Appendix D).
While the overall project of growing this analysis approach into a
comprehensive Invasive Specie Forecasting System involves many elements,
kriging will continue to be the focus of our efforts to improve model
performance.

![](media/image18.jpeg){width="5.0in" height="2.7159722222222222in"}

![](media/image19.jpeg){width="5.04375in" height="2.754166666666667in"}

![](media/image20.jpeg){width="5.5in" height="2.988888888888889in"}

# Consideration of NASA Inputs

## Overview of NASA Inputs Matched Against Technical Requirements

## Detailed Description and Analysis of NASA Inputs to Improve DST

## *Engineering*

{I think most of this material can come from our [Software Design
Document]{.underline}.}

## *High-End Computing*

{This also is from our Baseline Software Design document. Shows our
approach to algorithm speed ups -- but the language/timings, etc. need
to be updated and contextualized for the evaluation report. We're
progressed through this almost to completion and can update this or
other benchmarking documents with improvement results.}

[Parallelization Strategy]{.underline} -- The overall goal for code
improvement is to reduce processing times and increase the amount of
data handled by the model. As described below, increasing the amount of
data handled by the model translates into either increasing
spatiotemporal resolution or increasing coverage (Table 6.1). We first
wish to accomplish quantitative improvements in the underlying model
that have been agreed upon by the user community as minimal advances
needed to improve core capabilities. These goals were driven by the
knowledge that in Year 3 of the project we will build a 32-node cluster
in the USGS facility. We refer to these as "Community Improvement
Goals." The ESTO/CT program, however, provides access to even greater
computational capabilities that can be used to apply this modeling
approach to some important and challenging problems that have heretofore
been unapproachable. We would therefore like to use CT\'s clusters to
attain more challenging performance improvement goals at the same time
we are accommodating basic needs. We refer to these complimentary
challenges as "Advanced Improvement Goals." Each will be described
below.

![](media/image21.png){width="5.623611111111111in"
height="4.519444444444445in"}

[Milestone 1 -- First Code Improvement (Parallelization):]{.underline}
*Improve implementation of the model to deliver canonical products
[m]{.underline}X faster than the baseline implementation.*

For Milestone 1, we will focus on the community goal of increasing speed
by at least a factor of 25. This corresponds to "*m* = 25" in our
Milestone F Goal Statement. Our approach will be to develop parallel,
Message Passing Interface (MPI) versions of the model. In the model, the
kriged estimate for each pixel of the output surface does not depend on
neighboring pixels. The kriging process is thus spatially independent,
and we expect to achieve good scaling performance by domain
decomposition. A parallel implementation of kriging would have three
steps. The first step would broadcast input data from a "control node"
to all other nodes of a cluster. The second step would have each node
compute its piece of the kriged surface. For example, node 1 could
process rows 1 -- 16, node 2 could process rows 17 -- 32, etc. The final
step would assemble the entire kriged surface on the control node. The
parallel second step is maximally efficient; the first and last serial
steps are the overhead costs of a parallel implementation.

Achieving a 25x speed up will require approximately 75% scaling
efficiency on a 32 node cluster, or its equivalent on a larger
cluster[^1]. This efficiency requires that we krig the CGFS residuals in
about 3.0 minutes. A completely efficient implementation would give a
run time of about 2.5 minutes. We must therefore incur no more than 30
seconds of parallel overhead to meet our goal. We feel this is
attainable since even a serial broadcast of the input data using 'scp',
followed by a serial gather of the individual pieces of the kriged
surface, can be accomplished in this time on a 10base-T network. A
cluster built with 2Gbps Myrinet and using much lower overhead
communications calls should be able to perform substantially better.
Further, it should be possible to hide nearly all the step three
communications overhead by overlapping the computation of the current
row with the sending of results from the previous row to the control
node. In this approach each processor would calculate the first row,
asynchronously send a message to the control node with the results of
the first row, then immediately continue processing the second row.
Synchronization would only be required at the end of the second row to
check that the results of the first row had been received successfully.
This approach should be efficient since only a small number (\~31) of
fairly large messages (\~2.6 kB) need to be sent each time a row has
been calculated (\~5.7 seconds).

[Milestone 2 -- Second Code Improvement (Adaptive Kriging):]{.underline}
*Improve implementation of to accommodate 10X more input data over
baseline condition and Milestone 1 at [n]{.underline}X the time required
in the baseline implementation. (Depending on the science problem, this
enhanced capability may be used to increase spatial resolution, temporal
resolution, or coverage.*

For Milestone 2, we will focus on the community goal of increasing by a
factor of 10 the amount of data that the model is able to process. There
are two ways that the model can be construed to accommodate more data:
the model can either handle more data points in its kriging routine or
it can krig data over a larger area. For both cases, we propose to
combine the improvements achieved through parallelization with an
adaptive approach to sub-sampling datasets.

In natural systems, spatial processes have a finite range of influence.
As the number of sample data points grows (say from 79, in the case of
CGFS, to 10x = 790), we do not necessarily wish to scale the kriging
procedure to the full set of sampled data. (In fact, if the entire data
set is used, computing the covariance matrix need only be done once,
vastly reducing the computational complexity of the kriging task.)
Ideally, we would rather scale on the basis of the number of
spatially-relevant nearest neighbors in the region we wish to estimate.
In many cases, this will involve only a fraction of a larger dataset. In
the original applications of the model, USGS arbitrarily kriged using
only 18 of 79 sample points based on a simple notion of nearest neighbor
Euclidian distance of sample points from points of the output surface.
We propose to refine these techniques by using a "sliding window" that
at each point of the output surface adaptively selects a subset of
sample points for kriging based on statistics or on an understood
spatial influence on the physics or biology of the dependant variables
being examined. We refer to this approach as "adaptive kriging."

Intelligent minimization of the number of nearest neighbors can
significantly improve the model's ability to handle larger datasets.
Perhaps more important, it provides a context for exploring mechanistic
aspects of the modeling problem that are elided by the overarching
statistical approach. Assuming that adaptive kriging uses 10% of the
total number of sampled points, a 10x increase in dataset size
essentially returns us to the baseline condition where we would expect a
25x speed-up through the Milestone F parallelization. Since the problem
scales linearly with area, performing adaptive kriging over a ten-fold
larger output area would result in an expected speed-up of approximately
2.5x. These cases correspond respectively to "*n* = 25" and "*n* = 2.5"
in our Milestone G Goal Statement.

It is important to note that adaptive kriging using more data over a
fixed-size area essentially increases model spatial resolution; adaptive
kriging using a fixed-size dataset over a larger output area increases
model spatial coverage; in both cases, repeated runs of adaptive kriging
using time-series data increases the model's temporal resolution. All
three classes of improvement have been identified by the research
community as needed enhancements to the science and technology
underlying biotic prediction.

[Advanced Applications]{.underline} -- The scenarios described here
represent only a starting point and a means of documenting quantitative
improvements in the underlying model. Along the way, and not explicitly
a part of the project's formal milestones, we also hope to apply the
model to some important and challenging problems that have heretofore
been unapproachable. For example, an important component of the model is
its consideration of the source of error in the application of the
model. We would like to use Monte Carlo simulation methods to examine
the effects of error propagation, which should become possible with next
year's work. We are also developing scenarios where the system can be
used for interactive "what-if" explorations of large datasets. Finally,
and perhaps most important, we are positioning the project to
accommodate a new suite of data sources that will become available soon.
County-level plot data are currently being gathered from a range of
sources for the entire state of Colorado (66.7 million acres), and
decadal-scale time series data on the spread of several invasive species
in the Southwestern US (500 million acres) are being assembled. Within
the next year, we believe it will be possible to use 512-node and larger
clusters running this modeling system to take on problems three to four
orders of magnitude larger than the canonical cases.

![](media/image22.png){width="6.124305555555556in"
height="2.4243055555555557in"}

## *Data Services, Measurements, and Models*

{Much of this can come from out Software Design Document, the Schnase's
ISDS proposal, and Morisette's MODIS Time-Series proposal.}

## Description of Past and Current NASA Contributions

{Basically, it goes like this:

[Engineering]{.underline}

-   Predicting Regional-Scale Exotic Plant Invasions in Grand
    Staircase-Escalante National Monument (Schnase, NASA Carbon NRA,
    10/01-9/04, \$715K)

[High-End Computing]{.underline}

-   Biotic Prediction: Building the Computational Technology
    Infrastructure for Public Health and Environmental Forecasting
    (Schnase, NASA ESTO/CT CAN, 2/02-2/05, \$1.1M)

[Data Services, Measurements, and Models]{.underline}

-   The Invasive Specie Data Service: Towards Operational Use of ESE
    Data in the USGS Invasive Species Decision Support System (Schnase,
    NASA REASoN CAN, 12/03-12/08, \$3.1M)

-   Value Added Products from MODIS Time-Series Data Sets to Support
    DOI/USGS Invasive Species Management (Morisette, NASA EOS NRA,
    12/03-12/06, \$250K) }

# Gaps in Meeting DST Needs

# Conclusions and Recommendations

The National Invasive Species Council has noted, "no comprehensive
national system is in place for detecting and responding to incipient
invasions." Yet the threat of invasive species is perhaps our most
urgent economic and conservation challenge. There is a growing sense
among land management agencies that a national assessment of native and
non-indigenous plant diversity needs to be completed on all public
lands. Especially high on the agenda are issues such as detecting the
loss of native plant diversity caused by non-indigenous plant species,
predicting where non-indigenous species are most likely to damage native
diversity so that management can be targeted at the most vulnerable
areas, and developing a science-based long-term monitoring plan for
vegetation and soil resources. The work proposed here addresses concerns
common to most land managers, and NASA has an essential role to play
helping USGS build this capacity.

## Findings

## Recommendations

## Next Steps

## Technology Plan

##  {#section-1 .unnumbered}

[^1]: As part of this project, we will build and deploy a cluster of at
    least 32 nodes at USGS specifically to support these modeling
    activities.