# --  abundance weighted plot trait means
# !!! make sure all traits are saved as numbers instead of factors !!!
setwd("C:/Documents and Settings/kreitler/My Documents/FIA")
data1=read.csv("plot_spp_abun3-28.csv", header=TRUE)  #plot spcd #individuals
traits=read.csv("species.traits.csv",header=TRUE)    #check column names -- they change, which messes up code
data2=merge(data1,traits,all.x=TRUE)
trees=read.csv("trees.i.csv",skip=4,header=TRUE)
#write.csv(data2,"data2.csv")
env=read.csv("plots_env.csv", header=TRUE)
plotList=unique(data1$PLOT)

head(data1)
head(traits)
head(data2)

# make the matrix before hand so it doesn't bomb
data2$relativeAbund<-NA
plotseedmass<-c()
plotwoodspecgrav<-c()
plottreetype<-c()
fireresist<-c()
resprout<-c()
sla<-c()
maxheight<-c()
#firetolerance<-c()
#serotiny<-c()

# -- plot level trait weighted mean loop ----
for(i in 1:length(plotList)   )
{
      temp<- subset(data2, data2$PLOT==plotList[i])      # loop through each plot [i]
      indivs<- sum(temp$CountOfCN)                       # number of individuals at each plot
      data_rows<-as.integer(row.names(temp))             # ?
      relAbund <-temp$CountOfCN/indivs                   # relative abundance
      data2$relativeAbund[data_rows] <-  relAbund        # write relative abundance to data2
      plotseedmass[i]<- weighted.mean(temp$seedmass_mg, temp$relativeAbund,na.rm = TRUE) #make sure column names are right
      plotwoodspecgrav[i]<- weighted.mean(temp$woodspecgrav, temp$relativeAbund,na.rm=TRUE)
      resprout[i]<- weighted.mean(temp$resprout, temp$relativeAbund, na.rm=TRUE)
      sla[i]<- weighted.mean(temp$sla,temp$relativeAbund,na.rm=TRUE)
      plottreetype[i]<- weighted.mean(temp$treetype, temp$relativeAbund,na.rm=TRUE)
      maxheight[i]<- weighted.mean(temp$maxheight,temp$relativeAbund,na.rm=TRUE)

          
}
result<-data.frame(plotList,plotseedmass,plotwoodspecgrav,resprout,sla,plottreetype,maxheight)
write.csv(result,"plottraitmean_5-15.csv")
head(result)

# correlation matrix -- PTMs, env variables
envplot=merge(env,result, by.x = "PLOT", by.y = "plotList")
head(envplot)
#write.csv(envplot,"envplot.csv")  
envplot2=data.frame(envplot$LAT,envplot$ELEV,envplot$CA_PPT,envplot$CA_TMAX,envplot$CA_TMIN,envplot$AvgOfVOLGCS,envplot$Carbon,envplot$plotseemass,envplot$plotwoodspecgrav,envplot$resprout,envplot$sla,envplot$plottreetype,envplot$maxheight)
head(envplot2)
plot(envplot2)

# CART models
library(rpart)
attach(envplot)

plot(envplot$plotwoodspecgrav~envplot$ELEV)
fit1=rpart(plotwoodspecgrav~ELEV+CA_PPT+CA_TMAX+CA_TMIN+LAT, data=envplot)
plot(fit1)
text(fit1)
lm1=lm(plotwoodspecgrav~ELEV+CA_PPT+CA_TMAX+CA_TMIN+LAT, data=envplot)
summary(lm1)

fit2=rpart(log(plotseedmass)~ELEV+CA_PPT+CA_TMAX+CA_TMIN+LAT, data=envplot)
plot(fit2)
text(fit2)
lm2=lm(log(envplot$plotseedmass)~envplot$CA_TMAX)

hist(log(AvgOfVOLCFGRS))
fit3=rpart(log(AvgOfVOLCFGRS)~ELEV+LAT+CA_PPT+CA_TMAX+CA_TMIN)
plot(fit3)
text(fit3)
lm3=lm(log(AvgOfVOLCFGRS)~CA_PPT)
plot(log(AvgOfVOLCFGRS)~CA_PPT)
abline(lm3)

detach(envplot)