I've created a base map as seen in the document provided. In addition to this I've also made a combined raster file with Fs values that I wish to overlay over the basemap.
I was wondering how I can overlay one on top of the other. Both have the same spatial coordinates.
Do i need to merge the raw data or can i simply plot one ontop of the other.
library(raster)
library(exactextractr)
library(mapview)
library(tmap)
library(sf)
library(ggplot2)
theme_set(theme_bw())
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)
library(rgdal)
library(classInt)
library(RColorBrewer)
library(basemaps)
library(mapedit)
library(terra)
library(tidyterra)
library(dplyr)
library(scales)
library(geodata)
library(showtext)
library(maptools)
library(sp)
library(rgeos)
library(GISTools)
indir <-'C:/Assessment2'
rast <- raster('C:/Assessment2/Coweeta_slope.tif')
shplist <- list.files(indir, pattern = '*.shp',all.files = TRUE, full.names = TRUE, recursive = TRUE)
x <- lapply(shplist, shapefile)
m <- do.call(rbind, x)
m$slope <- exact_extract(rast, m, 'max')
spatial.variables <- as.data.frame(t(rbind(m$cohesion, m$soil_depth, m$soil_dens, m$hydro_cons, m$slope, m$frict_ang, m$slope)))
m$Fs = round(as.numeric(as.numeric(m$cohesion)/(as.numeric(m$soil_depth)*as.numeric(m$soil_dens)*9.81))+((as.numeric(m$hydro_cons)*cos(as.numeric(m$slope))*tan(as.numeric(m$frict_ang)))/sin(as.numeric(m$slope))),digits =2)
{r include=FALSE}
tmap_mode("view")
tm_shape(m) +
tm_polygons(col = "Fs")
# Preparting workspace
m
setwd('C:/Assessment2')
getwd()
dir.create('mydir')
# Converting polygons to .shp file
writeOGR(obj=m, dsn='mydir', layer= 'mshp', driver='ESRI Shapefile')
# Raster conversion process
Fshp <- readOGR(dsn='C:/Assessment2/mydir', layer='mshp')
# Creating raster skeleton, setting coordinate system and projection
Ftemp <- raster(ncols=8, nrows=15,xmn=275309 ,xmx=276835, ymn=3880601, ymx=3881923, crs='+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs')
# Populating Ftemp
Frast1 <- rasterize(Fshp, Ftemp, getcover=TRUE, progress = 'text')
Frast1
Frast <- as(Frast1, "SpatRaster")
r_init <- 'C:/Assessment2/DEM.tif'
map <- terra::rast(r_init)
names(map) <- "alt"
r_init
map
autoplot(map) +
theme_minimal()
## Creating hillshade effect
slope <- terrain(map, 'slope', unit = 'radians')
aspect <- terrain(map, 'aspect', unit = 'radians')
hill <- shade(slope, aspect, 30, 270)
# Normalising names
names(hill) <- 'shades'
# Hillshading, but we need a palette
pal_greys <- hcl.colors(1000, 'Grays')
ggplot() +
geom_spatraster(data = hill) +
scale_fill_gradientn(colors = pal_greys, na.value = NA)
# Using a vector of colors
index <- hill %>%
mutate(index_col = rescale(shades, to = c(1, length(pal_greys)))) %>%
mutate(index_col = round(index_col)) %>%
pull(index_col)
# Getting columns
vector_cols <- pal_greys[index]
# Avoiding resampling and not using aes
hill_plot <- ggplot() +
geom_spatraster(
data = hill, fill = vector_cols, maxcell = Inf,
alpha = 1
)
hill_plot
# Hypso gradient
grad_hypso <- hypso.colors2(10, 'dem_poster')
autoplot(map) +
scale_fill_gradientn(colours = grad_hypso, na.value = NA)
{r}
map_limits <- minmax(map) %>% as.vector()
# Rounded to lower and upper 500
map_limits <- c(floor(map_limits[1] / 500), ceiling(map_limits[2] / 500)) * 500
# Making min value to 0.
map_limits <- pmax(map_limits, 0)
# Comparing
minmax(map) %>% as.vector()
map_limits
{r}
base_plot <- hill_plot +
# Avoiding resampling with maxcell
geom_spatraster(data = map, maxcell = Inf) +
scale_fill_hypso_tint_c(
limits = map_limits,
palette = 'dem_poster',
alpha = 0.4,
labels = label_comma(),
# For the legend I use custom breaks
breaks = c(
seq(0, 500, 100),
seq(750, 1500, 250),
2000
)
) + ggplot(Frast,
add = TRUE,
alpha = .5)
base_plot
# Loading in custom fonts
myload_fonts <- function(fontname, family,
fontdir = tempdir()) {
fontname_url <- utils::URLencode(fontname)
fontzip <- tempfile(fileext = '.zip')
download.file(paste0('https://fonts.google.com/download?family=', fontname_url),
fontzip,
quiet = TRUE,
mode = 'wb'
)
unzip(fontzip,
exdir = fontdir,
junkpaths = TRUE
)
# Loading fonts
paths <- list(
regular = 'Regular.ttf',
bold = 'Bold.ttf',
italic = 'Italic.ttf',
bolditalic = 'BoldItalic.ttf'
)
namefile <- gsub(' ', '', fontname)
paths_end <- file.path(
fontdir,
paste(namefile, paths, sep = '-')
)
names(paths_end) <- names(paths)
sysfonts::font_add(family,
regular = paths_end['regular'],
bold = paths_end['bold'],
italic = paths_end['italic'],
bolditalic = paths_end['bolditalic']
)
return(invisible())
}
{r}
# Theming
showtext::showtext_auto()
# Adjusting text size
base_text_size <- 10
for1_plot <- base_plot +
# Changing guide
guides(fill = guide_legend(
title = 'm',
direction = 'horizontal',
nrow = 1,
keywidth = 1,
keyheight = 0.3,
label.position = 'bottom',
title.position = 'right',
override.aes = list(alpha = 1)
)) +
labs(
title = 'Coweeta Experimental Catchment',
subtitle = 'Factors of Safety for Landslide Scars',
caption = paste0(
'Produced by Aleksander Zarebski \n Topographic Data: Shuttle Radar Topography Mission (SRTM)\n Landslide Variables: Coweeta Experimental Catchment Observatory'
)
) +
theme_minimal(base_family = 'sans') +
theme(
plot.background = element_rect('grey97', colour = NA),
plot.caption = element_text(size = base_text_size - 6),
plot.title = element_text(face = 'bold', size = base_text_size),
plot.subtitle = element_text(
margin = margin(b = 10),
size = base_text_size - 2
),
axis.text = element_text(size = base_text_size - 4),
legend.position = 'bottom',
legend.title = element_text(size = base_text_size - 6),
legend.text = element_text(size = base_text_size - 6),
legend.key = element_rect('grey50'),
legend.spacing.x = unit(0, 'pt')
)
for1_plot
