TrajectoryAnalysisTutorial.Rmd

---
title: "A How-to Guide for Trajectory Analysis"
author: "Rufai Omowunmi Balogun,  Harrison Leduc, Daniel Neau, Sunita Phuyal"
date: "2024-11-13"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# A tutorial on how to run an end to end trajectory analysis on categorical remote sensing data

## General Overview
This document showcases with examples, the different steps a user need to take to effectively implement the Trajectory Analysis of long time series of categorical datasets.

The popular method in land change modelling typically focuses on capturing changes between time points with uniform durations. These kinds of modelling captures are limited 
because they are not able to capture short term changes within shorter period. The trajectory analysis is a method developed to capture the changes in a category through 
time across different time intervals and sites whether the duration and sites are uniform or not. These changes in categories are captured in terms of the unified size on 
each site, the alternation and exchange between categories.

We show a use case of the evolution of Mangroves Landcover in the Mahakam Delta, Indonesia. This location is a hotspot for Shrimp Pond Aquaculture and hence the transitions 
(exchange) exists between the Mangrove Land Cover Category and the Pond Aquaculture Category. This document mostly build on the [timeseriesTrajectories](https://github.com/bilintoh/timeseriesTrajectories) 
R package and the concepts developed by Bilintoh et al in "Bilintoh, T.M., Pontius Jr, R,.G., & Zhang, A., (2024). Methods to compare sites concerning a Category’s change during various time intervals. 
GIScience & Remote Sensing, 61 (1), 14. https://doi.org/10.1080/15481603.2024.2409484"

In addition to this document, we have prepared an [explainer video](ADD LINK TO EXPLAINER VIDEO HERE) that contains the key concepts necessary for understanding trajectory analysis for remote sensing 
applications as well as interpreting the output charts produced by the Trajectory Analysis R package. We anticipate this tutorial to be useful for 1.) students (undergraduates and graduates), 
2.) Researchers, 3.) Professionals, and 4.) Anyone curious about the concepts and applications of trajectory analysis.


### Requirements
Since this tutorial was developed in R programming language, we expect the users to have familiarity with the programming language and tools like RTools and RStudio for setting up the development environment.For a quick refresher on the R programming language, see this [crash course on R programming language](https://www.youtube.com/watch?v=ZYdXI1GteDE) and the [Extra Resources](#extra-resources) section of this document. 

## Step 1: Setting up your development environment
### I. Download and install R and RStudio
Start by installing the [R programming](https://posit.co/download/rstudio-desktop/) language on your workstation along with the [RStudio IDE](https://posit.co/download/rstudio-desktop/). Select the installers that matches your Operating System. When working with R, you can think of RStudio as the interface and R as the bones.

### II. Download and install RTools
Select and download the RTools that matches the version of R that you installed in the previous step. See the [Extra Resources](#extra-resources) section of this document for 
The RTools is required to properly install the R packages used in this tutorial. Here is the [link](https://cran.r-project.org/bin/windows/Rtools/) to the RTools page: <https://cran.r-project.org/bin/windows/Rtools/>
Click the blue text 'Rtools44 installer'. After the installer has been downloaded click it and go through the wizard. We suggest leaving everything as it is, this will give you a simple base install.

### III. Install the required packages

Run the following:

``` 
install.packages('devtools')
devtools::install_github('bilintoh/timeseriesTrajectories')
```
The second line will give a prompt with additional directions: when prompted type '1' and hit enter, this will update all necessary packages that [timeseriesTrajectories](https://github.com/bilintoh/timeseriesTrajectories) uses.

We also use the `terra` package, run: 

```
install.packages('terra')
```

## Step 2: Data Preparation
This sections showcases the steps needed to go from a multi-class categorical data to a binary data format required to run a trajectory analysis.

## Background on our example data
We are reviewing the Mahakam Delta in Indonesia. Specifically, how mangroves transition to Shrimp ponds. Having trajectories of the Mahakam Mangroves can help give insight to what the future land may look like, impacting shrimp farmers, locals, and the ecosystem as a whole. The [timeseriesTrajectories](https://github.com/bilintoh/timeseriesTrajectories) R package example has pre-made binary files, which are already in TIF format. However, users do not usually start with a well-prepared datasets as in the example (where we started with a multiclass dataset in a .rst format). Hence, we showcase a simple process for preparing the datasets to capture our category of interest and into the expected data format for  [timeseriesTrajectories](https://github.com/bilintoh/timeseriesTrajectories) 

```{r echo=TRUE}
## load packages
library(terra)

# four categories in our datasets::: 1. Mangrove, 2. Pond Aquaculture, 3. Water, and 4. Other LandCover
## load datasets
ras1999 <- rast("data/Mahakam_1999.rst")
ras2014 <- rast("data/Mahakam_2014.rst")
ras2018 <- rast("data//Mahakam_2018.rst")
ras2020 <- rast("data//Mahakam_2020.rst")
ras2022 <- rast("data//Mahakam_2022.rst")
```

### We start with these Multiclass datasets across 5 timepoints

```{r echo=FALSE}
mahakam_rasters <- list(ras1999, ras2014, ras2018, ras2020, ras2022)
years <- c(1999, 2014, 2018, 2020, 2022)
par(mfrow = c(2, 3))
for (i in seq_along(mahakam_rasters)) {
  coltab(mahakam_rasters[[i]]) <- data.frame(values=c(1, 3, 4, 5), cols=rainbow(4))
  plot(mahakam_rasters[[i]], main = paste("Year:", years[i]))
}
par(mfrow = c(1, 1))

```

### I. Ensure all the datasets are in the right projection

In order to get the best representation of the category change, it is recommended to use Equal Area Projection such as the Cylindrical Equal Area when re-projecting the data.
The Cylindrical Equal-Area Projection is a flexible projection that requires defining the central meridian and standard parallels to minimize distortion for the region of 
interest.

```{r echo=TRUE}
mahakam_rasters <- list(ras1999, ras2014, ras2018, ras2020, ras2022)
for (i in seq_along(mahakam_rasters)) {
    print(crs(mahakam_rasters[[i]]))
}
```
It is not clear from this what the Coordinate Reference System is, hence we will assign the CRS for the Mahakam Delta and reproject to the to each of our data input to ensure they overlap properly.

```{r echo=TRUE}
cyl_equal_area_crs <- "+proj=cea +lon_0=120 +lat_ts=0 +datum=WGS84 +units=m +no_defs"
for (i in seq_along(mahakam_rasters)) {
    crs(mahakam_rasters[[i]]) <- cyl_equal_area_crs
}

```

```{r echo=TRUE}
mahakam_rasters_reprojected <- lapply(mahakam_rasters, function(raster) {
    project(raster, cyl_equal_area_crs)
})

```
### II. Binarize the yearly images
```{r echo=TRUE}
# extract categories
categories1999 <- unique(values(mahakam_rasters_reprojected[[1]])) # 1999
categories1999 <- categories1999[!is.na(categories1999)]

categories2014 <- unique(values(mahakam_rasters_reprojected[[2]])) # 2014
categories2014 <- categories2014[!is.na(categories2014)]

categories2018 <- unique(values(mahakam_rasters_reprojected[[3]])) # 2018
categories2018 <- categories2018[!is.na(categories2018)]

categories2020 <- unique(values(mahakam_rasters_reprojected[[4]])) # 2020
categories2020 <- categories2020[!is.na(categories2020)]

categories2022 <- unique(values(mahakam_rasters_reprojected[[5]])) # 2022
categories2022 <- categories2022[!is.na(categories2022)]

# generate binary maps per year
category_layers2022 <- lapply(categories2022, function(cat) {
  layer <- mahakam_rasters_reprojected[[5]] == cat
  names(layer) <- paste0("category2022_", cat)
  return(layer)
})

category_layers2020 <- lapply(categories2020, function(cat) {
  layer <- mahakam_rasters_reprojected[[4]] == cat
  names(layer) <- paste0("category2020_", cat)
  return(layer)
})

category_layers2018 <- lapply(categories2018, function(cat) {
  layer <- mahakam_rasters_reprojected[[3]] == cat
  names(layer) <- paste0("category2018_", cat)
  return(layer)
})

category_layers2014 <- lapply(categories2014, function(cat) {
  layer <- mahakam_rasters_reprojected[[2]] == cat
  names(layer) <- paste0("category2014_", cat)
  return(layer)
})

category_layers1999 <- lapply(categories1999, function(cat) {
  layer <- mahakam_rasters_reprojected[[1]] == cat
  names(layer) <- paste0("category1999_", cat)
  return(layer)
})
```



### III. Stack the category of interest and save
This step saves the yearly binary images for mangrove into a stack and save as a TIFF file, ready for analysis using the [timeseriesTrajectories](https://github.com/bilintoh/timeseriesTrajectories) R package. 

```{r echo=TRUE}

# stack mangroves categories across five time points and save to tiff
mangroves_stack <- c(category_layers1999[[2]], 
                     category_layers2014[[3]], # mangrove category was stored as 3
                     category_layers2018[[2]],
                     category_layers2020[[2]],
                     category_layers2022[[2]]
                     )

writeRaster(mangroves_stack,
          "data/mangroves_stacked.tif",
          overwrite = TRUE)

```

```{r echo=FALSE}
years <- c(1999, 2014, 2018, 2020, 2022)
names(mangroves_stack) <- years

mangrove_colors <- c("blue", "red")

par(mfrow = c(2, 3),  mar = c(4, 4, 4, 12))
for (i in 1:nlyr(mangroves_stack)) {
  plot(mangroves_stack[[i]], col = mangrove_colors, 
       main = paste("Mangrove Presence in", years[i]),
       axes = TRUE,
  )
}
```


## Step 3: Run trajectory Analysis
After properly preparing our datasets, we can now run trajectory analysis for our class of interest -- Mangroves and see its evolution across time and space.

### I. Load the R package

```{r echo=TRUE}
library("timeseriesTrajectories")
```

### II. Load stacked binarized mangrove dataset

```{r echo=TRUE}
rastStack <- rast("data/mangroves_stacked.tif")
```

### Compute the number of presence of change and plot

```{r echo=TRUE}
vert_units <- "number of pixels"
dataPreview(rastStack,
            timepoints = years,
            vertunits = vert_units,
            xAngle = 0)

num_pres_change <- presenceData(rastStack, nodata = -999)

presencePlot(num_pres_change,
             pltunit = "m",
             dataEpsg = 2933,
             scalePos = "bottomright",
             narrowPos = "topleft",
             narrowSize = 1,
             categoryName = "Mangroves",
             xAxis = "Longitude (m)",
             yAxis = "Latitude (m)",
             axisText = 0.7,
             axisLabel = 0.7,
             plotTitle = 0.8)

```

### Compute and visualize mangrove Trajectories 

```{r echo=TRUE}
traj_data <- rastertrajData(rastStack,
                            zeroabsence = 'yes')
```


```{r echo=TRUE}
## visualize the trajectories
trajPlot(traj_data,
         axisShow = "no",
         categoryName = "Mangroves",
         narrowPos = NA,
         scalePos = NA,
         scaleSize = 1.2,
         axisText = 1.2,
         axisLabel = 1.2,
         plotTitle = 1.2,
         legendTex = 0.8,
         xAxis = "Longitude (m)",
         yAxis = "Latitude (m)",
         downsample = TRUE)
```


### Compute the annual change in Mangrove presence defined by Unified size

```{r echo=TRUE}
stackbar_data <- rasterstackData(rastStack,
                                 timePoints = years,
                                 spatialextent = 'unified',
                                 zeroabsence = 'yes',
                                 annualchange = 'yes',
                                 categoryName = 'variable',
                                 regionName = 'Mahakam Delta',
                                 varUnits = "(Square Kilometers)",
                                 constant  = 1)

stackbarPlot(stackbar_data,
             axisSize = 10,
             lbAxSize = 10,
             lgSize = 7.5,
             titleSize = 12,
             datbreaks = "no",
             upperlym = 35,
             lowerlym = - 50,
             lymby = 5,
             upperlym2 = 0.5,
             lymby2 = 0.1,
             xAngle = 0)
```



## Open Questions
- Evaluation of trajectory analysis,
- Multi-class trajectory analysis,

## Extra Resources {#extra-resources}
1. [R Programming Crash Course](https://www.youtube.com/watch?v=ZYdXI1GteDE) 
2. [How to Install R and RStudio on Windows](https://www.youtube.com/watch?v=9SzKJH93t5o)
3. [How to Install R and RStudio on Mac](https://www.youtube.com/watch?v=I5WIMX4LK8M)
4. [How to Install RTools](https://www.youtube.com/watch?v=VEvw43iF6rY)


## Contributors
1. [Rufai Omowunmi Balogun](https://github.com/Ruphai), Graduate School of Geography, Clark University
2. [Harrison Leduc](https://github.com/HLed12), Department of Geography, Clark University
3. [Daniel Neau](https://github.com/danieljneau), Department of Geography, Clark University
4. [Sunita Phuyal](https://github.com/MeSunita), Graduate School of Geography, Clark University


## References
1. Bilintoh, T.M., Pontius Jr, R,.G., & Zhang, A., (2024). Methods to compare sites concerning a Category’s change during various time intervals. GIScience & Remote Sensing, 61 (1), 14. https://doi.org/10.1080/15481603.2024.2409484
2. Center for Geospatial Analytics, [Coastal Habitat Mapping Portal](https://www.clarku.edu/centers/geospatial-analytics/projects/pond-aquaculture-and-its-impact-on-mangroves-and-other-coastal-wetlands/)
3. Explainer Video for Bilintoh et al. 2024. (by Mangrove3)