\begin{article} %% DOCSTART

\title{Using GRASS and R for Landscape Regionalization through PAM Cluster Analysis}
\author{Allan D. Hollander}
% Keep this next line to generate the title page.  Don't worry, it will fit into our master document well, even if it creates a separate title page at the moment.  A date will automatically be shown as well, don't worry about it.
\maketitle

\section{Abstract}
Landscape regionalization is a frequently encountered need in the geographical
sciences, having applications ranging from sampling design to conservation
prioritization. One technique for partitioning the landscape is to use cluster
analysis of GIS layers describing the area under study. Here I present a GIS
technique that uses partitioning around medoids as its clustering algorithm.
Partitioning around medoids (PAM) is a non-hierarchical clustering algorithm
that is related to the commonly-used k-means clustering technique. PAM differs
from the k-means algorithm in that a) PAM assigns cluster centroids to actual
data observations, rather than using values averaged over subsets of the entire
dataset and b) PAM accepts categorical data as input in addition to numerical
data. These properties of PAM make it useful for landscape regionalization
because often one wishes to incorporate categorical variables such as vegetation
class or soil types in the regionalization. I illustrate the PAM technique with
an example of sampling design for a local-scale analysis of agroecosystems in
Northern California. For this work I use GRASS and R, generating in GRASS a set
of random points covering the study area and attributing these points with
values from raster and vector layers of interest, importing this data table into
R for the PAM cluster analysis, and exporting the resulting clusters back to
GRASS for geographic visualization.

\section{Background} % Sections automatically end when a new one is started.
A long-standing concern in geography has been regionalization, or how to define
areas in space that share similar characteristics according to some theme of
interest. Quantitative work on regionalization methods has been undertaken for a
considerable period of time (e.g. \cite{spence1970}, \cite{johnston1970} and
continues to the present (see \cite{duque2007} for a recent
review). In the earth and environmental sciences, regionalization finds
application in areas as diverse as hydrology (\cite{wiltshire1986}), biogeography
(e.g. \cite{patten2008}, \cite{proche2005}), and climatology (\cite{stooksbury1991}). 
Although no quantitative method can provide a single definitive regionalization suitable 
for all applications, it is useful to have a suite of
formal methods available when a problem calling for defining a new set of
regions arises. In what follows I present a simple method for generating regions
for use in a GIS analysis, working with the open source GIS software GRASS
(\cite{grass2009}) and the statistical package R (\cite{r2009}).

The particular focus of this approach is to use cluster analysis methods to
support studies in landscape ecology and in conservation. In conservation, one
class of problems that can be addressed by cluster analyses is reserve design,
where a classical problem is to ensure that a reserve network contains
sufficient representation of all elements of conservation interest (e.g.
\cite{vane1991}). A cluster analysis approach can serve to capture
dissimilarities in the environment by placing these in different clusters and
hence provide a means to adequately represent the biological diversity of a
landscape (e.g. \cite{trakhtenbrot2005}). Another class of problems is
sampling design. In ecological inventories, the number of sampling plots one can
place is always quite limited. A usual strategy to increase the efficiency of
one's sampling process is stratified design: that is, rather than placing
sampling points randomly across the study region, one partitions the landscape
across environmental variables of interest, and subsequently places similar
numbers of sampling points within each partition. Cluster analysis lends itself
well to defining the strata, since it provides a formal method to separate the
variability in the environment into different groups.

Cluster analysis algorithms can be divided into hierarchical algorithms and
partitional algorithms (\cite{kaufman2005}). In hierarchical algorithms,
clusters are constructed sequentially to form a nested hierarchy of cluster
types. In partitional algorithms, the set of clusters is constructed all at
once, often requiring establishment of the desired number of clusters
beforehand. In general, cluster analysis proceeds in two steps, the first being
creation of a dissimilarity matrix, where one uses a metric to compute the
distance in attribute space between a pair of observations and assigns this
value to the element of the matrix where its row and column is indexed by the id
numbers of the observation pair. The second step is to take the distances in the
dissimilarity matrix and apply an iterative algorithm to sort the observations
into different clusters. 

In landscape analysis, often the environmental properties of interest are
categorical variables. Examples of these include habitat types or soil
classifications. The fact that some of the variables are categorical poses
problems for many clustering algorithms, for instance the most commonly-used
partitional clustering algorithm, k-means clustering. In k-means clustering, the
number of clusters is set before running the algorithm. and the algorithm is
initialized by randomly placing that number of cluster seeds among the
multidimensional environmental attribute space. The algorithm then adjusts the
position of these cluster seeds so as to minimize the dissimilarities within a
cluster and maximize the dissimilarities between clusters. The k-means algorithm
works only with continuous variables, and the cluster centroids are defined to
be the mean values of the attributes of all the observations in the cluster.
Since this algorithm relies upon computation of mean values, it fails to work
with categorical variables, and an alternative algorithmic strategy needs to be
taken if one wishes to incorporate categorical variables.

There are two elements in this strategy. The first is to use a dissimilarity
metric that works with both categorical and numerical variables. One such metric
is Gower's similarity coefficient (\cite{gower1971}). This coefficient computes the
distance between two observations by averaging a difference metric over both the
numerical variables and the categorical variables. This difference metric ranges
from 0 to 1 for both types of variables. For the numerical variables, the
difference metric is the absolute value of the ratio of the difference of the
values of the two observations and the range of the variable across all
observations. For categorical variables, the difference metric is 0 if the two
observations share the same categorical value and 1 if they differ. The second
element in this strategy is using a clustering algorithm that works with mixed
variables. The algorithm used here is partitioning around medoids (PAM) (\cite{kaufman2005}). 
This algorithm is similar to k-means clustering but differs
from it in that the cluster centers (the ``medoids'') are actual observation
points rather than centroids averaged over a number of observations. The medoids
therefore provide the sample point that is the exemplar for each cluster.

This approach can easily be carried out using the open source software programs
R and GRASS. R is excellent for testing novel statistical techniques since
academic statisticians nowadays develop many new algorithms in the R statistical
environment. GRASS provides excellent integration of raster and vector GIS
capabilities and scripting abilities.

\section{Methodology}
Figure \ref{fig:workflow} illustrates the workflow in this analysis. The diagram at upper left
shows the stack of GIS layers together with the random points to be overlaid on
the GIS stack. This creates a data table which is fed into the R statistical
environment at upper right. At lower right the clusters are constructed, which
are then mapped in the GIS at lower right. Below I detail the steps one takes to
carry out this workflow.

\begin{figure*}
  \centering
   \includegraphics[scale=1.0]{hollander/pamworkflow2} % Point to the graphics file
   \caption{PAM cluster workflow} % Provide a caption for the figure
   \label{fig:workflow} % change figure_label to a short description that you will use later on
\end{figure*}

\begin{enumerate}

\item The first step is to assemble the data layers to be used in the analysis in
GRASS. Both raster and vector data can be used. It is important to recognize
which datasets are categorical (e.g. land cover) and which are numeric (e.g.
soil percent organic matter content), since treatment of these will vary in the
statistical analysis. For vector data there may be several fields of interest
within the attribute table of a single data layer.

\item The sampling density needs to be determined. There are no fast rules for
doing this, but one needs to pay attention to the spatial resolution of the
datasets being queried. Given a range of spatial resolutions of the datasets,
there is no point in sampling at the finest resolution, but sampling only at the
coarsest resolution will by definition lose a lot of the detail in the rest of
the datasets. A choice somewhere in the middle of the range may be the most
satisfying choice. For instance, in an application where many of the data
variables come from a 1-kilometer gridded soils map, one might choose a sampling
density of 4 points per square kilometer.

\item The routine for generating random points in GRASS, v.random, asks for a fixed
number of points to be generated within the current geographic region. Computing
that value from a desired sampling density requires calculating the area of the
region from its bounding box and dividing through by the sampling density.

\item Often one wishes to restrict sampling to a subregion within the rectangular
geographic region. In GRASS this may be accomplished by defining the subregion
as a vector polygon layer, and then selecting the points that fall within that
polygon using the \texttt{v.select} command.

\item Next one makes a list of the variables to be sampled and adds these as
columns to the attribute table of the sample point layer using the command
\texttt{v.addtable}, paying attention to the data types of the variables when one creates
the columns.

\item These columns then get instantiated with data from the layers of interest,
using the commands \texttt{v.what.rast} and \texttt{v.what.vect}. Running these commands over a
large stack of layers can get tedious but this routine can easily be scripted.

\item The attribute table for the sample points is then exported so that it can be
read into the R statistical environment. There are several options for export
formats; of these a good choice is dbf since it directly preserves data type
information of the data table. To export dbf, one selects the dbf option from
the \texttt{db.out.ogr} command.

\item One then switches over to the the R environment. Though support for spatial
datasets in R is well-developed,\footnote{See
\url{http://cran.r-project.org/web/views/Spatial.html}} which includes a direct
interface to GRASS through the spgrass6 package), the steps that follow do not
depend on this interface, and hence one does not need to understand how R treats
spatial data in this workflow. The first step in the R environment is to read
the sample point data using the \texttt{read.dbf()} function which is part of the foreign
package.

\item One then checks to see that the variables in the data frame have the correct
data type. Categorical variables may have been represented as integers in this
table, especially if they have been sampled from a raster, for example one
displaying discrete land cover classes. Such variables need to be reassigned to
categorical type using the \texttt{factor()} function in R. Conversely, sometimes string
variables in the data table actually represent ordinal values, one example being
soil drainage being labelled as ``poorly drained'', ``well-drained'', or
``excessively drained''. For analysis, these variables should be manually
reassigned to integer ranking.

\item When one reads in the sample point data table into R with the \texttt{read.dbf()}
function, the first column of the data frame contains the category id value from
the vector sample point layer. This id value obviously is not useful in
multivariate cluster analysis and needs to be omitted. This is easily
accomplished by constructing a new data frame from the original data frame
(function \texttt{data.frame()}), leaving out the category id variable as well as any
other variables that one now wishes to omit. One then assigns the category id
values to the row names with the \texttt{row.names()} function to preserve identification
of the sample points through the cluster analysis.

\item There are two steps to the cluster analysis, creating a dissimilarity matrix
and then running a clustering algorithm on the resulting matrix. In R, the
library that does both functions is \texttt{cluster}.

\item The routine to construct the dissimilarity matrix is called \texttt{daisy()}. It has
options for several dissimilarity metrics but defaults to using Gower's
similarity metric if the inputted data frame has both numeric and categorical
variables. Applying \texttt{daisy()} to the data frame above creates a dissimilarity
matrix where each entry represents the environmental attribute dissimilarity
between a pair of sample points.

\item This dissimilarity matrix is then fed into the \texttt{pam()} clustering routine. A
key parameter that is also input into \texttt{pam()} is the number of clusters to be
generated by the algorithm. Because in this set of studies the clusters were to
be used in a sampling design containing relatively few plots, few clusters were
desired, so the starting point for this parameter was 5 to 7.

\item The R cluster library provides a visualization aid called silhouette
plotting (\cite{rousseeuw1987}) in the function \texttt{silhouette()}. This tool helps identify
clusters which are poorly defined in the sense that many of their members could
easily be placed in another cluster on the basis of dissimilarity distance. By
using this tool while iterating through a broad range of numbers of clusters in
a set of \texttt{pam()} runs, one can identify the choice of the number of clusters that
provides for best cluster discrimination. This choice of the parameter is then
used in selecting the final cluster analysis.

\item Because the PAM algorithm selects particular sample points that are most
representative of each cluster (the ``medoids'') it is easy to inspect the values
at these sample points by extracting rows from the original data frame by the
medoid id number. One can then give the different clusters interpretive labels
e.g. ``Mixed hardwood upland''.

\item Next one brings the result of the cluster analysis back into GRASS. There
are several steps to this process. First, one exports a table from R giving the
id values for the sample points in one column and the cluster analysis
assignments in the second column. (The function \texttt{write.csv()} is useful for this).
Next one needs to assign the cluster labels to the sample points. This is most
easily done by adding a column to the sample point vector layer with the command
\texttt{v.db.addcol} to contain the integer cluster labels. This column can then be 
filled with the cluster labels by running a series of update SQL commands that
assign an integer cluster label based upon the category id value of the sample
points. Shell commands such as awk and sed can be used to generate a file
containing the update SQL commands from the R export table.

\item The original sample points can then be plotted on a map labeled with the
cluster id values. Depending upon the application, it may be desirable to assign
cluster values to the entire surface of the map region, and not just refer to
the labeled sample points. One way to accomplish this is to construct a Voronoi
tessellation using the command \texttt{v.voronoi} and assign each of the demarcated
polygons the cluster value of its enclosed sample point.

\end{enumerate}

\section{Case Study}

As an illustration of this method I present an example from a study of carbon
storage in a mixed agricultural and natural landscape in Northern California
(\cite{williams2011}). Estimation of carbon stocks on the landscape has become
an important concern in light of a desire to mitigate greenhouse gas emissions
through carbon sequestration (e.g. \cite{vankooten2009}). Methodologies for
estimating carbon stocks are only just now being developed, an example being the
California Climate Action Registry Forest Protocol (\cite{climate2009}),
 and these methodologies need to be expanded to a broader suite
of landscapes. This study examined a set of vine tracts and the adjoining
woodland ecosystems from an organic vineyard in Mendocino County (the Bonterra
label of the Fetzer Vineyards). The aim of the study was to estimate carbon
stocks across both the vineyards and wildlands using field sampling and model
building, to compare the aboveground stocks with the soil carbon stocks, and to
develop a methodology that can use these models across a broader landscape for
accurate carbon accounting.

Over the 1149 hectare landscape of the study area, it was only feasible to place
93 vegetation plots and dig 44 soil pits, and to maximize the efficiency of
placing these samples, we followed the PAM cluster analysis methodology outlined
above to classify the landscape for sampling. In this procedure, we used four
different data layers in the GIS data stack. These included vegetation type — a
categorical variable describing the land cover using one of nine types in the
California Wildlife Habitat Relationships system, elevation, slope, and
clear-sky solar radiation summed over the entire year as modeled from the
elevation map. (calculated using the \texttt{r.sun} routine in GRASS). We sampled these
four layers using 1013 points distributed randomly over the Fetzer lands. With
the data matrix from this random sampling, we produced a matrix giving the
dissimilarity between each of the points. This dissimilarity matrix, computed
using Gower's coefficient, was used as the input into the PAM clustering
routine. Based on the distinctiveness of the groupings produced when several
different numbers of clusters were requested, we selected seven as the final
number of clusters. Table \ref{tab:cluster} presents the values of these environmental variables
for the 7 medoid points. These clusters could be characterized as following:two
vineyard types, one on flat ground at low elevations and one upland type; an
annual grassland type at moderate elevation; a mixed hardwood-conifer type at
relatively high elevation and steep slopes; and two mixed hardwood types, one at
relatively high elevation, steep slopes, and lower solar radiation, and the
other at more moderate slopes and elevation. Figure 2 plots the location of the
different clusters on the 5 tracts of the Fetzer lands. 
\begin{table*}\footnotesize
\centering
        \begin{tabular}{|p{1.5cm}|p{1.5cm}|p{2cm}|p{1.5cm}|p{1.5cm}|p{1.5cm}|}
        \hline
        Number of Vegetation Plots & Number of Soil Pits & Habitat Type & 
            Elevation (meters) & Slope (degrees) & Solar Radiation ($w \; hr/m^2/yr$) \\ \hline
        0 & 19 & Vineyard I & 173 & 0.7 & $2.18 \times 10^6$ \\ \hline
        0 & 6 & Vineyard II & 572 & 10.7 & $1.99 \times 10^6$ \\ \hline
        0 & 2 & Annual Grassland & 213 & 5.0 & $2.27 \times 10^6$ \\ \hline
        17 & 1 & Valley Riparian & 152 & 0.4 & $2.18 \times 10^6$ \\ \hline
        39 & 2 & Mixed Hardwoods I & 534 & 21.3 & $1.97 \times 10^6$ \\ \hline
        26 & 13 & Mixed Hardwoods II & 225 & 10.7 & $2.16 \times 10^6$ \\ \hline
        11 & 1 & Mixed Conifer-Hardwood & 407 & 25.8 & $2.06 \times 10^6$ \\ \hline        

        \end{tabular}
        \caption{Values of Environmental Variables for Cluster Medoids in Fetzer Study}\label{tab:cluster}
\end{table*}

Table \ref{tab:cluster} also shows the number of sample points that were assigned to each
cluster. As is evident, the clustering was only used as a rough guide to
distributing cluster points. No vegetation plots were assigned to vineyards or
annual grasslands, since the aim was to focus on woody vegetation types. The
soil pits were divided roughly equally into pits on vineyards and pits on
wildlands, but the wildland pits were predominantly assigned to the Mixed
Hardwoods II class. The vegetation plots and soil pits were assigned to the
cluster types in the field by reference to a map of the cluster points similar
to Figure \ref{fig:clustlabel}. 

\begin{figure}
   \centering
      \includegraphics[width=\linewidth]{hollander/fetzerclusters3} % Point to the graphics file
   \caption{Labeled Cluster Points for Fetzer Study Area} % Provide a caption for the figure
   \label{fig:clustlabel} % change figure_label to a short description that you will use later on
\end{figure}


Rather than estimating by eye what cluster point is nearest a given location, it
is possible to create a surface from the cluster points and assign a cluster to
a location directly. One method for doing this is the Voronoi tessellation
technique outlined. This was not done in the case study, but Figure \ref{fig:voronoi} shows the
result of doing this for the cluster points in the Fetzer study region.

\begin{figure}
   \centering
   \includegraphics[width=\linewidth]{hollander/fetzervoronoi2} % Point to the graphics file
   \caption{Voronoi Tessellation of Fetzer Cluster Points} % Provide a caption for the figure
   \label{fig:voronoi} % change figure_label to a short description that you will use later on
\end{figure}

\section{Conclusion}
This example illustrates how one can use the PAM cluster analysis technique with
a set of GIS layers to guide field sampling for research in landscape ecology.
We have used this technique to sample other agroecosystems, for example in a
study of soil ecology along streams and canals of an agricultural landscape in
Yolo County, California (\citet{culman2010}, \citet*{young-mathews2010}), and
in ongoing work examining soil ecology and gaseous exchange in vineyards in Napa
and San Joaquin counties in California. Many more sophisticated variations on
cluster analysis could be used: for example, there are a number of algorithms
that use spatial constraints directly as part of the clustering process (e.g.
\citet{guo2011}). But as the example shows, in many cases the results of the
clustering are only used as a rough guide for later work, and adopting a simple
algorithm is advantageous. An aim for later work will be to automate this
procedure. This is made possible by using the GRASS Python scripting API to
facilitate the interaction with R in performing the PAM analyses.

\section{Acknowledgements}
I would like to thank Louise E. Jackson, John N. Williams, Steven W. Culman, and
Kerri Steenwerth for the discussions and problem definition that led to the
development of the PAM regionalization algorithm. Support for the Fetzer study
was provided by Brown-Forman Corporation, Fetzer Vineyards, and the Kearney
Foundation of Soil Science with additional support from the USDA-ARS and the
University of California Agricultural Experiment Station.




\bibliographystyle{plainnat}
\bibliography{hollander/Rgrasspam}
		
\address{
Allan D. Hollander\\
Information Center for the Environment\\
University of California, Davis\\
\email{adhollander@ucdavis.edu}
}

% finish the document with this command and we're all done!
\end{article} %% DOCEND
