UsingRastersInR.Rmd

---
title: "Using Raster Data in R"
author: "G Perkins"
date: "28/08/2019"
output: ioslides_presentation
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## 1. Raster Basics 

- Loading in raster 
- CRS and projections 
- Raster Calculations
- Viewing and mapping Calculations 
- Exercise.....


## Packages history and tools : 

- raster (first released in 2010 Hijman)
- sp package also supports some special types of data SpatialGridDataFrame & SPatialpixalsDataFrame
- stars: Spatiotemporal Arrays, Raster and Vector Data Cubes (Pebesma)
- helper packages : fasterize, rgdal, mapview 
- Other tools and bridges (RSAGA, RQGIS, RPYGeon)


## Using Rasters in R : Work flow 

-  Many geospatial data work flows rely on raster data sets; this covers everything from elevation data, to remote sensed data, lidar and many others. 
- Often you need to work with both raster and vector data sets. 
- Today we will introduce how read in, check, visualise raster datasets and work through the common functions and how to use them. 


## Your Task: 
You are given co-ordinates for a study area - download the data and provide a summary of the elevation. 



## Load packages and set up data folder 

```{r set-up , echo = TRUE, message= FALSE}

library(raster)

# set up the data location 
data.dir <- "C:/Training/R_intermediate/data" ## this need to be fixed to relative path

# see what goodies are in the folder 
list.files(data.dir)

```

## Loading Rasters 

```{r loading rasters, echo = TRUE}

# read in a single raster 
dem <- raster(file.path(data.dir, "DEM.tif"))

# look at the fundamental raster information : 
# dimensions, cell size, crs , min and max
dem

# check resolution 
res(dem)


# check a summary of the data distribution 
summary(dem) # throws an error 
summary(dem, maxsamp = ncell(dem)) # force to check all raster 

# explore the data a little - we have lots of NAs
head(dem) 
values(dem)
unique(values(dem))


# dealing with NA data
dem <- setMinMax(dem) 

# set NA values 
NAvalue(dem) = -9999

# set NA values 

minValue(dem)
maxValue(dem)

# subset to only include the peaks 

dem.peaks <- clamp(dem, lower = 1000, useValues = TRUE )

#slightly different to 

dem.peaks  <- reclassify(dem , cbind(0 , 1000 , NA), right=FALSE)



```

## Visualise rasters 

```{r Visualise rasters, echo = TRUE}
#quick check of the raster 
plot(dem)

# more fancy exploring 
library(rasterVis)
histogram(dem)
plot3D(dem)

library(mapview)
mapview(dem)

# convert to data frame and plot 
library(ggplot2)
dem_df <- as.data.frame(dem, xy = TRUE)
str(dem_df)

# ggspatial # layer_spatial 

ggplot() +
    geom_raster(data = dem_df , aes(x = x, y = y, fill = DEM)) +
    scale_fill_viridis_c() +
    coord_quickmap()

#...... much more time consuming

```


## Check projections, assign and transform projections

```{r checking projections, echo = TRUE }

# Coordinate Reference System (CRS)
# we can have CRS in three formats (ESPG, Proj4string, WKT) 
# in proj4string the "+" seperates each element (similar to csv is seperated with a comma)

dem
crs(dem)

# Questions: what units are the data in ? 
# units = m (meters)

#proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 + x_0=1000000 + y_0=0 
#ellps=GRS80 
#towgs84=0,0,0,0,0,0,0
#units=m 
#no_defs 

# we can check the information on our ESPG:3005
CRS("+init=epsg:3005")

projection(dem)  # gives us the string version proj4string

# what about projected data ? 
# lets read in another raster 

ndvi <- raster(file.path(data.dir, "NDVI.tif"))
ndvi 

# what is the cell size of the raster? 
# what is the projection? 

projection(ndvi) # what is the crs for this raster - is it the same as the dem layer?

#[1] "+proj=utm +zone=9 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

# transform a projection to match our dem 
ndvi.3005 <- projectRaster(ndvi, dem) 

ndvi.3005

summary(ndvi.3005)
#              NDVI
#Min.    -0.2097639
#1st Qu.  0.5884543
#Median   0.6328611
#3rd Qu.  0.6891333
#Max.     0.8716997
#NA's     0.0000000


# What if we want to reclassify our rater and summarise the results 
library(dplyr)

dem.class <- reclassify (dem, c(-Inf, 800, 1,
                                700, 1000, 2, 
                                1000, 1250, 3, 
                                1250, Inf, 4))

plot(dem.class)


dem.class <- as.data.frame(dem.class) %>%
  filter(!is.na(DEM))

dem.class <- dem.class %>%
  group_by(DEM) %>%
  summarise(count = n())


```


# Raster Calculation

```{r}

# Raster Calculations 

# for small sized rasters 
rasters.out <- dem - ndvi.3005

# for larger sized rasters (processing intensive)
overlay() # for single layers
calc() # for stacks 


outputRaster <- overlay(dem, ndvi.3000, fun = functionName)


CHM_ov_HARV <- overlay(dem,
                       ndvi.3000,
                       fun = function(r1, r2) { return( r1 - r2) })


# write raster out 
# write out geotiff. 

writeRaster(CHM_ov_HARV, "CHM_HARV.tiff",
            format="GTiff",
            overwrite=TRUE,
            NAflag=-9999)

```


# Exercise 1: 

DEM hillshade and differences 

# read in a slope and aspect raster and create a hillshade 
# check the crs and cell size 


# bonus question: generate a contour map 

```{r}

slope <- raster(file.path(data.dir, "slope.tif")) 

aspect <- raster(file.path(data.dir, "aspect.tif"))

# create a hillshade 
dem.hill <- hillShade(slope, aspect, 
                      angle=40, 
                      direction=270)

plot(dem.hill,
     col=grey.colors(100, start=0, end=1),
     legend=F)


# export geotiff 
writeRaster(dem.hill,
            filename="demhill.tif",
            format="GTiff",
            options="COMPRESS=LZW",
            overwrite = TRUE,
            NAflag = -9999)

## bonus 

dem.contour <- rasterToContour(dem, maxpixels=100000)

plot(dem.contour, add = T)
plot(dem.contour)


```




## 2. Manipulating Rasters 


# What about if we want to stack multiple raster layers and extract values ? 
# How to we match up rasters with different extents and sizes? 




## Aggregate/disaggregate 

```{r matching multiple raters for use in raster calculater}



## STILL TO TIDY UP :


# Oftern we are looking at multiple raster : 
# many rasters from different datasets
# working with a stack or brick 

# for example lets 

# we now know our raster objects are in the same crs
crs(dem)
crs(ndvi.3005)

# what about the extent? 
extent(dem)
extent(ndvi.3005)

# all good matching extents!







################# UP TO HERE ####################################

 r <- raster(file.path(lidar.dir,layers.list[i])) 
  crs(r) <- "+init=epsg:3005"
  re <- crop(r, rast2.5)
  re <- disaggregate(re, fac = 4) # Factor 10m / 4 = 2.5m.
  res <- resample(re, rast2.5)
  #stack(res,rast2.5) 
  writeRaster(res, paste("D:/PEM_DATA/Data/Layers/Dec_2.5m", "/",
                            gsub(".asc",".tif",outList[i]), sep = ""), 
              overwrite = TRUE)
  print(paste("Downscale complete for: ",  outFolder, "/", outList[i], sep = ""))


# what if we add another layer 
twi <- raster(file.path(data.dir, "twi.tif"))

crs(twi)
extent(twi)

# same CRS
# different call size
# different extent

# lets convert to the same resolution 
twi.25 <- aggregate(twi, fact = 5, fun = mean, expand = TRUE)

# check resolution and extent
res(twi.25)
extent(twi.25)

r2 = merge(twi.25, dem) # throws an error as not aligned


# as the rasters do not match up we need to project values of twi onto a raster with format of dem 

twi.2 <- projectRaster(twi, dem, res = 25, method = "bilinear") 


# UP TO HERE -------------------------------
# need to match the extents 

dem.t <-trim(dem)

#Modifying a raster object (spatial extent) 

crop # geographic sibset of a larger raster object 
      # by extent of by another raster 

trim :  crops the outerlayesr containing NAs 
extend : opposite of trim - pads out the 


```


# Cropping Rasters (by raster or by vector) 

```{r}

# Crop rasters to an area of Interest (raster) 

library(sf)
library(dplyr)

aoi<- st_read(dsn = data.dir,layer = "CanyonCreekBoundary", quiet = TRUE) %>%
 st_transform(3005) 

plot(aoi)

# crop a single raster 
dem.c <- crop(dem, aoi)

# lets see what is looks like 
plot(dem)
plot(dem.c, add = T, col = "red")

```

# rasterise values 

```{r}
# convert polygon to raster and crop 

roads <- st_read(dsn = data.dir,layer = "Deception_Atlas_Roads_Cleaned", quiet = TRUE) 

roads <- roads %>%
  st_transform(3005) %>%  # read in classified roads
  select (ROAD_SURFA)

plot(roads)
roads <- st_buffer(roads, dist = 25) %>% st_cast("MULTIPOLYGON")

# convert this to a raster 
library(fasterize)
rroads <- fasterize(roads, dem, field = "ROAD_SURFA")###convert to raster

# crop to the aoi or could crop as a vector 
roads.c <- crop(rroads, aoi)

plot(st_geometry(aoi))
plot(rroads, add = T) # extends past the linear aoi as a raster 


# or 


library(bcmaps)

x <- available_layers()
bec <- get_layer("bec", class = "sf")

st_crs(bec)
st_crs(aoi)


# Crop to area of interest 
bec_aoi <- st_intersection(bec, aoi)

plot(bec_aoi)

# create a key to reclassify ##STILL TO DO .....

map.key <- as.data.frame(bec_aoi)


# convert to a raster 
bec_aoi <- fasterize(bec_aoi, dem.c, field = "MAP_LABEL") ###convert to raster



```



```{r}


Mask to value 






```






## Stack Rasters 

```{r}

# Raster stacks and Bricks 

raster.s <- stack(dem.c,  roads.c)

nlayers(raster.s)
nlayers(dem.c)

#Make a Threat brick for analysis
ThreatBrick <- stack(HumanDensityR,LivestockDensityR,HuntDDensR,SecureR,FrontCountryR,RdDensR)
names(ThreatBrick) <- c('HumanDensity','LivestockDensity','HunterDensity','SecureHabitat','FrontCountry','RoadDensity')
Threat_file <- file.path("tmp/ThreatBrick")
saveRDS(ThreatBrick, file = Threat_file)

```

# Extract values 

```{r}

# extract values 


# read in csv and convert to spatial points 
sites <- read.csv(file.path(data.dir, "Deception_samples.csv"))
# convert to spatial points 


raster::extract(sites, raster.s )



# some plots ?????

```



## Exercise 2: 





Exercise 2 (30 min) Gen 
Stack land cover, BEC and DEM 









# multiband rasters 

```{r}

# modisTools or MODIStsl

install.packages("MODISTools") or MODIS 

library(MODISTools)

mt_products()
mt_bands()


bands <- mt_bands(product = "MOD11A2")
head(bands)

mt_bands(product= Daymet)

mt_dates
mt_sites()

all.sites <- mt_sites()

all.sites %>% 
  filter(country == "Canada")

all.can <- all.sites[country == "Canada",]


```





```{r}

# sentinel imagery 


remotes::install_github("ranghetti/sen2r")

library(sen2r)

check_sen2r_deps()
sen2r()

s2_list(spatial_extent = aoi, time_interval = c(Sys.Date() - 5, Sys.Date()))


s2_list(spatial_extent = NULL, tile = NULL, orbit = NULL,
  time_interval = c(Sys.Date() - 10, Sys.Date()), time_period = "full",
  level = "auto", apihub = NA, max_cloud = 100,
  output_type = "vector")
Arguments



## lidar 
install.packages("lidR")

library(lidR)

las <- readLAS(file.path(data.dir,"Lidar", "bcalb_20160813.las"))

plot(las)

thr <- c(0,2,5,10,15)
edg <- c(0, 1.5)
chm <- grid_canopy(las, 1, pitfree(thr, edg))

plot(chm)














getwd()
download.file("https://pub.data.gov.bc.ca/datasets/177864/tif/bcalb/103k/103K016.zip", destfile = "temp")
unzip("temp")
works!


```

```{r}

install.packages("elevatr", dep = T)
library(elevatr)

library(sf)

aoi<- st_read(dsn = data.dir,layer = "CanyonCreekBoundary", quiet = TRUE) 
plot(aoi)

aoi <- 


aoi.sp <- as(aoi, 'Spatial')

strm <- get_elev_raster(aoi, z = 1) 



```





---

## References: 

* https://rspatial.org

* https://datacarpentry.org/r-raster-vector-geospatial/01-raster-structure/

* https://geocompr.robinlovelace.net/geometric-operations.html#geo-ras

* https://www.youtube.com/watch?v=yhpkx_xO-LE

* https://csgillespie.github.io/efficientR/



https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf

---