demo/pub_code/R.Figures/figure.trend.empir.vs.mod.R

In USGS-R/mda.lakes: Tools for combining models, data, and processing for lakes

## Figure to compare empirical vs modeled trends
library(plyr)
library(nlme)

better.JAS = function(df){
  if(nrow(df) > 4){
    return(mean(df$V1))
  }else{
    return(NA)
  }
}

better.trend = function(df){
  if(nrow(df) > 10 & abs(nrow(df) - (max(df$year) - min(df$year))) < 3){
    return(data.frame(slope=lm(df$V1 ~ df$year)$coeff[2], syear=min(df$year)))
  }else{
    return(data.frame(slope=NA, syear=NA))
  }
}

wtr = read.table('../supporting files/wtemp.obs.tsv', sep='\t', header=TRUE)

wtr.near.surf = ddply(wtr, c("WBIC", "DATETIME"), function(df) mean(df$WTEMP[df$DEPTH<=5], na.rm=TRUE))

mons = as.POSIXlt(wtr.near.surf$DATETIME)$mon+1

wtr.near.surf = wtr.near.surf[mons >=7 & mons <= 9,]
wtr.near.surf$year = as.POSIXlt(wtr.near.surf$DATETIME)$year+1900

JAS.mean = ddply(wtr.near.surf, c("WBIC", "year"), better.JAS)
JAS.mean = JAS.mean[!is.na(JAS.mean$V1),]

slopes = ddply(JAS.mean, c("WBIC"), better.trend)
slopes = slopes[!is.na(slopes$slope), ]

data = read.table('../Output/omg.huge.output.tsv', sep='\t', header=TRUE)


lter.wbic = read.table('../supporting files/lter.south.lakes.tsv', header=TRUE, sep='\t')
lter.wbic$region = "South"
epi.temps = read.table('../Output/KdScenarios/stable.metrics.out.tsv', sep='\t', header=TRUE)

#lter.mod = data[data$lakeid %in% lter.wbic$WBIC, ]
lter.mod = epi.temps[epi.temps$lakeid %in% lter.wbic$WBIC,]
lter.mod = lter.mod[lter.mod$year >= 1995 & lter.mod$year <= 2011, ]

lter.data = JAS.mean[JAS.mean$WBIC %in% lter.wbic$WBIC,]
lter.data = lter.data[lter.data$year >= 1995 & lter.data$year <= 2011, ]

lter.mod.slopes = ddply(lter.mod, c("lakeid"), function(df) lm(df$mean_epi_temp ~ df$year)$coeff[2])
lter.mod.slopes.sens = ddply(lter.mod, c("lakeid"), function(df) zyp.sen(mean_epi_temp ~ year, df)$coeff[2])

lter.slopes = ddply(lter.data, c("WBIC"), function(df) lm(df$V1 ~ df$year)$coeff[2])
lter.slopes.sens = ddply(lter.data, c("WBIC"), function(df) zyp.sen(V1 ~ year, df)$coeff[2])


names(lter.mod.slopes.sens) = c("WBIC", 'slope.mod')
names(lter.slopes.sens) = c('WBIC', 'slope.empir')


lter.empir.mod = join(lter.mod.slopes.sens, lter.slopes.sens)
lter.slopes.south = join(lter.empir.mod, lter.wbic)

lter.wbic = read.table('../supporting files/lter.north.lakes.tsv', header=TRUE, sep='\t')
lter.wbic$region = "North"
epi.temps = read.table('../Output/KdScenarios/stable.metrics.out.tsv', sep='\t', header=TRUE)

#lter.mod = data[data$lakeid %in% lter.wbic$WBIC, ]
lter.mod = epi.temps[epi.temps$lakeid %in% lter.wbic$WBIC,]
lter.mod = lter.mod[lter.mod$year >= 1985 & lter.mod$year <= 2011, ]

lter.data = JAS.mean[JAS.mean$WBIC %in% lter.wbic$WBIC,]
lter.data = lter.data[lter.data$year >= 1985 & lter.data$year <= 2011, ]

lter.mod.slopes = ddply(lter.mod, c("lakeid"), function(df) lm(df$mean_epi_temp ~ df$year)$coeff[2])
lter.mod.slopes.sens = ddply(lter.mod, c("lakeid"), function(df) zyp.sen(mean_epi_temp ~ year, df)$coeff[2])

lter.slopes = ddply(lter.data, c("WBIC"), function(df) lm(df$V1 ~ df$year)$coeff[2])
lter.slopes.sens = ddply(lter.data, c("WBIC"), function(df) zyp.sen(V1 ~ year, df)$coeff[2])


names(lter.mod.slopes.sens) = c("WBIC", 'slope.mod')
names(lter.slopes.sens) = c('WBIC', 'slope.empir')


lter.empir.mod = join(lter.mod.slopes.sens, lter.slopes.sens)
lter.slopes.north = join(lter.empir.mod, lter.wbic)

lter.slopes = rbind(lter.slopes.north, lter.slopes.south)


tiff('NTL.slopes.tiff', width=1200, height=1200, res=300, compression='lzw')
boxplot(slope.empir ~ region, data=lter.slopes, lwd=2, ylab='JAS Trend', xlab='NTL Region')
dev.off()

tiff('mod.v.empir.trends.tiff', width=1200, height=1200, res=300, compression='lzw')
plot(lter.slopes$slope.empir, lter.slopes$slope.mod, 
     xlab="Empirical trends", ylab="Modeled trends", pch=16, xlim=c(-0.03, 0.18), ylim=c(-0.03, 0.18))
lines(c(-1,1), c(-1,1), lwd=2)
dev.off()