Due: noon Wed, Jan 21 2015 via GauchoSpace


Scale is an ever present and important topic in ecology. In this lab, you will explore landcover changes in space and time for a city in Santa Barbara County.

Forecasting Change with a Markov Model

Source: Urban, D.L and D.O. Wallin (2001) Introduction to Markov Models. In: Learning Landscape Ecology: A Practical Guide to Concepts and Techniques.

One way to summarize landscape change over time is to tally the instances, on a cell-by-cell basis, of changes from cover type in one period to that in another. For \(m\) cover types, a raw tally matrix then becomes an \(m x m\) matrix. Each element \(n_ij\) in the matrix then corresponds to the number of cells that changed from the cover type \(i\) to cover type \(j\).

A raw tally matrix can then be converted into proportions by dividing the row totals to generate a transition matrix \(P\). The elements, \(p_ij\), of \(P\) are the proportions of cells from the original cover type \(i\) that converted to cover type \(j\). The diagonal elements, \(p_ii\) are the proportions of cells that did not change.

A first-order Markov model (Useher, 1992) assumes that to predict the state of the system at a future time step \(t+1\), one need only know the state of the systme at the current time \(t\). The heart of a Markov model is the transition matrix \(P\), which summarizes the probability that a cell in cover type \(i\) will change to cover type \(j\) during a single time step.

\[ x_{t+1} = x_t P \]

One convenient features is that the next projection can be multiplied to arrive at \(t+2\):

\[ x_{t+2} = x_{t+1} P = x_{t}PP = x_tP^2 \]

so that the state of the system at some future interval \(k\) can be generally be calcualted with:

\[ x_{t+k} = x_t P^k \]

It is worth pointing out that any detected transition from one cover type to another could’ve been interceeded by other transitions at finer temporal transitions.

For more on Markov models, read about Markov Chains.


For your writeup, you’ll just strip away this Rmarkdown document to the necessary code, make requested modifications, answer the questions and render as HTML to turn in via GauchoSpace. Be sure to update your name at the top.

Why R?

Most of this lab could’ve been accomplished using ArcGIS, however when it comes to manipulation of data and modeling, R offers more direct methods. R was originally developed for exploratory data analysis by Bell Labs in the 1970’s, then known as “S” (for statistics) before it evolved into the commercial S-Plus, which prompted a return movement to it’s non-commercial predecessor, hence “R”. R is open source and cross-platform with a wide variety of packages ready to apply (see CRAN Task Views, eg for Spatial).

Open Rmarkdown in RStudio

Normally R launches just a command window. You’ll instead use the friendlier RStudio interface.

  1. Extract lab2_scale.z to your H:215. Copy the following file:


    into your home directory for the class H:\esm215, right-click and unzip there. The remainder of this lab will assume your working directory is H:\esm215\lab2_scale. If you use a different working directory, just set that at the beginning of the script below.

  2. Open lab2_scale.Rmd. The default editor for Rmarkdown files (explained later) is Tinn-R. We’ll use the more capable RStudio editor instead, which you can invoke from Windows Explorer by right-clicking on the lab2_scale.Rmd file -> Open with -> RStudio. May be slow to launch.

  3. Update raster package. The raster package on the lab machines needs to get updated to the latest version to avoid producing an error message . Go to the Packages pane and click the Update button. This might take a couple minutes to get the full listing. Scroll way down to raster, and tick just the box for updating the raster package, and click Install Updates button. (Do not update all packages because this will take a long time and hang up your RStudio session.) If all goes well, you should be able to type “raster”" in the Packages search bar (after a Refresh) and see that it’s registered as Version 2.3-12 like so:

    • select with mouse, Ctrl+Enter to run

You’ll notice several panes:

For help on function, type help(raster) or shortcut ?raster into the Console or simply place cursor on function in a script of Editor pane and hit the F1 key. You can also search the help for functions across available packages like so: help_search('zonal') or shortcut ??zonal.

What is R Markdown (*.Rmd)?

R Markdown is a file format allowing you to weave easy formatting of text (markdown) with code (R). For more, see:

R Basics

Here are a few basics:

  • Commments are preceeded with # and show up as light green in the RStudio Editor and light grey in the rendered Rmarkdown document.

  • White aspace does not matter, but indentation is a good idea for readability.

  • Functions have arguments: function(argument1, argument2, ...).

  • Variable assignment(= or <-) to 'string' or '"string"' vs number 1 vs boolean TRUE/FALSE or T,F.


Read Data

# load necessary libraries having useful functions, suppressing startup messages and warnings
  library(reshape) # reshaping functions: rename
  library(dplyr)   # dataframe functions: select
  library(rgdal)   # read/write raster/vector data
  library(raster)  # raster functions
  library(rgeos)   # vector functions
  library(knitr)   # knitting functions: kable

# Set your working directory

# read in landcover rasters, relative to working directory
nlcd_2001 = raster('raster/nlcd_2001.tif')
nlcd_2006 = raster('raster/nlcd_2006.tif')
nlcd_2011 = raster('raster/nlcd_2011.tif')

# read in county and city boundary vectors
# note arguments to readOGR function for ESRI Shapefile: 
#  data source name (dsn) = directory
#  layer = shapefile filename without shp extension
cities = readOGR('vector', 'city_boundary_no_ocean')
## OGR data source with driver: ESRI Shapefile 
## Source: "vector", layer: "city_boundary_no_ocean"
## with 8 features and 5 fields
## Feature type: wkbPolygon with 2 dimensions

Plot and Transform

Try plotting the cities on top of landcover.

# try plotting landcover and cities on top
plot(cities, border='black', lwd=2, add=T)