This document is the result of a collaborative effort between Donal O’Leary, Katherine Rosa, and Trevor Bloom in support of our graduate educations here at Western Washington University.

Introduction

Trevor Bloom is currently pursuing an M.S. in Biology studying the effects of wildfire and climate change on the distribution of the flower species Saxifraga austromontana. More about his research may be found on his website here

For this project Trevor has modeled a climactic niche that these plants occupy. Katie and I then took those models and extracted the elevation and latitude for each pixel in the DEM that has its centroid within one of these suitable habitat polygons.

Here this image shows the current SDM (Species Distribution Model) and the future SDM for Saxifraga austromontana. Using these polygons we extracted the elevation and latitude values from the DEM file provided. Map by Trevor Bloom

R Code

First, let’s do a little housekeeping.

## This code produced by Donal O'Leary and Katherine Rosa
library(sp)
library(spatstat)
## 
## spatstat 1.41-1       (nickname: 'Ides of March') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.41-1 is out of date by more than 3 months; we recommend upgrading to the latest version.
library(ncf)
library(spdep)
## Loading required package: Matrix
#install.packages("sp")
#install.packages("spatstat")
#install.packages("ncf")
#install.packages("spdep")

Present Suitable Habitat

# Import csvs of lat and dem information
present.rawlat=read.csv("/Users/bgatts/Desktop/WWU/Katie_R_Help/extract_present_lat.txt", stringsAsFactors=F)
present.rawdem=read.csv("/Users/bgatts/Desktop/WWU/Katie_R_Help/extract_present_dem_corrected.csv", stringsAsFactors=F)

dem=present.rawdem$VALUE
dem=data.frame(dem)
dem$lat=present.rawlat$VALUE

For reasons that may be apparent later on, we need to scale these values from zero to one for both axes so that it is easy to plot them in a square window. To do so we create a new data.frame.

# Scale dem and lat
df=data.frame(dem$dem/4500)
df$lat=dem$lat/70
names(df)=c("dem","lat")

Here we define the “window” required for the ppp object.

# Define Window
window=owin(c(0,1),c(0,1))

# Create ppp object
ppp1=as.ppp(df, W=window)
## Warning in ppp(X[, 1], X[, 2], window = win, marks = marx, check = check):
## data contain duplicated points
# plot density
plot(density(ppp1))

Next we adjust the image to remove the lower 30 degrees of latitude which are not represented in this dataset, giving us a more clear image.

#re-scale latitude to get rid of the bottom half which is unused
df$lat=(df$lat-.42857)/.57143

# Create ppp object
ppp2=as.ppp(df, W=window)
## Warning in ppp(X[, 1], X[, 2], window = win, marks = marx, check = check):
## data contain duplicated points
# plot density
plot(density(ppp2))

Let’s briefly compare these plots:

#Set up the two plots per image
par(mfrow=c(1,2))

plot(density(ppp1))
plot(density(ppp2))

Future Suitable Habitat

#First, import the CSV.
#Remember, the DEM needed to be corrected to read in properly.
future.rawlat=read.csv("/Users/bgatts/Desktop/WWU/Katie_R_Help/extract_future_lat.txt", stringsAsFactors=F)
future.rawdem=read.csv("/Users/bgatts/Desktop/WWU/Katie_R_Help/extract_future_dem_corrected.csv", stringsAsFactors=F)

#Create new Object
dem.fut=future.rawdem$VALUE
summary(dem.fut)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      10    1419    1679    1841    2099    4244
#Turn it into a data.frame
dem.fut=data.frame(dem.fut)
head(dem.fut)
##   dem.fut
## 1    1340
## 2    1203
## 3    1181
## 4    1019
## 5    1200
## 6    1186
#Add latitude
dem.fut$lat=future.rawlat$VALUE
head(dem.fut)
##   dem.fut      lat
## 1    1340 65.22613
## 2    1203 65.18728
## 3    1181 65.17433
## 4    1019 65.13548
## 5    1200 65.16138
## 6    1186 65.04482
# 1x1 window, requires data to be scaled.
window2=owin(c(0,1),c(0,1))

# Scale variables
# Will want to add some code here to properly scale the variables, might want to also adjust the window
# For these max values, enter in the actual upper limit (or like 4,500 meters for DEM and 70 deg north for lat)
df.fut=data.frame(dem.fut$dem/4500)
df.fut$lat=dem.fut$lat/70
names(df.fut)=c("dem","lat")


# Create ppp object
ppp.fut1=as.ppp(df.fut, W=window2)
## Warning in ppp(X[, 1], X[, 2], window = win, marks = marx, check = check):
## data contain duplicated points
# plot density
plot(density(ppp.fut1))

#re-scale latitude to get rid of the bottom half which is unused
df.fut$lat=(df.fut$lat-.42857)/.57143

# Create ppp object
ppp.fut2=as.ppp(df.fut, W=window2)
## Warning in ppp(X[, 1], X[, 2], window = win, marks = marx, check = check):
## data contain duplicated points
# plot density
plot(density(ppp.fut2))

#Set up the two plots per image
par(mfrow=c(1,2))

plot(density(ppp.fut1))
plot(density(ppp.fut2))

Final Plots

Let’s clean up these plots some - nice work Katie

library(colorRamps)
par(mfrow=c(1,2))
plot(density(ppp2), col=blue2red(10), xlab="Elevation",ylab="Latitude", zlim=c(0,1900000))
title(main="Present Distribution Density", xlab= "Elevation (m)", ylab= "Latitude (degrees)")
axis(side=1, at=seq(0,1, .11111), labels=seq(0,4500,500))
axis(side=2, at=seq(0,1, .25), labels=seq(30,70,10))
plot(density(ppp.fut2), col=blue2red(10), xlab="Elevation",ylab="Latitude", zlim=c(0,1900000))
title(main="2080 Distribution Density", xlab= "Elevation (m)", ylab= "Latitude (degrees)")
axis(side=1, at=seq(0,1, .11111), labels=seq(0,4500,500))
axis(side=2, at=seq(0,1, .25), labels=seq(30,70,10))

Elevation and Latitude Center of Distribution

Here we calculate the center of the density models shown above. The code is a bit opaque, but we essentially called for the highest value in the density plot and then back-tracked its row and column from its position within the 128X128 matrix. There is almost certainly an easier way to do this, but this technique appears to work.

# center of density x & y of max density
present.density=density(ppp1)
future.density=density(ppp.fut1)
#which.max returns the position in the number array $v which is basically the pixel value of the density image
# this position must then be interpreted in terms of rows/columns
present.maxval=which.max(present.density$v)
#4934/128
#[1] 38.54688
#> 128*38
#[1] 4864
#> 4934-4864
#[1] 70
future.maxval=which.max(future.density$v)
#5967/128
#[1] 46.61719
#> 128*46
#[1] 5888
#> 5967-5888
#[1] 79

## Basically we found that the mad density for the present is at (column 38, row 70 [counted from the bottom])
## Basically we found that the mad density for the future is at (column 46, row 79 [counted from the bottom])
# we then convert those to lat and elevation
#present:
present.max.density.lat=40*70/128+30
present.max.density.elev=4500*38/128
## Present max density lat=51.875 Elev=1336
future.max.density.lat=40*79/128+30
future.max.density.elev=4500*46/128
## Future max density lat=54.6875 Elev=1617

Histograms!

Here is another look at the data. These histograms show the current and future distributions captured by these SDMs (species distribution models).

## Histograms
## Two per plot
#par(mfrow=c(1,2))
## Elevation
#hist(dem$dem, main="Current Pixel Distribution - Elevation", xlab="Elevation (meters)")
#hist(dem.fut$dem, main="Future Pixel Distribution - Elevation", xlab="Elevation (meters)")

hist(dem$dem, col=rgb(1,0,0,0.5), main="Present vs Future Elevation", xlab="Elevation(meters)", ylab="Number of Pixels")
hist(dem.fut$dem, col=rgb(0,0,1,0.5), add=T)
legend(3000, 40000, c("Present","Future"), fill=c("red","blue"))
abline(v=mean(dem$dem), col="red", lwd=3)
abline(v=mean(dem.fut$dem), col="blue", lwd=3)

## Latitude

hist(dem$lat, col=rgb(1,0,0,0.5), main="Present vs Future Latitude", xlab="Latitude(Degrees)", ylab="Number of Pixels")
hist(dem.fut$lat, col=rgb(0,0,1,0.5), add=T)
legend(32, 53000, c("Present","Future"), fill=c("red","blue"))
abline(v=mean(dem$lat), col="red", lwd=3)
abline(v=mean(dem.fut$lat), col="blue", lwd=3)

The lines plotted here show the mean value for the present and future distributions.