## unholtify(); Clarke Cooper, 2019.
## Adapted from Z.Lin's answer to me at
## https://stackoverflow.com/questions/54695468/
# r-translating-recursively-divided-areas-into-xy-coordinates-for-ggplot2-chlor

library(dplyr)
library(ggplot2)

## The Holt Research Forest study area is divided into 100 x 100m 
## blocks, which are divided into 50 x 50m quadrats and 25 x 25 subquadrats.
## unholtify() converts HRF block-quad-subquad locations into [x,y]
## coordinates for center of each square so ggplot can make maps of
## the data -- or any other reason you might want a more rastery arrangement.  
## Not entirely bombproof but does its job.  Examples and two accessory 
## functions are below.
unholtify <- function(df) {
	offset <- c(50,25,12.5)
	stopifnot('BLOCK' %in% names(df)) 
	df %>%
		separate(BLOCK, into=c('X','Y'), sep=1, remove=FALSE) %>%

	  ## convert the chars of BLOCK into [x,y] coordinates 
		## for the center of the area
	  mutate(X = factor(X, levels = as.character(3:8)),
	         Y = factor(Y, levels = LETTERS[10:3]),
	  			 WIDTH = 100) %>% 
	  mutate_at(vars(X, Y),
	            function(i) as.integer(i) * 100 - offset[1]) %>%

	  ## adjust coordinates for QUAD, if applicable
	  {if('QUAD' %in% names(.)) 
	  mutate(., X = case_when(QUAD %in% c("1", "3") ~ X - offset[2],
	                       QUAD %in% c("2", "4") ~ X + offset[2],
	                       TRUE ~ X),
	         Y = case_when(QUAD %in% c("1", "2") ~ Y + offset[2],
	                       QUAD %in% c("3", "4") ~ Y - offset[2],
	                       TRUE ~ Y),
	  			 WIDTH = ifelse(is.na(QUAD), WIDTH, 50))
	  else . } %>%
	
	  ## further adjust coordinates for SUBQUAD, if applicable
	  {if('SUBQUAD' %in% names(.)) 
	  mutate(., X = case_when(SUBQUAD %in% c("1", "3") ~ X - offset[3],
	                       SUBQUAD %in% c("2", "4") ~ X + offset[3],
	                       TRUE ~ X),
	         Y = case_when(SUBQUAD %in% c("1", "2") ~ Y + offset[3],
	                       SUBQUAD %in% c("3", "4") ~ Y - offset[3],
	                       TRUE ~ Y),
	  			 WIDTH = ifelse(is.na(SUBQUAD), WIDTH, 25))
	  else . } 
}

## This is a set of specs to add to a ggplot so the axes get
## marked with HRF study area row and column designations.
holtTicks <- list(
	  scale_x_continuous(breaks = seq(50, 550, by = 100),
	                     labels = as.character(3:8),
	                     expand = c(0, 0)) ,
	  scale_y_continuous(breaks = seq(50, 750, by = 100),
	                     labels = LETTERS[10:3],
	                     expand = c(0, 0)) ,
	  coord_equal(xlim = c(0, 600), ylim = c(0, 800))
)

## Beyond this point is all examples and helper functions.
{ stop("unholtify(); see unholtifier.R for examples.")
### Playing with Timber Inventory
## To run the examples download the 1996 HRF Timber Inventory data from 
## https://www.uvm.edu/femc/data/archive/project/hrftimberinv
## Quadratic mean and basal area functions are below the plot examples.
TInv_1996 <- read_csv("TInv_1996.csv",
	col_types = cols(COND = col_character(), DBH = col_integer(), 
		FIASP = col_character(), QUAD = col_character(), SPEC = col_character(),
		SUBQUAD = col_character(), TREENO = col_character(), YEAR = col_integer()))

## Mean DBH
TInv_1996 %>% 
	unholtify() %>%
	group_by(BLOCK, QUAD, SUBQUAD) %>%
	mutate(plotStat = mean(DBH, na.rm=TRUE)) %>%
	ggplot(aes(x = X, y = Y, fill = plotStat)) +
		geom_tile(aes(height = WIDTH, width = WIDTH)) +
		labs(title='Mean DBH', fill='cm') +
		holtTicks

## Quad Mean DBH
TInv_1996 %>% 
	unholtify() %>%
	group_by(BLOCK, QUAD, SUBQUAD) %>%
	mutate(plotStat = rms(DBH)) %>%
	ggplot(aes(x = X, y = Y, fill = plotStat)) +
		geom_tile(aes(height = WIDTH, width = WIDTH)) +
		labs(title='Quadratic Mean DBH', fill='cm') +
		holtTicks

## BA
TInv_1996 %>% 
	unholtify() %>%
	mutate(tBA = bArea(DBH)) %>%
	group_by(BLOCK, QUAD, SUBQUAD) %>%
	mutate(plotStat = sum(tBA, na.rm=TRUE)) %>%
	ggplot(aes(x = X, y = Y, fill = plotStat)) +
		geom_tile(aes(height = WIDTH, width = WIDTH)) +
		labs(title='Basal Area', fill='m^2') +
		holtTicks


### Accessories
## Basal area
bArea <- function (dbh, sys='m', ...) {
	fc <- ifelse(sys == 'm', .00007854, .005454154)
	dbh^2 * fc
}
# bArea(14, 'e')

## Quadratic mean
rms <- function(x) {	sqrt(mean(x^2, na.rm=TRUE)) }
}