PontiR: What the Package Does (One Line, Title Case)

# Load devtools
library(devtools)

# Install Thomas's package from github repo
# install_github('bilintoh/timecomptest', build_vignettes = TRUE)

# Load other libraries
library(timecomptest)
library(rgdal)
library(dplyr)
library(magrittr)
library(raster)
library(reshape)
library(ggnewscale)
library(graphics)

# set path to working directory
wd <- setwd(
  "F:/DATA/MApComp/working"
  ) # make this relative and point to MBM and MODIS data

# import all raster files from the work directory as a list
data <- list.files(
  path = wd,
  pattern = '_cy.rst$',
  all.files = TRUE,
  full.names = FALSE
  )

mbm <- data[1:19]
modis <- data[20:38]

raster_test <- raster(data[1])
crs(raster_test) <- 4326

ext <- extent(-6000000, -5100000, -1700000, -880000)

cropped_test <- crop(raster_test, ext)
# cropped_rst <- c()

cropped_mbm <- sapply(mbm, function(x) {
  # print(x)
  r <- raster(x)
  crs(r) <- 4326
  crop(r, ext)
  })

cropped_modis <- sapply(modis, function(x) {
  # print(x)
  r <- raster(x)
  crs(r) <- 4326
  crop(r, ext)
})

# list(cbind(
#   x = c(-7000000, -6100000, -6100000, -7000000), 
  # y = c(0, 900000, 900000, 0))))
# pol <- st_geometry(pol) %>% st_as_sf() %>% cbind(ID = 1:nrow(.), .)

# stack raster files
mbmStack <- stack(cropped_mbm)
modisStack <- stack(cropped_modis)

# make a data frame from the raster brick
df_mbm <- as.data.frame(mbmStack, xy = TRUE) # anything beyond 2000 rows and 2000 cols would be slow to run
df_modis <- as.data.frame(modisStack, xy = TRUE)

# make a list of data frames using the catTrajectory function
?timecomptest::catTrajectory # look at documentation

# mbm
traj_data_mbm <- timecomptest::catTrajectory(
  dfRaster = df_mbm,
  noData = 0, # default
  category = 1 # default
)

# modis
traj_data_modis <- timecomptest::catTrajectory(
  dfRaster = df_modis,
  noData = 0, # default
  category = 1 # default
)


# Trajectories:
# 0. Mask
# 1. Absence
# 2. Presence
# 3. Gain
# 4. Loss

# This function is used to generate a map of trajectories created in first data
# frame of list trajectory_data
?mapTrajectories

# Create a vector of pixel resolution
pix_size_mbm <- c(yres(cropped_mbm[[1]]), xres(cropped_mbm[[1]])) # change this to reflect pixel size
pix_size_modis <- c(yres(cropped_modis[[1]]), xres(cropped_modis[[1]])) # change this to reflect pixel size


# Create the map using mapTrajectories
mapTrajectories(
  rasList = traj_data_mbm,
  pixelRes = pix_size_mbm,
  mapTitle = 'MAP OF CHANGE: MAPBIOMAS',
  labx = 'X(m)',
  laby = 'Y(m)',
  axisFont = 1.5,
  titleFont = 1.3,
  legendFont = 1
)

# Use the stackedBarData function to create data for creating the stacked bar chart
stk_data_mbm <- stackedBarData(
  cattrajectory = traj_data_mbm,
  timePoints = seq.int(2001,2019),
  pixelSize = pix_size_mbm
)

stk_data_modis <- stackedBarData(
  cattrajectory = traj_data_modis,
  timePoints = seq.int(2001,2019),
  pixelSize = pix_size_modis
)

# Use the StackedTrajectories function to create the stacked bars
StackedTrajectories(
  stk_data_mbm,
  Category = 'FOREST',
  xAxisText = 12,
  yAxisText = 12,
  axisTitle = 14,
  plotTitle = 16,
  legendText = 12
) + 
  ggplot2::labs(x = 'Time Interval', y = 'Change (thousand km^2)') + 
  ggplot2::guides(x = ggplot2::guide_axis(angle = 90)) +
  ggplot2::ggtitle('Change involving Forest during 18 time intervals from 2001 to 2019 (MapBiomas)')

StackedTrajectories(
  stk_data_modis,
  Category = 'FOREST',
  xAxisText = 12,
  yAxisText = 12,
  axisTitle = 14,
  plotTitle = 16,
  legendText = 12
) + 
  ggplot2::labs(x = 'Time Interval', y = 'Change (thousand km^2)') + 
  ggplot2::guides(x = ggplot2::guide_axis(angle = 90)) +
  ggplot2::ggtitle('Change involving Forest during 18 time intervals from 2001 to 2019 (MODIS)')

# browseVignettes('timecomptest')