

## ---- Load packages and functions ----

require(data.table)
source('SnowpackVisFunctions.R')



## ---- Load sample data ----

## Profile date
ProfileDate <- '2018-01-08'

## Read data
Profiles <- readRDS('SampleProfiles.rds')
ProfileTimeline <- readRDS('SampleProfileTimeline.rds')



## ---- Prepare data for visualizations ----

## Colour palettes
colIACS <- data.frame(grain_type = c('SH', 'DH', 'PP', 'DF', 'RG', 'FCxr', 'FC', 'MFcr', 'MF', 'IF'),
                      colour = c('#FF00FF', '#0000FF', '#00FF00', '#228B22', '#FFB6C1', '#ADD8E6', '#ADD8E6', '#FF0000', '#FF0000','#00FFFF'), stringsAsFactors = F)
colNEW <- data.frame(grain_type = c('SH', 'DH', 'PP', 'DF', 'RG', 'FCxr', 'FC', 'MFcr', 'MF', 'IF'),
                     colour = c('#ee3a1d', '#4678e8', '#ffde00', '#f1f501', '#ffccd9', '#dacef4', '#b2edff', '#addd8e', '#d5ebb5', '#a3ddbb'), stringsAsFactors = F)
colNEW4 <- data.frame(grain_type = c('SH', 'DH', 'PP', 'DF', 'RG', 'FCxr', 'FC', 'MFcr', 'MF', 'IF'),
                      colour = c('#95258f', '#95258f', '#ffde00', '#ffde00', '#dacef4', '#dacef4', '#dacef4', '#d5ebb5', '#d5ebb5', '#d5ebb5'), stringsAsFactors = F)
grainDict <- colNEW

## Number of profiles
nProfiles <- length(Profiles)

## rbind into large table where each row is a layer
ProfileTable <- proRbind(Profiles)

## Count grain types by age
DistinctStnGrains <- unique(ProfileTable[c('station_id','deposition_date','grain_type')])
Grains <- table(DistinctStnGrains$grain_type, DistinctStnGrains$deposition_date) / nProfiles * 100
NonPrstGrain <- -Grains[!(rownames(Grains) %in% c('SH','DH')),]
PrstGrain <- Grains[rownames(Grains) %in% c('SH','DH'),]

## Factor into terrain bins
TerrainBins <- c('ALP N', 'ALP E', 'ALP S', 'ALP W', 'TL N', 'TL E', 'TL S', 'TL W', 'BTL N', 'BTL E', 'BTL S', 'BTL W')
ProfileTable$terrain <- factor(paste(ProfileTable$band, ProfileTable$dir), ordered = T, levels = TerrainBins )

## Derive relative measure for sensitivity to triggers
ProfileTable$sensitivity <- exp(-ProfileTable$ssi)

## Re-order to put important grains on top of plots
ProfileTable <- ProfileTable[order(factor(ProfileTable$grain_type, ordered = T, levels = c('RG','MF','MFcr','FCxr','FC','DF','PP','DH','SH'))),] ## put important grains on top

## Assign grain colours
ProfileTable$col <- getColoursGrainType(ProfileTable$grain_type)

## Rank deposition dates for plotting
ProfileTable$ageRank <- data.table::frank(as.numeric(as.Date(ProfileTable$deposition_date)), ties.method = 'dense')




## ---- Fig 2: Compare colour palettes ----

ProfileDay <- ProfileTimeline[[which(proSummary(ProfileTimeline)$date == ProfileDate)]]

layout(matrix(c(1,1,2,1,1,2,3,3,4,3,3,4,5,5,6,5,5,6), nrow = 6, ncol = 3, byrow = T))
par(mar = c(3,5,1,1), las = T)
ylim <- c(0, 110)

grainDict <- colIACS
plotTSSnowProfile(ProfileTimeline, DateEnd = ProfileDate, ylab = 'Height (cm)')
legend('topleft', legend = grainDict$grain_type, col = grainDict$colour, pch = 15, bty = 'n', cex = 0.6)
plot(ProfileDay, ylim = ylim, TempProfile = F); box()

grainDict <- colNEW
plotTSSnowProfile(ProfileTimeline, DateEnd = ProfileDate, ylab = 'Height (cm)')
legend('topleft', legend = grainDict$grain_type, col = grainDict$colour, pch = 15, bty = 'n', cex = 0.6)
plot(ProfileDay, ylim = ylim, TempProfile = F); box()

grainDict <- colNEW4
plotTSSnowProfile(ProfileTimeline, DateEnd = ProfileDate, ylab = 'Height (cm)')
legend('topleft', legend = grainDict$grain_type, col = grainDict$colour, pch = 15, bty = 'n', cex = 0.6)
plot(ProfileDay, ylim = ylim, TempProfile = F); box()
grainDict <- colNEW




## ---- Fig 3: Side-by-side profiles ----

par(mfrow = c(1,1), las = T)
plotSideBySideProfiles(Profiles, SortMethod = 'hs', TopDown = T, ylab = 'Depth (cm)')
legend('bottom', legend = colNEW$grain_type, col = colNEW$colour, pch = 15, horiz = T, cex = 0.7, bty = 'n')




## ---- Fig 4: Profile aggregate ----

## Draw barplots
par(mfrow = c(1,1), las = T, mar = c(4,6,1,5))
barplot(NonPrstGrain, col = getColoursGrainType(rownames(NonPrstGrain)), horiz = T, border = NA, space = 0,
        xlim = c(min(colSums(NonPrstGrain)), max(colSums(PrstGrain))),
        xlab = 'Percentage of profiles', xaxt = 'n', yaxt = 'n')
barplot(PrstGrain, col = getColoursGrainType(rownames(PrstGrain)), horiz = T, border = NA, space = 0, add = T, yaxt = 'n', xaxt = 'n', xlab = '')

## Add x-axis with positive number labels on both sides
Xlabels <- pretty(c(min(colSums(NonPrstGrain)), max(colSums(PrstGrain))), n = 4)
axis(1, at = Xlabels, labels = abs(Xlabels))
abline(v = 0)
box()

## Add deposition date labels
DateArray <- as.Date(colnames(Grains))
Saturdays <- which(weekdays(DateArray, abbreviate = F) == 'Saturday')
SaturdayLabels <- format(DateArray[weekdays(DateArray, abbreviate = F) == 'Saturday'], '%b %d')
Saturdays <- Saturdays[2:length(Saturdays)] ## drop first one for clean plot
SaturdayLabels <- SaturdayLabels[2:length(SaturdayLabels)] ## drop first one for clean plot
axis(2, at = Saturdays, labels = SaturdayLabels)
abline(h = Saturdays, lty = 3, col = 'grey')
title(ylab = 'Deposition date', mgp = c(5,1,0))

## Add age labels
Age <- as.Date(ProfileDate) - DateArray
SaturdayAge <- Age[Saturdays]
axis(4, at = Saturdays, labels = SaturdayAge)
text(par("usr")[2]*1.5, mean(par("usr")[3:4]), "Age (days)", srt = 90, xpd = TRUE, pos = 4)

legend('bottomleft', legend = colNEW$grain_type, col = colNEW$colour, pch = 15, cex = 0.7, bty = 'n')


## ---- Fig. 5: Terrain bin jitter plots ----

par(mfrow = c(1,1), las = T)
plot(-depth ~ jitter(as.numeric(terrain), amount = 0.5), ProfileTable, col = ProfileTable$col,
     xaxt = 'n', xlab = '', yaxt = 'n', ylab = 'Depth (cm)',
     pch = 20, cex = 0.4)
axis(1, at = 1:length(TerrainBins), labels = TerrainBins, tick = F, cex.axis = 0.6)
abline(h = 0, v = seq(0.5, length(TerrainBins) + 0.5, by = 1), col = 'gray')
DepthGrid <- pretty(c(-max(ProfileTable$hs), 0), n = 4)
axis(2, at = DepthGrid, labels = -DepthGrid)
abline(h = DepthGrid, lty = 3, col = 'dark grey')
legend('bottom', legend = colNEW$grain_type, col = colNEW$colour, pch = 15, horiz = T, cex = 0.7, bg = 'white')



## ---- Fig. 6: Terrain bins scaled by sensitivity ----

par(mfrow = c(1,1), las = T)
plot(-depth ~ jitter(as.numeric(terrain), amount = 0.5), ProfileTable, col = ProfileTable$col,
     cex = ProfileTable$sensitivity*2,
     pch = 16, xaxt = 'n', xlab = '', yaxt = 'n', ylab = 'Depth (cm)')
axis(1, at = 1:length(TerrainBins), labels = TerrainBins, tick = F, cex.axis = 0.6)
abline(h = 0, v = seq(0.5, length(TerrainBins) + 0.5, by = 1), col = 'gray')
DepthGrid <- pretty(c(-max(ProfileTable$hs), 0), n = 4)
axis(2, at = DepthGrid, labels = -DepthGrid)
abline(h = DepthGrid, lty = 3, col = 'dark grey')
legend('bottom', legend = colNEW$grain_type, col = colNEW$colour, pch = 15, horiz = T, cex = 0.7, bg = 'white')



## ---- Fig. 7: Likelihood plots ----

par(mfrow = c(1,2), mar = c(4,5,0,0))

## Spatial distribution
barplot(NonPrstGrain, col = getColoursGrainType(rownames(NonPrstGrain)), horiz = T, border = NA, space = 0,
        xlim = c(min(colSums(NonPrstGrain)), max(colSums(PrstGrain))),
        xlab = 'Percentage of profiles', xaxt = 'n', yaxt = 'n')
barplot(PrstGrain, col = getColoursGrainType(rownames(PrstGrain)), horiz = T, border = NA, space = 0, add = T, yaxt = 'n', xaxt = 'n', xlab = '')
Xlabels <- pretty(c(min(colSums(NonPrstGrain)), max(colSums(PrstGrain))), n = 4)
axis(1, at = Xlabels, labels = abs(Xlabels))
abline(v = 0)
box()
DateArray <- as.Date(colnames(Grains))
Saturdays <- which(weekdays(DateArray, abbreviate = F) == 'Saturday')
SaturdayLabels <- format(DateArray[weekdays(DateArray, abbreviate = F) == 'Saturday'], '%b %d')
axis(2, at = Saturdays, labels = SaturdayLabels)
abline(h = Saturdays, lty = 3, col = 'grey')
title(ylab = 'Deposition date', mgp = c(4,1,0))

## Sensitivity
plot(ageRank ~ sensitivity, ProfileTable,
     col = ProfileTable$col, pch = 16, cex = 0.2, xlim = c(0, 0.5),
     yaxt = 'n', ylab = '', xlab = 'Relative sensitivity to triggers', bty = 'n')
abline(h = Saturdays, lty = 3, col = 'grey')
box()
legend('bottomright', legend = colNEW$grain_type, col = colNEW$colour, pch = 15, cex = 0.7, bty = 'n')
