R/04_plot-species-distributions.R

In phytomosaic/pnw: Tree functional strategies in the US Pacific Northwest

######################################################################
#
#  Plot species distributions on PNW map
#
#  Rob Smith, phytomosaic@gmail.com, 30 Mar 2020
#
##  CC-BY-SA 4.0 License (Creative Commons Attribution-ShareAlike 4.0)


### preamble
rm(list=ls())
# devtools::install_github('phytomosaic/ecole')
require(ecole)
require(viridis)
require(rgdal)
cc <- viridis::viridis(99)
getwd()
load('./data/pnw_tree.rda', verbose=T)
d <- pnw_tree  ;  rm(pnw_tree)

### first, plot ecosubregions
xtabs(~ewh+ew, data=d)
# png('./fig/fig_00_eastwest.png',
#     wid=4.5, hei=4.5, uni='in', bg='transparent', res=700)
set_par(1)
plot(d$lon, d$lat, col=colvec(as.factor(d$ewh)),
     pch=16, cex=0.4, asp=1.4, ylab='Latitude', xlab='Longitude')
legend('bottomright', legend=levels(as.factor(d$ewh)),
       fill=colvec(1:4, alpha=1), bty='n')
# dev.off()

### second, plot species distributions
# aggregate: periodic ann increment for each species
frq  <- aggregate(d$spp, list(d$spp), length)
(frq <- frq[rev(order(frq$x)),])
nm   <- as.character(frq[,1])
# rescale PAI, just for plotting
d$ann_incr <- c(d$dia_y2 - d$dia_y1) ^ 0.7
# fix one species that had NO survivors, just for plotting
d$ann_incr[d$spp=='malus_fusca'] <- runif(sum(d$spp=='malus_fusca'))
# aggregate mean per plot
p <- aggregate(list(ann_incr = d$ann_incr),
               list(plt_cn=d$plt_cn, spp=d$spp),
               mean, na.rm=T, drop=F)
p <- cbind(p,d[match(p$plt_cn, d$plt_cn),
               c('lat','lon','elev')])
# shapefiles
trg_prj <- CRS('+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-95
           +x_0=0 +y_0=0 +datum=NAD83 +units=km +no_defs')
# read PCT shapefiles
w <- rgdal::readOGR(dsn='./data_raw/gis/PacificCrestTrail',
                    'PacificCrestTrail')
w <- spTransform(w, trg_prj)    # reproject Lambert Conformal Conic
# read state shapefiles
s <- tigris::states(cb=T, resolution='20m') # tigris states map
s <- s[s@data$STUSPS %in% c('CA','OR','WA'),]
s <- spTransform(s, trg_prj)    # reproject Lambert Conformal Conic
# reproject points
`make_xy` <- function(xy, crs, ...){
    xy <- as.data.frame(xy)
    xy <- xy[!(is.na(xy$lon) | is.na(xy$lat)),] # rm NAs
    coordinates(xy) <- ~lon+lat
    proj4string(xy) <- CRS('+init=epsg:4269')  # NAD83 dec deg
    return(spTransform(xy, crs))
}
p <- make_xy(p, crs=trg_prj)   # reproject Lambert Conformal Conic
# function to plot
`map_points` <- function (p, state=s, pct=w, main='', ...) {
    plot(state, lwd=.5, lty=1, col='grey70', axes=F, xlab='', ylab='',
         xlim=c(-2471.7435, -1626.256), ylim=c(300, 1468.024))
    plot(p, pch=16, col=ecole::colvec(p$ann_incr,alpha=0.95),
         cex=0.6, add=T)
    plot(pct, lwd=1, col='gold', add=T) # add PCT
    title(main=main, line=-2, cex=0.9)
}
### loop over all 57 species
n <- 12 # number of panels per figure
for (k in 1:5) {
    # k <- 1
    png(paste0('C:/Users/Rob/Desktop/fig_0',k,'_PAI.png'),
        wid=9, hei=5, unit='in', bg='transparent', res=700)
    layout(matrix(1:n, nrow=2, ncol=6, byrow=T))
    par(oma=c(0,0,0,0), mar=c(.1,0,.1,.1))
    sk <- ((k*n)-(n-1)):(n*k)
    sk <- sk[sk <= 57] # not to exceed 57 species
    for (i in sk) {
        map_points(p[which(p$spp == nm[i]),] , main=nm[i])
    }
    dev.off()
}
d$ann_incr <- NULL  # remove temporary values
####    END    ####