Time-Series Matrix (TSMx): A visualization tool for plotting multiscale temporal trends

<b>Time-Series Matrix (TSMx): A visualization tool for plotting multiscale temporal trends</b>
<br><br>
TSMx is an R script that was developed to facilitate multi-temporal-scale visualizations of time-series data. The script requires only a two-column CSV of years and values to plot the slope of the linear regression line for all possible year combinations from the supplied temporal range. The outputs include a time-series matrix showing slope direction based on the linear regression, slope values plotted with colors indicating magnitude, and results of a Mann-Kendall test. The start year is indicated on the <i>y</i>-axis and the end year is indicated on the <i>x</i>-axis. In the example below, the cell in the top-right corner is the direction of the slope for the temporal range 2001–2019. The red line corresponds with the temporal range 2010–2019 and an arrow is drawn from the cell that represents that range. One cell is highlighted with a black border to demonstrate how to read the chart—that cell represents the slope for the temporal range 2004–2014. This publication entry also includes an excel template that produces the same visualizations without a need to interact with any code, though minor modifications will need to be made to accommodate year ranges other than what is provided.
<br><br>
<i>TSMx for R was developed by Georgios Boumis; TSMx was originally conceptualized and created by Brad G. Peter in Microsoft Excel.</i>
<br><br>
<i>Please refer to the associated publication</i>:
<br>
Peter, B.G., Messina, J.P., Breeze, V., Fung, C.Y., Kapoor, A. and Fan, P., 2024. Perspectives on modifiable spatiotemporal unit problems in remote sensing of agriculture: evaluating rice production in Vietnam and tools for analysis. <i>Frontiers in Remote Sensing</i>, 5, p.1042624.
<br>
<a href="https://www.frontiersin.org/journals/remote-sensing/articles/10.3389/frsen.2024.1042624">https://www.frontiersin.org/journals/remote-sensing/articles/10.3389/frsen.2024.1042624</a>
<br><br>
<img src ="https://github.com/cartoscience/seagul/blob/main/tsmx/TSMx_Sample-Chart.png?raw=true" width="1000" height="auto"</img>
TSMx sample chart from the supplied Excel template. Data represent the productivity of rice agriculture in Vietnam as measured via EVI (enhanced vegetation index) from the NASA MODIS data product (MOD13Q1.V006).
<br><br>
<b>TSMx</b> R script:
<pre>
# import packages

library(dplyr)
library(readr)
library(ggplot2)
library(tibble)
library(tidyr)
library(forcats)
library(Kendall)

options(warn = -1) # disable warnings

# read data (.csv file with "Year" and "Value" columns)

data <- read_csv("EVI.csv")

# prepare row/column names for output matrices

years <- data %>% pull("Year")
r.names <- years[-length(years)]
c.names <- years[-1]
years <- years[-length(years)]

# initialize output matrices

sign.matrix <-
matrix(data = NA,
nrow = length(years),
ncol = length(years))
pval.matrix <-
matrix(data = NA,
nrow = length(years),
ncol = length(years))
slope.matrix <-
matrix(data = NA,
nrow = length(years),
ncol = length(years))

# function to return remaining years given a start year

getRemain <- function(start.year) {
years <- data %>% pull("Year")
start.ind <- which(data[["Year"]] == start.year) + 1
remain <- years[start.ind:length(years)]
return (remain)
}

# function to subset data for a start/end year combination

splitData <- function(end.year, start.year) {
keep <-
which(data[['Year']] >= start.year & data[['Year']] <= end.year)
batch <- data[keep,]
return(batch)
}

# function to fit linear regression and return slope direction

fitReg <- function(batch) {
trend <- lm(Value ~ Year, data = batch)
slope <- coefficients(trend)[[2]]
return(sign(slope))
}

# function to fit linear regression and return slope magnitude

fitRegv2 <- function(batch) {
trend <- lm(Value ~ Year, data = batch)
slope <- coefficients(trend)[[2]]
return(slope)
}

# function to implement Mann-Kendall (MK) trend test and return significance
# the test is implemented only for n>=8

getMann <- function(batch) {
if (nrow(batch) >= 8) {
mk <- MannKendall(batch[['Value']])
pval <- mk[['sl']]
} else {
pval <- NA
}
return(pval)
}

# function to return slope direction for all combinations given a start year

getSign <- function(start.year) {
remaining <- getRemain(start.year)
combs <- lapply(remaining, splitData, start.year = start.year)
signs <- lapply(combs, fitReg)
return(signs)
}

# function to return MK significance for all combinations given a start year

getPval <- function(start.year) {
remaining <- getRemain(start.year)
combs <- lapply(remaining, splitData, start.year = start.year)
pvals <- lapply(combs, getMann)
return(pvals)
}

# function to return slope magnitude for all combinations given a start year

getMagn <- function(start.year) {
remaining <- getRemain(start.year)
combs <- lapply(remaining, splitData, start.year = start.year)
magns <- lapply(combs, fitRegv2)
return(magns)
}

# retrieve slope direction, MK significance, and slope magnitude

signs <- lapply(years, getSign)
pvals <- lapply(years, getPval)
magns <- lapply(years, getMagn)

# fill-in output matrices

dimension <- nrow(sign.matrix)

for (i in 1:dimension) {
sign.matrix[i, i:dimension] <- unlist(signs[i])
pval.matrix[i, i:dimension] <- unlist(pvals[i])
slope.matrix[i, i:dimension] <- unlist(magns[i])
}

sign.matrix <- data.frame(sign.matrix)
pval.matrix <- data.frame(pval.matrix)
slope.matrix <- data.frame(slope.matrix)

# rename rows/columns of output matrices

names(sign.matrix) <- c.names
row.names(sign.matrix) <- r.names
names(pval.matrix) <- c.names
row.names(pval.matrix) <- r.names
names(slope.matrix) <- c.names
row.names(slope.matrix) <- r.names

# pretty-print slope direction output matrix

sign.matrix <-
sign.matrix %>% rownames_to_column() %>% gather(colname, value, -rowname)
names(sign.matrix) <- c("Start", "End", "Sign")
sign.matrix[['Color']] <- NA
sign.matrix [['Trend']] <- NA
pos.ind <- which(sign.matrix[['Sign']] > 0)
neg.ind <- which(sign.matrix[['Sign']] < 0)
sign.matrix[['Color']][pos.ind] <- "#99CCFF"
sign.matrix[['Color']][neg.ind] <- "#FF9966"
sign.matrix[['Trend']][pos.ind] <- "+"
sign.matrix[['Trend']][neg.ind] <- "-"
sign.matrix[['Trend']][is.na(sign.matrix[['Trend']])] <- ""

ggplot(sign.matrix,
aes(x = End,
fct_rev(Start),
col = Color,)) +
scale_color_identity() +
scale_fill_identity() +
coord_fixed() +
geom_tile(aes(fill = Color), color = 'black') + xlab("end year") + ylab("start year") +
labs(title = "Slope direction", caption = "Boumis and Peter (2021)") +
theme(
axis.text.x = element_text(angle = 15),
axis.text.y = element_text(angle = 15),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face = "italic"),
plot.caption = element_text(face = "italic")
) + scale_x_discrete(position = "top") +
geom_text(aes(label = Trend), color = "black")

# save .png of slope direction output matrix

ggsave(
"Slope-Direction.png",
device = "png",
unit = "in",
height = 7,
width = 12,
dpi = 300
)

# pretty-print MK significance output matrix

pval.matrix <-
pval.matrix %>% rownames_to_column() %>% gather(colname, value, -rowname)
names(pval.matrix) <- c("Start", "End", "Pval")
pval.matrix[['Color']] <- NA
pval.matrix [['Sign']] <- NA
sign.ind <- which(pval.matrix[['Pval']] < 0.05)
nsign.ind <- which(pval.matrix[['Pval']] > 0.05)
pval.matrix[['Color']][sign.ind] <- "#99FFFF"
pval.matrix[['Color']][nsign.ind] <- "#CC33CC"
pval.matrix[['Sign']][sign.ind] <- "<0.05"
pval.matrix[['Sign']][nsign.ind] <- ">0.05"
pval.matrix[['Sign']][is.na(pval.matrix[['Pval']])] <- ""

ggplot(pval.matrix,
aes(x = End,
fct_rev(Start),
col = Color,)) +
scale_color_identity() +
scale_fill_identity() +
coord_fixed() +
geom_tile(aes(fill = Color), color = "black") + xlab("end year") + ylab("start year") +
labs(title = "Mann-Kendall trend significance", caption = "Boumis and Peter (2021)") +
theme(
axis.text.x = element_text(angle = 15),
axis.text.y = element_text(angle = 15),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face = "italic"),
plot.caption = element_text(face = "italic")
) + scale_x_discrete(position = "top") +
geom_text(aes(label = Sign), color = "black", size = 2)

# save .png of MK significance output matrix

ggsave(
"MK-Significance.png",
device = "png",
unit = "in",
height = 7,
width = 12,
dpi = 300
)

# pretty-print slope magnitude output matrix

slope.matrix <-
slope.matrix %>% rownames_to_column() %>% gather(colname, value, -rowname)
names(slope.matrix) <- c("Start", "End", "Magnitude")

ggplot(slope.matrix, aes(End, fct_rev(Start), fill = Magnitude)) +
xlab("end year") + ylab("start year") +
geom_tile(color = "black") +
coord_fixed() +
labs(title = "Slope magnitude", caption = "Boumis and Peter (2021)") +
scale_fill_gradient2(
low = "#FF3300",
mid = "#FFFFFF",
high = "#0033FF",
midpoint = 0.0
) +
theme(
axis.text.x = element_text(angle = 15),
axis.text.y = element_text(angle = 15),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face = "italic"),
plot.caption = element_text(face = "italic")
) + scale_x_discrete(position = "top") +
guides(fill = guide_colourbar(
title = "",
barwidth = 0.5,
barheight = 7.5
))

# save .png of slope magnitude output matrix

ggsave(
"Slope-Magnitude.png",
device = "png",
unit = "in",
height = 7,
width = 12,
dpi = 300
)
</pre>
<b>Sample outputs:</b>
<img src ="https://github.com/cartoscience/seagul/blob/main/tsmx/TSMx_Slope-Direction_Small.png?raw=true" width="1000" height="auto"</img>
<img src ="https://github.com/cartoscience/seagul/blob/main/tsmx/TSMx_Slope-Magnitude_Small.png?raw=true" width="1000" height="auto"</img>
<img src ="https://github.com/cartoscience/seagul/blob/main/tsmx/TSMx_MK-Significance_Small.png?raw=true" width="1000" height="auto"</img>
<br>
<b>Project information:</b><br>
SEAGUL: Southeast Asia Globalization, Urbanization, Land and Environment Changes
<br>
<a href="http://seagul.info/">http://seagul.info/</a><br>
<a href="https://lcluc.umd.edu/projects/divergent-local-responses-globalization-urbanization-land-transition-and-environmental">https://lcluc.umd.edu/projects/divergent-local-responses-globalization-urbanization-land-transition-and-environmental</a>
<br>
This project was made possible by the the NASA Land-Cover/Land-Use Change Program (Grant #: 80NSSC20K0740)
<br><br>
Another publication developed by the SEAGUL research group that is related to the kind of modifiable temporal unit problems (MTUP) presented here can be viewed at: <a href="https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/M4ZGXP">MSZSI: Multi-Scale Zonal Statistics [AgriClimate] Inventory</a>.