R/02_data-processing-competition-coef.R

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

######################################################################
#
#  competition coefficients and stem-maps based on xy coordinates
#
#  Rob Smith, phytomosaic@gmail.com, 01 Apr 2020
#
##  CC-BY-SA 4.0 License (Creative Commons Attribution-ShareAlike 4.0)


### load data
rm(list=ls())
require(data.table)
load('./data/pnw_tree.rda', verbose=T)
d <- pnw_tree  ;  rm(pnw_tree)

### tre_cn is indeed a unique tree identifier
anyDuplicated(d$tre_cn) # expect 0

### FIA plot layout diagram
`plot_stems` <- function(d, pltcn) {
        if(!any(d$plt_cn == pltcn)) stop('`plt_cn` not in data')
        cx <- cy <- 120        # origin's xy coords at 120-ft
        r  <- 120              # radius for placing subplots
        a  <- c(0,120,240)+90  # vector of angles (degrees)
        `rad` <- function(deg) {deg*pi/180} # function for radians
        x1 <- c(cx, cx + r * cos(rad(a)))   # x for subplots 2,3,4
        y1 <- c(cy, cy + r * sin(rad(a)))   # y for subplots 2,3,4
        z <- d[plt_cn == pltcn,]
        z <- z[order(z$subp, z$dist),]
        if(all(is.na(z$x)) | all(is.na(z$y)))
                stop('no coords for this plot')
        # PNW 58.9-ft macroplots
        symbols(x=x1,y=y1,circ=rep(58.9,length(x1)),inches=F,xaxt='n',
                yaxt='n', bty='n', ann=F, asp=1, ylim=c(0,300))
        # 6-ft buffer
        symbols(x=x1,y=y1,circ=rep(52.9,length(x1)),inches=F,fg=2,add=T)
        # FIA 24-ft subplots
        symbols(x=x1,y=y1,circ=rep(24, length(x1)),inches=F,add=T)
        # 6-ft buffer
        symbols(x=x1,y=y1,circ=rep(18,length(x1)),inches=F,fg=2,add=T)
        for(i in 1:4) { # translate stem positions in each subplot
                # not within 6-ft of edge
                points(x1[i] + z$x[z$subp==i & !z$is_edge],
                       y1[i] + z$y[z$subp==i & !z$is_edge],
                       pch=16, cex=0.5, col='#000000')
                # is within 6-ft of edge
                points(x1[i] + z$x[z$subp==i & z$is_edge],
                       y1[i] + z$y[z$subp==i & z$is_edge],
                       pch=16, cex=0.5, col='#FF0000')
        }
}
### view a single plot
plot_stems(d, pltcn=unique(d$plt_cn)[8688])
### view a single plot that should contain seedlings
plot_stems(d, pltcn=unique(d[is_seed==TRUE,plt_cn])[8])
### view several
par(mfrow=c(2,2), pty='s')
for (i in 1:4) {
        plot_stems(d, pltcn=unique(d$plt_cn)[2020+i])
}

### Jenny needs a file where:
#     fid    = focal tree id
#     fspp   = focal tree species
#     nspp   = neighbor tree species
#     nsize  = neighbor tree size (inches) from first measurement
#     year   = year of first measurement
#     subpid = subplot id (following convention 'plot_subplot')
### function to make such a file:
`f` <- function(d, dist_cutoff=18.0, size_cutoff=5.0) {
        # omit if xy are missing
        d <- d[!(is.na(x) & is.na(y)),]
        # omit if BOTH diameters missing
        d <- d[!(is.na(dia_y1) & is.na(dia_y2)),]
        # establish ONE DBH for each (regardless survivor/recruit)
        tmp  <- cbind(dia_y1=d$dia_y1, dia_y2=d$dia_y2)
        isna <- is.na(tmp[,1])
        tmp[isna,1] <- tmp[isna,2]
        d$DBH       <- tmp[,1] # new column DBH
        rm(tmp)
        if (is.numeric(dist_cutoff)) { # keep only trees in buffer
                d <- d[d$dist > dist_cutoff,]
        }
        if (is.numeric(size_cutoff)) { # keep only non-seedlings
                d <- d[d$DBH > size_cutoff,]
        } # whole dataset has 170914 individuals remaining....
        ### begin loop for each subplot ! ! ! TIMEWARN ! ! !
        out <- data.frame(matrix(NA, nrow=0, ncol=6)) # collector
        nsubpid <- length(unique(d$subpid))
        iter <- 0
        pb   <- txtProgressBar(min = 0, max = nsubpid, style = 3)
        for (k in unique(d$subpid)) { # loop for each subplot
                # cat((iter <- iter + 1) / nsubpid, '...')
                z   <- d[d$subpid==k,]       # trees in single subplot
                if (NROW(z)==1) { # if no neighb, then competition = NA
                        m <- data.frame(
                                z[,c('tre_cn','spp','spp',
                                     'DBH','measyr_1','subpid')])
                        names(m) <- paste0('X',1:6)
                        m[,c(3:4)] <- NA
                } else {
                        idx <- t(combn(1:NROW(z),2)) # index row-pairs
                        m <- data.frame(matrix(NA,nrow=NROW(z),ncol=6))
                        for (i in 1:NROW(idx)) { # for each pair
                                m[i,] <- c(
                                        # focal
                                        z[idx[i,1], c('tre_cn','spp')],
                                        # neighbor
                                        z[idx[i,2], c('spp','DBH',
                                                      'measyr_1',
                                                      'subpid')])
                        }
                }
                out <- rbind(out, m)
                setTxtProgressBar(pb, (iter <- iter + 1))
        }
        close(pb)
        names(out) <- c('fid','fspp','nspp','nsize','year','subpid')
        return(out)
}
a <- f(d, dist_cutoff=18.0, size_cutoff=5.0)
# save(a, file='C:/Users/Rob/Desktop/a.rda')

### TODO: calculate competition coefficients from `a`

# ### TODO: Jenny to calculate neighborhood competition coefficients
# `calc_comp` <- function(d, dist_cutoff=18.0) {
#         d <- d[d$dist > dist_cutoff,] # keep only trees inside buffer
#         for (k in unique(d$subpid)) { # loop for each subplot
#                 z <- d[d$subpid==k,]
#                 ### TODO: calculate competition coefficients in 'z'
#         }
#         return('TODO!')
# }

####    END    ####