demo/pub_code/R.Figures/grl.figure3.R

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

###

library(plotrix)
library(data.table)
library(plyr)
source('sens.confint.mod.R')

## Just a try
all.med.area = c()
all.med.slope = c()
all.index = c()
source('grl.resample.calcs.R')

for(l in 1:1000){
	all.slopes = data.table(grl.resample.calcs())
	#all.slopes = fread('downsampled.slopes.csv')
	
	#meandt = ddply(all.slopes,'wbic',function(df)mean(df$dt))
	#meandt = meandt#[meandt$V1 > 5,]
	all.slopes = all.slopes[wbic%in%sample(unique(all.slopes$wbic),70), ]
	#nrow(meandt)
	
	#all.slopes = all.slopes[wbic%in%meandt$wbic, ]
	
	wbic.size = unique(all.slopes[,list(wbic,area)])
	wbic.size = wbic.size[order(area),]
	
	n = nrow(wbic.size)-30
	med.area = rep(NA,n)
	med.slope = rep(NA,n)
	
	for(i in 1:n){
		
		bot = wbic.size[i:(i+29),]
		
		med.slope[i] = median(all.slopes[wbic%in%bot$wbic,]$slopes)
		med.area[i] = median(bot$area)
		
	}
	all.med.area = c(all.med.area, med.area)
	all.med.slope = c(all.med.slope, med.slope)
	all.index = c(all.index, med.slope = rep(l,n))
	cat(l,'\n')
}

tiff('figure.3.2.tiff', width=2000, height=1500, res=300, compression='lzw')

plot(all.med.area,all.med.slope, log='x', xaxt='n', pch=20, bg=rgb(0.5,0.5,0.5, 0.025), col=rgb(0.5,0.5,0.5, 0.05),
		 ylab=expression(Median~Trend~(degree*C~yr^-1)), xlab=expression(Median~Lake~Area~(km^2)), cex=0.5)

tmp = loess(all.med.slope~all.med.area)
lines(sort(all.med.area), predict(tmp, sort(all.med.area)), lwd=2, col='black')
axis(1, at=c(5e5, 1e6, 2e6, 5e6, 1e7), 
		 labels=c("0.5", "1", "2", "5", "10"))
		 #labels=c(expression(''*5*x*10^5),expression(''*1*x*10^6),expression(''*2*x*10^6),expression(''*5*x*10^6), expression(''*1*x*10^7)))

dev.off()