# ------------------------------------------------------#
#                                                       #
#  GIS Using R: terra package                           #
#  Dennis Milechin                                      #
#  2/26/2023                                            #
#                                                       #
#  Evaluation Link:                                     #
#  http://scv.bu.edu/survey/tutorial_evaluation.html    #
#                                                       #
# ------------------------------------------------------#


# Terra is a replacement for an existing package called raster
# Benefits of terra
#   1. functions are usually more computationally efficient
#   2. lets you work on large raster datasets that are too 
#      large to fit into the main memory.



########################################################
## Suggested Reading
##
## Discussion of terra and raster.
## https://geocompr.robinlovelace.net/spatial-class.html#an-introduction-to-terra
########################################################

#----------------------------------------#
#  Install required packages          ####
#----------------------------------------#

# Install sf package 
install.packages("terra")
install.packages("spDataLarge", repos = "https://nowosad.r-universe.dev")
install.packages("tidyterra")
install.packages("spData")
install.packages("tidyvers")

#----------------------------------------#
#  Expectations for this tutorial     ####
#----------------------------------------#

# There are many functions in this package, 
# we are unable to cover them all.  List of 
# functions available:
# https://rspatial.github.io/terra/reference/index.html


#----------------------------------------#
#  Let's explore terra                   ####
#----------------------------------------#
library(terra) 
library(spData)
library(spDataLarge)
library(tidyterra)
library(ggplot2)


##############################
# RASTER Loading and Plotting
##############################

raster_filepath = system.file( "raster/nz_elev.tif", package="spDataLarge" )

# First lets look at the meta data, or the header, of the raster data
describe( raster_filepath )

# To load data or create a new raster data set, we
# use the rast() function.
my_rast = rast( raster_filepath )

class( my_rast )

# Let's look at the my_rast object
my_rast

# We can plot the data
plot( my_rast )

# If we have "tidyterra" installed we can use ggplot
ggplot( data=my_rast ) +
  geom_spatraster()

# This doesn't work, we need to specify data in the geom_spatraster()
# function.
ggplot() +
  geom_spatraster( data=my_rast )

# Let's take a look at the arguments for geom_spatraster
?geom_spatraster



ggplot() +
  geom_spatraster( 
                data=my_rast, 
                #show.legend = NA,
                #maxcell=50
                )

# If we want to adjust the colors we use one
# of the scale fill functions.  Tidyterra provides us with some
# fill functions to use: https://dieghernan.github.io/tidyterra/reference/index.html#scales

# Functions are groups into three categories
#  scale_*_terrain_d(): For discrete values.
#  scale_*_terrain_c(): For continuous values.
#  scale_*_terrain_b(): For binning continuous values.


ggplot() +
  geom_spatraster( data=my_rast ) +
  scale_fill_terrain_c() +
  theme_bw()

?scale_fill_terrain_c


ggplot() +
  geom_spatraster( data=my_rast ) +
  scale_fill_hypso_tint_c() +
  theme_bw()

?scale_fill_hypso_tint_c

ggplot() +
  geom_spatraster( data=my_rast ) +
  scale_fill_hypso_tint_c( 
     #palette="moon",
     #direction=-1,
     #alpha=0.5,
     #na.value="black",
     #limits=c(100, 1000),
     #breaks=c(0, 1000, 2000, 3000, 4000, 5000, 6000, 7000),
     #guide = guide_colorbar(
     #    direction="horizontal",
     #    title.position = "top",
     #     barwidth = 15
     #  ),
    ) +
  # labs(fill = "elevation (m)") +
  # theme(legend.position = "bottom") +
  theme_bw() 
   


# We can also create contour plots
# Examples from: https://dieghernan.github.io/tidyterra/reference/geom_spat_contour.html
volcanoe_path <- system.file("extdata/volcano2.tif", package = "tidyterra")
r_volc <- rast(volcanoe_path)


ggplot() +
  geom_spatraster_contour( data=r_volc )

# Again we can fine tune things with function attributes
ggplot() +
  geom_spatraster_contour(
    data = r_volc, aes(color = after_stat(level)),
    binwidth = 1,
    linewidth = 0.4
  ) +
  scale_color_gradientn(
    colours = hcl.colors(20, "Inferno"),
    guide = guide_coloursteps()
  ) +
  theme_minimal()

# There is also a fill variation of this function:
ggplot() +
  geom_spatraster_contour_filled(
    data = r_volc, 
    breaks = seq(80, 200, 10)) +
  scale_fill_hypso_d()

# One can also combine both of them in one plot
ggplot() +
  geom_spatraster_contour_filled(
    data = r_volc, breaks = seq(80, 200, 10),
    alpha = .7
  ) +
  geom_spatraster_contour(
    data = r_volc, breaks = seq(80, 200, 2.5),
    color = "grey30",
    linewidth = 0.1
  ) +
  scale_fill_hypso_d()


#################################
# Working with Multi-Band Raster
#################################

# Let's look at a multi raster file
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")

# Let's look at the header for this file.
describe(multi_raster_file)

multi_rast = rast(multi_raster_file)

class(multi_rast)

multi_rast

# Check how many layers exist in the raster file:
nlyr(multi_rast)

# And their names
names(multi_rast)

# Let's run it through the plot
plot(multi_rast)

# Use ggplot
ggplot() +
  geom_spatraster(data=multi_rast)

# Notice by default it will plot all layers at once.
# But the warning message tells us we can plot by
# 1 layer, or to use facet_wrap function.

ggplot() +
  geom_spatraster( data=multi_rast, aes( fill=landsat_1 ))

# We can try facet_wrap next

ggplot() +
  geom_spatraster( data=multi_rast) +
  facet_wrap(~lyr, ncol = 2)


# If we want to extract a specific layer
# we can use the subset() function.  We 
# can use an index value or the name of the layer.

multi_rast3 <- subset(multi_rast, 3)
multi_rast3

multi_rast4 <- subset(multi_rast, "landsat_4")
multi_rast4

# Or we can use the "$" method
multi_rast1 <- multi_rast$landsat_1


# We can also combine layers into one object:
multi_rast341 <- c(multi_rast3, multi_rast4, multi_rast1)
multi_rast341


# In the previous examples, the dimensions and resolution were identical.
# What if they are not?

# Let's resample the data into a new raster with different dimensions.

# First lets define the extent, or area of interest
# based on the existing raster data.
AOI <- ext(multi_rast3)

# We also need to assign the coordinate reference system to 
# the new raster.  We can extract that using crs() function
crs(multi_rast3)



# Let's create a template raster that we will resample into.
templ_rast <- rast( 
                   nrows=25, 
                   ncols=30, 
                   extent=AOI, 
                   crs=crs(multi_rast3) 
                   )

# Could we use SF for extracting and assign CRS?
library(sf)
st_crs(multi_rast3)

rast( 
  nrows=25, 
  ncols=30, 
  extent=AOI, 
  crs=st_crs(multi_rast3)
  ) 

# Unload the "sf" package.
detach("package:sf", unload=TRUE)
st_crs()

# Its not in the correct format for the rast function.
# So it is best to use the terra provided CRS functions.


resampled_rast <- resample(multi_rast3, templ_rast)

# Let's visually compare the rasters

ggplot() +
  geom_spatraster( data=multi_rast3) +
  ggtitle("Original")

ggplot() +
  geom_spatraster( data=resampled_rast) +
  ggtitle("Resampled")

# Let's compare the header information
multi_rast3
resampled_rast

# Let's try to combine
combined <- c(multi_rast3, resampled_rast)

# It complains the dimensions don't match.
# Lesson learned, the dimensions need to match to combine!


###########################
# A closer look at rast()
###########################

?rast

# Let's create some raster files from scratch.

# Specify number of columns and rows
rast1 <- rast(
               nrows=10,
               ncols=10,
               vals=runif(100, 0, 800)
             ) 

rast1

# Notice this defaults to an extent for the entire world.

ggplot() +
  geom_spatraster( data=rast1) +
  # geom_sf( data=world ) +
  ggtitle("Rast1")


# Let's limit the extent by using xmin and ymin.
rast2 <- rast(
            nrows=10,
            ncols=10,
            xmin=0,
            ymin=0,
            vals=runif(100, 0, 80)
          ) 

rast2

ggplot() +
  geom_spatraster( data=rast2) +
  geom_sf( data=world ) +
  ggtitle("Rast2")

# Alternatively we can specify a resolution and terra
# will fill in the area with appropriate number of
# columns and rows
rast3 <- rast(
  resolution=c(10,10),
  vals=runif(100, 0, 800)
)

rast3


# Without a constraint on the extent, it will do this for the entire
# world.
ggplot() +
  geom_spatraster( data=rast3) +
  geom_sf( data=world ) +
  ggtitle("Rast3")

# We can constrain it again by specifying extent.
rast4 <- rast(
  resolution=c(10,10),
  xmin=-50,
  ymin=-50,
  xmax=0,
  ymax=0,
  vals=runif(10, 0, 800)
) 

rast4

ggplot() +
  geom_spatraster( data=rast4) +
  geom_sf( data=world ) +
  ggtitle("Rast4")


# We can also extract an extent of an existing object
# that may represent the area of interest, like we did before
# with multi_rast3 object

nz_elev_path <- system.file( "raster/nz_elev.tif", package="spDataLarge" )
nz_elev <- rast(nz_elev_path)
AOI <- ext(nz_elev)

rast5 <- rast(
              nrows=10,
              ncols=10,
              extent=AOI,
              vals=runif(50, 0, 800),
              crs= crs(nz_elev)
              )

rast5

# Since the Area of Interest is New Zealand, lets
# plot only the AOI.
x_extent <- c(xmin(rast5), xmax(rast5))
y_extent <- c(ymin(rast5), ymax(rast5))


ggplot() +
  geom_spatraster( data=rast5) +
  geom_sf(data=world, fill=NA, color="orange", linewidth=1) +
  coord_sf( xlim = x_extent, ylim = y_extent) +
  ggtitle("Rast5")

# One of the data conditioning steps of raster data analysis is
# re sampling all the raster datasets onto the same raster grid with 
# the same CRS.  
# You can use properties of an existing one as a template:

template_rast <- rast(
                    res= res(rast5),
                    extent= ext(rast5),
                    crs= crs(rast5)
                   )

nz_elev_resampled <- resample(nz_elev, template_rast)

# We can plot the result
ggplot() +
  geom_spatraster( data=nz_elev_resampled) +
  geom_sf(data=world, fill=NA, color="orange", linewidth=1) +
  coord_sf( xlim = x_extent, ylim = y_extent) +
  ggtitle("Rast5")

# Now that the layers are aligned, we can stack them.
nz_stack <- c(nz_elev_resampled, rast5)

ggplot() +
  geom_spatraster( data=nz_stack) +
  geom_sf(data=world, fill=NA, color="orange", linewidth=1) +
  coord_sf( xlim = x_extent, ylim = y_extent) +
  facet_wrap(~lyr, ncol = 2)

# The resample() has some noteworthy arguments.
# method - to choose the resampling method.
# filename - Output result to a file rather then memory
?terra::resample

# Let's look at the "Source" property of nz_elev_resampled
nz_elev_resampled

# Notice its "memory"
# No resample with a "filename" defined
nz_elev_resampled2 <- resample(
                              nz_elev,
                              template_rast,
                              filename="test.tif",
                              overwrite=TRUE 
                              )

nz_elev_resampled2
# Notice now the data lives on the disk rather then in memory.

###########################
# RASTER CRS Transformation
###########################

# The approach to transforming raster data into another CRS
# is more complicated, as described in the link below:
# https://rspatial.org/spatial/6-crs.html#transforming-raster-data


###########################
# RASTER Data Exploration
###########################
# We can gets stats on this raster layer using summary function.

summary(nz_elev)

?summary

summary(nz_elev, size=1000)

# We can also apply other summarizing functions, such as:
# max() - standard deviation
# min() - minimum value
# sum() - maximum value
# range() - min and max value
# Median() - Median value

my_rast_max <- max(nz_elev)
my_rast_max

# Take a look at the dimensions, why do we still have several rows and columns?


# We need to use the global() function to return a single value. Otherwise
# it the function will be applied to each individual cell in the raster.
global(nz_elev, max)

# Now we get an NA value, this suggests there are NA values in the raster
global(nz_elev, max, na.rm=TRUE)

# We can try out some other aggregation functions:
global(nz_elev, min, na.rm=TRUE)
global(nz_elev, sum, na.rm=TRUE)
global(nz_elev, range, na.rm=TRUE)


# We can also create common statistics plots of the values
hist(nz_elev)
density(nz_elev)
boxplot(nz_elev)


###########################
# Accessing Cell Values
###########################

sample_rast <- rast(
  nrows=10,
  ncols=10,
  vals=1:10,
)

ggplot() +
  geom_spatraster( data=sample_rast)


# You can extract all values using values() function
values(sample_rast)

# Can also use indexing techniques
sample_rast[1,] # Row indexing
sample_rast[2,]
sample_rast[,1] # Column indexing

sample_rast[2:4, 2:4] # Block of values

# We can also use cell numbers to extract values

cells <- cellFromRowCol(sample_rast, col=2, row=c(2,3,4))
cells

sample_rast[cells]

# You can also extract xy coordinates from cell numbers:
xy <- xyFromCell(sample_rast, cells)
xy

# You can also extract values based on xy coordinates using extrac()
# function
extract(sample_rast, xy)

###########################
# Updating cell values
###########################

# You can update existing cell values with assignment operator
sample_rast[2:4, 2:4] <- 100

sample_rast
values(sample_rast)

ggplot() +
  geom_spatraster( data=sample_rast)

# We can also specify logical conditions on how we want to update 
# the cells.  Here I am seeting all cells with value of 100 to "NA"
sample_rast[sample_rast == 100] <- NA

sample_rast
values(sample_rast)

# We can plot this to see a gray box representing NA values.
ggplot() +
  geom_spatraster( data=sample_rast)



###########################
# Raster Algebra  
###########################

# Let's create a raster with values of 1
rast_val_1 <- rast(
  nrows=5,
  ncols=5,
  vals=1,
  names="layer1"
)

values(rast_val_1)

# Let's create another raster with values of 2
rast_val_2 <- rast(
  nrows=5,
  ncols=5,
  vals=2,
  names="layer2"
)

values(rast_val_2)

# Let's add the the rasters together
rast_result1 <- rast_val_1 + rast_val_2
values(rast_result1)

# We can create more complicated equations
rast_result2 <- rast_val_1 / (rast_val_2 + 100)
values(rast_result2)

# If we have a stack of rasters, we can still apply raster math
rast_stack <- c(rast_val_1, rast_val_2)
rast_stack

rast_rasult3 <- rast_stack$layer1 * 5 + rast_stack$layer2 
values(rast_rasult3)

# We can rename the layer to something more meaningful
names(rast_rasult3) <- "result3"
rast_rasult3





##############################
# Terra Supports Vector Data  
##############################

# Learn more about SpatVector here: https://rspatial.org/spatial/7-vectmanip.html
# We can read in a vector data set using vect() function

f <- system.file("ex/lux.shp", package="terra")
p <- vect(f)
p

class(p)

ggplot()+
  geom_spatvector(data=p, aes(fill=NAME_2))


# Extract specific column
name_2 <- p$NAME_2
name_2
class(name_2)  # Notice it loses the SpatVector class

name_2 <- p[, "NAME_2"]
name_2
class(name_2)  # Notice this time SpatVector class is retained

# We can also convert SpatVector to an "sf" object
library(sf)
sf_p <- st_as_sf(p)
sf_p
class(sf_p)

ggplot()+
  geom_sf(data=sf_p)

# We can also convert it back to SpatVector
spatvect_p <- vect(sf_p)
spatvect_p
class(spatvect_p)

# Terra includes a collection of overlay functions as well, found here:
#   https://rspatial.org/spatial/7-vectmanip.html#overlay
#

# One can convert a vector data to raster using terra
# An example taken from: ?terra::rasterize

f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)

ggplot() +
  geom_spatvector( data=v, aes(fill=NAME_2))

r <- rast(v, ncols=75, nrows=100)
z <- rasterize(v, r, "NAME_2")

z
values(z)

ggplot() +
  geom_spatraster( data=z, aes(fill=NAME_2) )

# We can also convert this back to 
z_poly <- as.polygons(z)
z_poly
ggplot() +
  geom_spatvector( data=z_poly, aes(fill=NAME_2))

# We can use vector to crop a raster data set to a certain extent
elev_file <- system.file("ex/elev.tif", package="terra")
rast_elev <- rast(elev_file)

border_file <- system.file("ex/lux.shp", package="terra")
vect_border <- vect(border_file)

ggplot() +
  geom_spatraster( data=rast_elev) +
  geom_spatvector( data=vect_border, fill=NA, linewidth=1, color="orange" )

# Let's say we are only interested in specific areas as defined
# by the subset command below:
vect_subset <- vect_border[9:12,]

ggplot() +
  geom_spatraster( data=rast_elev) +
  geom_spatvector( data=vect_subset, fill=NA, linewidth=1, color="orange" )

# Apply crop
rast_crop <- crop(rast_elev, v[9:12,], mask=TRUE)

ggplot() +
  geom_spatraster( data=rast_crop) +
  geom_spatvector( data=vect_border, fill=NA, linewidth=1, color="orange" )


###### Helpful resources:
#  Spatial Data Science with R and "terra"
#  https://rspatial.org/index.html
#
#  "Geocomputation with R".
#  https://geocompr.robinlovelace.net/index.html
#

####################
## QUESTIONS?
####################
##
##  Evaluation form: http://rcs.bu.edu/eval
##
##  contact: help@scc.bu.edu
##