##=====================================================
##=========== FOREST INVENTORY SCRIPT =================
##=====================================================
##====== NORTHWOODS NATURAL RESOURCES, LLC ============
##======    WWW.NORTHWOODSRESOURCES.COM    ============
##======     PREPARED BY: LEE MUELLER	   ============
##=====================================================
library('reshape')

## This script is designed to take options for various
## forest inventory systems to analyze and generate
## inventory reports for various forestry purposes.
## configure the file accordingly and pay attention to
## report generated locations.

## ===================================================== CONFIGURE
## Set directory in which the file is housed.
setwd("/home/lee/Documents/Northwoods/Inventories/")
## Your data will need to have the following variables:
name <- "West Fork Tree Farm" #Project Name
proj <- "McKarns" #Unique project ID
date <- "" #Inventory Date
inv <- read.csv("mckarn_inventory.csv") #File Name

## ===================================================== CODE
attach(inv) #format data
names(inv)

uniquestands <- unique(stand) #return unique stands
numstands <- length(uniquestands) #get number of stand

for ( i in 1:numstands) #open loop 1-number of stands
{

	## Separate trees by stand
	standsub <- subset(inv, stand==	i) #set subsample based on plot number. Note loop variable.
	colnames(standsub) <- names(inv)

	#Standard mathematics for samples. Based on 10 BAF point sample.
	standsub$numplots <- apply(standsub,1,function(row) length(unique(standsub$plot))) #add numplot column for number of plots per stand for later adjustment
	standsub$tpa <- (1/( (((standsub$dbh * 2.75)^2) * 3.14159)/43560))/standsub$numplots #add column for TPA of each sample

	#Need to begin creating a table which can be graphed. We need one which is TPA by species in each size class. Standard stand table.
	#Create appropriate dataframe and organize to be pivoted correctly.
	graph <- data.frame(standsub$species, standsub$dbh, standsub$tpa)
	names(graph) <- c("species", "dbh", "tpa")
	attach(graph)
	names(graph)

	#We want ACTUAL species names to show up in the graph. Let's do some replacing now.
	codes <- read.csv("speciescodes.csv") #Load the species codes file
	attach(codes)
	graph$species <- codes$text[match(graph$species, codes$code)] #match and replace

	# Make pivot table. TPA by DBH and species. TPA is the variable of interest being summated.
	pivot <- cast(graph, dbh ~ species, sum)

	#set dynamic filename
	filename <- paste( "plots/", proj, "-stand-", i, ".jpg", sep="") #if it doesn't work, this might be a problem.

	#open plot save function
	jpeg(filename)

	#Make a stacked barplot. X axis = DBH Y axis = TPA and stacks = species.
	par(xpd=T, mar=par()$mar+c(0,0,0,4))
	barplot(t(pivot), ylab="Trees Per Acre", xlab="Diameter at Breast Height (inches)", space=0.1, col=rainbow(12))

	#For some reason, "DBH" variable shows up as a label in the plot. This bit removes it.
	labels <- names(pivot)
	labels <- labels[2:length(labels)]

	#set legend with changed labels
	legend("topright", labels, fill=rainbow(12))

	# Restore default clipping rect
	par(mar=c(5, 4, 4, 2) + 0.1)

	#close save
	dev.off()
} #close loop