CourseR

# Install necessary packages and load them
if(!"terra" %in% installed.packages()){install.packages("terra")}
if(!"sf" %in% installed.packages()){install.packages("sf")}

library(terra)
library(sf)

# Create data directory if it does not yet exist
if (!dir.exists("data")) {
  dir.create("data")
}

# Download to data directory
download.file(url = 'https://github.com/GeoScripting-WUR/VectorRaster/releases/download/tutorial-data/landsat8.zip', destfile = 'data/landsat8.zip')

# Unzip to data directory
unzip('data/landsat8.zip', exdir = "data")

# Identify the right file
landsatPath <- list.files(path = "data/", pattern = glob2rx('LC8*.tif'), full.names = TRUE)
wagLandsat <- rast(landsatPath)

wagLandsat

names(wagLandsat)

wagLandsat[wagLandsat < 0] <- NA
# wagLandsat[[1]][wagLandsat[[1]] < 0] <- NA # 操作某层
# 对第1、2、3层进行操作
# for(i in 1:3) {
#   wagLandsat[[i]][wagLandsat[[i]] < 0] <- NA
# }

# Band names can be changed here
names(wagLandsat) <- c("band11", "band22", "band33", "band44", "band55", "band66", "band77")
names(wagLandsat)

# Select which bands to assign to the red, green, and blue colour channels
plotRGB(wagLandsat, 5, 4, 3)

# Download and unzip the Dutch municipalities boundaries
download.file(url = 'https://github.com/GeoScripting-WUR/Scripting4GeoIntro/releases/download/data/data.zip', destfile = 'data/GADM.zip')
unzip('data/GADM.zip', exdir = "data")

# Load level 2 dataset
nlCitySf <- st_read("data/gadm41_NLD_2.json")

names(nlCitySf)
print(head(nlCitySf))

dim(nlCitySf)

# Remove rows with NA
nlCitySf <- nlCitySf[!is.na(nlCitySf$NAME_2),]
dim(nlCitySf)

# Filter Wageningen
wagContour <- nlCitySf[nlCitySf$NAME_2 == 'Wageningen',]
dim(wagContour)

# Get the target CRS from the wagLandsat raster object
targetCRS <- st_crs(wagLandsat)
targetCRS

# Use sf to transform wagContour to the same projection as wagLandsat
wagContourUTM <- st_transform(wagContour, targetCRS)

wagLandsatCrop <- crop(wagLandsat, wagContourUTM)
wagLandsatSub <- mask(wagLandsat, wagContourUTM)

# Set graphical parameters (one row and two columns)
opar <- par(mfrow = c(1, 2))
plotRGB(wagLandsatCrop, 5, 4, 3)
plotRGB(wagLandsatSub, 5, 4, 3)
plot(st_geometry(wagContourUTM), add = TRUE, border = 'green', col = 'transparent', lwd = 3) # set fill colour to transparent

# Reset graphical parameters
par(opar)

download.file(url = 'https://github.com/GeoScripting-WUR/VectorRaster/releases/download/tutorial-data/wageningenWater.zip', destfile = 'data/wageningenWater.zip')
unzip('data/wageningenWater.zip', exdir = "data")

# Load the Water Shapefile directly as an sf object
water <- st_read('data/Water.shp')

names(water) # check attribute columns

waterUTM <- st_transform(water, targetCRS)

plot(waterUTM)

# Mask the Landsat data
wagLandsatSubW <- mask(wagLandsatSub, mask = waterUTM, inverse = TRUE)
# inverse = FALSE：遮蔽掩膜区域内的数据，保留区域外的数据。
# inverse = TRUE：遮蔽掩膜区域外的数据，保留区域内的数据。

# Plot
plotRGB(wagLandsatSubW, 5, 4, 3)
# plot(st_geometry(waterUTM), col = 'lightblue', add = TRUE, border = '#3A9AF0', lwd = 1) # use a site such as https://htmlcolorcodes.com/ to find the hexidecimal code for colours

# Create output directory if it does not yet exist
if (!dir.exists("output")) {
  dir.create("output")
}

# Export the simple feature to a KML. 
try(st_write(water, "output/water.kml", driver = "kml", delete_dsn = ifelse(file.exists("output/water.kml"), TRUE, FALSE)))

# Note that the code above checks for existence of the KML file and sets the delete_dsn option to true or false accordingly

# # Install and load library
# if(!"mapview" %in% installed.packages()){install.packages("mapview")}
# if(!"webshot" %in% installed.packages()){install.packages("webshot")}
# library(mapview)

# # Plot the water polygons
# map <- mapview(water)
# map

# # Save to file
# mapshot(map, url = 'mapview.html')

# Download to data directory
download.file(url = 'https://github.com/GeoScripting-WUR/VectorRaster/releases/download/tutorial-data/belgium_elev.zip', destfile = 'data/belgium_elev.zip')

# Unzip to data directory
unzip('data/belgium_elev.zip', exdir = "data")

# Download data
bel <- rast('data/belgium_elev.tif')

# Display metadata
bel

plot(bel)

# Example to perform system set-up checks
if(!"terra" %in% installed.packages()){install.packages("terra")}
library(terra)

gdal()

# Generate a SpatRaster
r <- rast(ncol = 40, nrow = 20)
class(r)

r

# Using the previously generated SpatRaster
# Let's first put some values in the cells of the layer
r[] <- rnorm(n = ncell(r))  # 将随机生成的值赋给 r 这个栅格图层中的每个单元格
# rnorm() 是 R 中的一个函数，用于生成正态分布的随机数。它会生成一组符合正态分布的随机数。
# ncell(r) 是一个函数，用来返回栅格对象 r 中的单元格数量。
# 即，这会告诉你栅格有多少个单元格（像素）
# rnorm(n = ncell(r)) 会生成与栅格单元格数量一样多的随机数。

# Create a SpatRaster with 3 layers
s <- c(r, r*2, r)
# 创建一个由三层组成的 SpatRaster 对象 s。c() 是 R 中的一个函数，用来连接多个对

# Let's look at the properties of the resulting object
s

s2 <- list(r, r*2, r)

# Let's look at the properties of the resulting object
s2

getwd()

# Create data directory if it does not yet exist
if (!dir.exists("data")) {
  dir.create("data")
}

# Download the data
# In case the download code doesn't work, use method = 'wget'
download.file(url = 'https://github.com/GeoScripting-WUR/IntroToRaster/releases/download/tahiti/gewata.zip', destfile = 'data/gewata.zip')

# Unpack the archive
unzip('data/gewata.zip', exdir = "data")

# When passed without arguments, list.files() returns a character vector, listing the content of the working directory
list.files()

# To get only the files with .tif extension in the data directory
list.files('data', pattern = glob2rx('*.tif'))

# Or if you are familiar with regular expressions
list.files('data', pattern = '^.*\\.tif$')

gewata <- rast('data/LE71700552001036SGS00_SR_Gewata_INT1U.tif')

gewata

# Again, make sure that your working directory is properly set
getwd()

# Download and unzip the data
download.file(url = 'https://github.com/GeoScripting-WUR/IntroToRaster/releases/download/tahiti/tura.zip', destfile = 'data/tura.zip')
unzip(zipfile = 'data/tura.zip', exdir = 'data')

# Retrieve the content of the tura sub-directory
turaList <- list.files(path = 'data/tura/', full.names = TRUE)

plot(rast(turaList[1]))

turaStack <- rast(turaList)
turaStack

ext(turaStack)

# Create output directory if it does not yet exist
if (!dir.exists("output")) {
  dir.create("output")
}

# Write this file to the 
writeRaster(x = turaStack, filename = 'output/turaStack.tif', names = list.files(path = 'data/tura/'), datatype = "INT2S")

ndvi <- (gewata[[4]] - gewata[[3]]) / (gewata[[4]] + gewata[[3]])

plot(ndvi)

# Define the function to calculate NDVI using app()
ndviApp <- function(x) {
    ndvi <- (x[[4]] - x[[3]]) / (x[[4]] + x[[3]])
    return(ndvi)
}
ndvi2 <- app(x = gewata, fun = ndviApp)

mean <- global(ndvi, fun = 'mean')
mean

# Get only the value from the dataframe
mean$mean

# Calculate the amount of NA values (which should be zero in this case)
global(ndvi, fun = 'isNA')

# Calculate the amount of non-NA values
global(ndvi, fun = 'notNA')

# Which, as we have no NA values, should be equal to
ncell(ndvi)

# One single line is sufficient to project any raster to any CRS
ndviLL <- project(ndvi, 'epsg:4326')

# Since this function will write a file to your working directory
# you want to make sure that it is set where you want the file to be written
# It can be changed using setwd()
getwd()

# Note that we are using the filename argument, contained in the ellipsis (...) of
# the function, since we want to write the output directly to file.
writeRaster(x = ndviLL, filename = 'output/gewataNDVI.tif')

download.file(url = 'https://github.com/GeoScripting-WUR/IntroToRaster/releases/download/tahiti/tahiti.zip', destfile = 'data/tahiti.zip')
unzip(zipfile = 'data/tahiti.zip', exdir = 'data')

# Load the data as a SpatRaster and investigate its contents
tahiti <- rast('data/LE70530722000126_sub.tif')
tahiti

plot(tahiti)

# Display names of each individual layer
names(tahiti)

plotRGB(tahiti, 3, 4, 5, stretch = "lin")  # 线性

plot(tahiti, 7)

# Extract cloud layer from the SpatRaster
cloud <- tahiti[[7]]

# Replace 'clear land' with 'NA'
cloud[cloud == 0] <- NA

# Plot the stack and the cloud mask on top of each other
plotRGB(tahiti, 3, 4, 5, stretch = "lin")
plot(cloud, add = TRUE, legend = FALSE)

# Extract cloud mask layer
fmask <- tahiti[[7]]

# Create a subset of the first six layers
tahiti6 <- subset(tahiti, 1:6)

# Perform value replacement
tahiti6[fmask != 0] <- NA

# First define a value replacement function
cloud2NA <- function(x) {
    x[1:6][x[7] != 0] <- NA
    return(x)
}

# Apply the function on the two SpatRasters
tahitiCloudFree <- app(x = tahiti, fun = cloud2NA)

# Visualize the output
plotRGB(tahitiCloudFree, 3, 4, 5, stretch = "lin")

if(!"sits" %in% installed.packages()){install.packages("sits")}
if(!"geojsonsf" %in% installed.packages()){install.packages("geojsonsf")}
if(!"stars" %in% installed.packages()){install.packages("stars")}
if(!"tmap" %in% installed.packages()){install.packages("tmap")}
if(!"httr2" %in% installed.packages()){install.packages("httr2")}
library(sits)

sits_list_collections(c("AWS", "MPC", "CDSE", "USGS"))

tahiti84 <- project(tahiti, "EPSG:4326")
(roi <- ext(tahiti84))