Section 22 Overview Maps
This code is based off a script provided by Christoph Schürz. It provides overview maps of our catchment meant for reporting. They are adjusted to fit our CS.
# -------------------------------------------------------------------------
# Script for catchment overview plots and buildR setup properties ---------
# Version 0.1
# Date: 2024-02-29
# Developer: Christoph Schürz christoph.schuerz@ufz.de
# -------------------------------------------------------------------------
require(tmap)
require(terra)
require(sf)
require(rcartocolor)
require(tidyverse)
# Shows all point shapes, adapted from the ggpubr package
show_point_shapes <- function () {
d <- data.frame(p = c(0:25))
p <- ggplot(data = d, aes(x = p%%6, y = p%/%6, shape = p, label = p)) +
scale_shape_identity() +
geom_point(size = 5, fill = "red1") +
geom_text(aes(y = p%/%6 + 0.25),
size = 3) +
scale_y_reverse() +
theme_void()
return(p)
}Settings:
# Settings ----------------------------------------------------------------
## Define path to your buildR project folder
#buildr_path <- 'Y:/Gruppen/optain/3. Case Studies/1 DEU - Schwarzer Schöps/SWAT_setup/schoeps_230207'
buildr_path <- "model_data/cs10_setup/optain-cs10/"
# Define the path where the output figure and the model setup properties should be saved
#output_path <- 'Y:/Gruppen/optain/3. Case Studies/1 DEU - Schwarzer Schöps/SWAT_setup/schoeps_230207'
output_path <- "model_data/results"
dir.create(output_path, showWarnings = F)
#point_source_path <- paste0(buildr_path, '/input_data/point_source/pnts_reichenbach.shp')
point_source_path <- 'model_data/input/pointsource/cs10-ps-data.shp'# Load and prepare vector layers from the buildR project ------------------
## Basin boundary
basin <- read_sf(paste0(buildr_path, '/data/vector/basin.shp'))
## Channel layer
channel <- read_sf(paste0(buildr_path, '/data/vector/cha.shp')) %>%
mutate(type = 'channel')
## Reservoir layer
res <- read_sf(paste0(buildr_path, '/data/vector/res.shp')) %>%
mutate(type = 'reservoir')
## HRU layer (which excludes already the water areas)
## A variable 'type_group' is generated which excludes the numbers at the
## endings of names, so that only the main group name remains
## (e.g. for 'field_1234' the 'type_group' is 'field')
hru <- read_sf(paste0(buildr_path, '/data/vector/hru.shp')) %>%
mutate(type_group = str_remove(type, '_[:digit:]+$'))
## Prepare the land layer for plotting ------------------------------------
## The reservoir layer is copied as a water layer and a type_group 'watr'
## is added
water <- res %>%
select() %>%
mutate(type_group = 'watr')
## The land layer is the combination the HRUs 'type_group's and the water
## areas
land <- bind_rows(hru, water)
## Check which 'type_group's are present in your setup. You have to define
## a label and a plot color for each.
unique(land$type_group)## [1] "frst" "rngb" "urml" "a_001f" "a_002f" "utrn" "a_003f" "a_004f"
## [9] "a_005f" "a_006f" "a_007f" "a_008f" "a_009f" "a_010f" "a_011f" "a_012f"
## [17] "past" "a_013f" "a_014f" "a_015f" "a_016f" "a_017f" "a_018f" "a_019f"
## [25] "a_020f" "a_021f" "a_022f" "a_023f" "a_024f" "a_025f" "a_026f" "a_027f"
## [33] "a_028f" "a_029f" "a_030f" "a_031f" "a_032f" "a_033f" "a_034f" "a_035f"
## [41] "a_036f" "a_037f" "a_038f" "a_039f" "a_040f" "a_041f" "a_042f" "a_043f"
## [49] "a_044f" "a_045f" "a_046f" "a_047f" "a_048f" "a_049f" "a_050f" "a_051f"
## [57] "a_052f" "a_053f" "a_054f" "a_055f" "a_056f" "a_057f" "a_058f" "a_059f"
## [65] "a_060f" "a_061f" "a_062f" "a_063f" "a_064f" "a_065f" "a_066f" "a_067f"
## [73] "a_068f" "a_069f" "a_070f" "a_071f" "a_072f" "a_073f" "a_074f" "a_075f"
## [81] "a_076f" "a_077f" "a_078f" "a_079f" "a_080f" "a_081f" "a_082f" "a_083f"
## [89] "a_084f" "a_085f" "a_086f" "a_087f" "a_088f" "a_089f" "a_090f" "a_091f"
## [97] "a_092f" "a_093f" "a_094f" "a_095f" "a_096f" "a_097f" "a_098f" "a_099f"
## [105] "a_100f" "a_101f" "a_102f" "a_103f" "a_104f" "a_105f" "a_106f" "a_107f"
## [113] "a_108f" "a_109f" "a_110f" "wetf" "a_111f" "a_112f" "a_113f" "a_114f"
## [121] "a_115f" "a_116f" "a_117f" "a_118f" "a_119f" "a_120f" "a_121f" "a_122f"
## [129] "a_123f" "a_124f" "a_125f" "a_126f" "a_127f" "a_128f" "a_129f" "a_130f"
## [137] "a_131f" "a_132f" "a_133f" "a_134f" "a_135f" "a_136f" "a_137f" "a_138f"
## [145] "a_139f" "a_140f" "a_141f" "a_142f" "a_143f" "a_144f" "a_145f" "a_146f"
## [153] "a_147f" "a_148f" "a_149f" "a_150f" "a_151f" "a_152f" "a_153f" "a_154f"
## [161] "a_155f" "a_156f" "a_157f" "a_158f" "a_159f" "a_160f" "a_161f" "a_162f"
## [169] "a_163f" "a_164f" "a_165f" "a_166f" "a_167f" "a_168f" "a_169f" "a_170f"
## [177] "a_171f" "a_172f" "a_173f" "a_174f" "a_175f" "a_176f" "a_177f" "a_178f"
## [185] "a_179f" "a_180f" "a_181f" "a_182f" "a_183f" "a_184f" "a_185f" "a_186f"
## [193] "a_187f" "a_188f" "a_189f" "a_190f" "a_191f" "a_192f" "a_193f" "a_194f"
## [201] "a_195f" "a_196f" "a_197f" "a_198f" "a_199f" "a_200f" "a_201f" "a_202f"
## [209] "a_203f" "a_204f" "a_205f" "a_206f" "a_207f" "a_208f" "a_209f" "a_210f"
## [217] "a_211f" "a_212f" "a_213f" "a_214f" "a_215f" "a_216f" "a_217f" "a_218f"
## [225] "a_219f" "a_220f" "a_221f" "a_222f" "watr"Clearly, we need to change our unique fields to one “class”
agrifield_index <- grepl(x=land$type_group, pattern = "a_") %>% which()
land$type_group[agrifield_index] <- "field"
unique(land$type_group)## [1] "frst" "rngb" "urml" "field" "utrn" "past" "wetf" "watr"## Define the plot colors and the labels in the map legend for each
## 'type_group'. To have comparable maps it is advisable to use the
## same colors for the respective land uses in each case study map.
land_colors <- tribble(~type_group, ~label, ~color,
'field', 'Cropland', '#ffe18c',
'past', 'Meadows', '#33a02c',
'frst', 'Forest', '#4e6c39',
'rngb', 'Shrubs, semi-natural', '#838245',
'urml', 'Urban, moderate density', '#e8718d',
'utrn', 'Urban, transport', '#c31619',
'wetf', 'Wetlands', '#d945b7',
'watr', 'Water areas', '#5394A9'
)
## Check if each 'type_group' was defined in the 'land_colors' table./
if(!all(unique(land$type_group %in% land_colors$type_group))) {
stop("Please define all 'type_group's from the layer 'land' in 'land_colors'!")
} else {
land$type_group <- factor(land$type_group,
levels = land_colors$type_group,
labels = land_colors$label)
}
# Load and prepare raster layers from the buildR setup --------------------
# The dem and soil raster layers can be rather large and would be downsampled
# during the plotting to 10^6 pixels by the plot function. Therefore, the
# layers are aggregated to about that number of pixels before plotting to
# save time in the plot process. If this was done once the aggregated layers
# are saved in the buildR project's raster data and are loaded if
if(!file.exists(paste0(buildr_path, '/data/raster/hillshade.tif'))) {
dem <- rast(paste0(buildr_path, '/data/raster/dem.tif'))
dim_dem <- dim(dem)
fct_agg <- round(sqrt(dim_dem[1]*dim_dem[2]/1e6))
dem_agg <- aggregate(dem, fct_agg) %>%
crop(., vect(basin), mask = TRUE)
slope <- terrain(dem_agg, "slope", unit="radians")
aspect <- terrain(dem_agg, "aspect", unit="radians")
hill <- shade(slope, aspect, 40, 270)
soil <- rast(paste0(buildr_path, '/data/raster/soil.tif'))
soil <- as.factor(soil)
soil_agg <- aggregate(soil, fct_agg) %>%
crop(., vect(basin), mask = TRUE)
writeRaster(hill, paste0(buildr_path, '/data/raster/hillshade.tif'))
writeRaster(dem_agg, paste0(buildr_path, '/data/raster/dem_agg.tif'))
writeRaster(soil_agg, paste0(buildr_path, '/data/raster/soil_agg.tif'))
} else {
hill <- rast(paste0(buildr_path, '/data/raster/hillshade.tif'))
dem_agg <- rast(paste0(buildr_path, '/data/raster/dem_agg.tif'))
soil_agg <- rast(paste0(buildr_path, '/data/raster/soil_agg.tif'))
}Prepping our MetNord Dataset for this script:
metadata <- read_csv("model_data/input/met/reanalysis3_metadata.csv", show_col_types = F)
projcrs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
stations <- sf::st_as_sf(x = metadata, coords = c("Long", "Lat"), crs = projcrs)
stations$type <- "placeholder"
write_sf(stations, "model_data/input/met/stations_spatial.shp")
slr_station <- data.frame(type = "slr", Lat = 59.6606132414636, Long = 10.7823462360971)
slr_station <- sf::st_as_sf(x = slr_station, coords = c("Long", "Lat"), crs = projcrs)
write_sf(slr_station, "model_data/input/met/slr_station.shp")
mapview::mapview(stations, zcol = "Elevation")+mapview::mapview(slr_station, zcol = "type")# Load additional input layers --------------------------------------------
#
# Weather data
# The example of weather data (how they are organized in CS1) shows how to
# prepare the data for plotting in the overview maps.
#
# The workflow is always as follows:
# - Load the individual point shape layers
# - Add a column 'type' which defines which variable it is
# - bind the point layers together
# - define a point settings table, to link the 'type' with a full label which
# is printed in the figures legend, the point shape, and point colors.
#
# You can customize the code to your data for the weather and the gauge data.
#
# You can get an overview of the available point shapes with the function
# show_point_shapes() to get their corresponding values.
# Load the individual point shape layers
## Read the shape file
pcp <- read_sf( "model_data/input/met/stations_spatial.shp")
# tmp <- read_sf(paste0(weather_path, '/', 'tmp.shp')) %>%
# mutate(type = 'tmp') %>%
# select(type)
# rh <- read_sf(paste0(weather_path, '/', 'rh.shp')) %>%
# mutate(type = 'rh') %>%
# select(type)
solar <- read_sf("model_data/input/met/slr_station.shp")
# wind <- read_sf(paste0(weather_path, '/', 'wind.shp')) %>%
# mutate(type = 'wind') %>%
# select(type)
# Bind the layers together
#weather <- bind_rows(pcp, tmp, rh, solar, wind)
weather <- bind_rows(pcp, solar)
weather_point_settings <- tribble( ~type, ~label, ~shape, ~color,
'placeholder', 'Precipitation,', 22, 'azure2',
'tmp', 'Temperature,', NA, 'green4',
'rh', 'Relative humidity,', NA, 'red3',
'wind', 'Wind speed', NA, 'pink3',
'slr', 'Solar radiation', 22, 'darkgoldenrod1'
)
if(!all(unique(weather$type %in% weather_point_settings$type))) {
stop("Please define all 'type's from the layer 'weather' in 'weather_point_settings'!")
} else {
weather$type <- factor(weather$type,
levels = weather_point_settings$type,
labels = weather_point_settings$label)
}
## Gauge data
## For the gauge data please follow the same approach to customize the plot and
## the plot legend.
flow <- read_sf("model_data/input/point/cs10_outlet_utm33n.shp")
flow$type = "flow"
gauges <- bind_rows(flow)#, ncon, pcon)
gauge_point_settings <- tribble(~type, ~label, ~shape, ~color,
'flow', 'Outlet', 25, 'royalblue3',
# 'ncon', 'Nitrate', 21, 'yellow',
# 'pcon', 'Phosphorous', 23, 'brown',
)
if(!all(unique(gauges$type %in% gauge_point_settings$type))) {
stop("Please define all 'type's from the layer 'weather' in 'gauge_point_settings'!")
} else {
gauges$type <- factor(gauges$type,
levels = gauge_point_settings$type,
labels = gauge_point_settings$label)
}
# Load point shape layer of point sources
point_sources <- read_sf(point_source_path) %>%
mutate(type = 'Point source') %>%
select(type)Generating the maps:
# Plotting maps -----------------------------------------------------------
## Plot options for tmap
## It is recommended to keep most of the options as they were defined.
## It may be necessary to change the legend positions or the aspect ratio (asp)
## if the legends overlap with the maps:
## - The legend position is defined by two arguments. The first is either
## 'left' or 'right', the second is 'top' or 'bottom'.
## Capital letters move the legend even further to the plot margin.
## - The default of asp is set to 1 to have square plot panels. Values > 1
## generate wider plot panels and can give the legend more space but generates
## a smaller map.
tmap_options(title.size = 1, # Size of a), b), c), d) (relative value)
legend.text.size = 0.5, # Size of labels in legend (relative value)
legend.title.size = 0.6, # Size of legend title (relative value)
legend.title.fontface = 'bold', # Legend titles are bold
legend.position = c('RIGHT', 'bottom'), # Legend position, here right and bottom (All caps e.g. push it to the VERY right margin)
frame = F, # Plot panels are plotted without frames
asp = 1, # Aspect ratio of plot panels, now plotted squared. Change if place for legend is too narrow. value > 1 = plot wider
max.categories = Inf) # Ignore this
## Define the intervals of the scale bar, depending on the catchment size.
scale_bar_intervals <- c(0, 2.5, 5) #km
## Plot for the locations of weather stations, gauges, and point source locations.
stat_map <- tm_shape(basin) +
tm_borders() +
tm_shape(channel) +
tm_lines(col = carto_pal(7, 'Teal')[7], legend.col.show = F) +
tm_shape(res) +
tm_polygons(col = carto_pal(7, 'Teal')[4], border.col = carto_pal(7, 'Teal')[7]) +
tm_shape(point_sources) +
tm_symbols(col = 'red2', shape = 4, size = 0.05, # border.col = 'yellow2',
# palette = gauge_point_settings$color,
# shapes = gauge_point_settings$shape,
just = c('center', 'center'),
legend.col.show = F,
legend.shape.show = F) +
tm_shape(gauges) +
tm_symbols(col = 'type', shape = 'type', border.col = 'black', size = 0.15,
palette = gauge_point_settings$color,
shapes = gauge_point_settings$shape,
just = c('center', 'center'),
legend.col.show = F,
legend.shape.show = F) +
tm_shape(weather) +
tm_symbols(col = 'type', shape = 'type', size = 0.05,
palette = weather_point_settings$color,
shapes = weather_point_settings$shape,
just = c('center', 'center'),
legend.col.show = F,
legend.shape.show = F) +
tm_add_legend(title = "")+
tm_add_legend(type = 'symbol',
col = weather_point_settings$color,
shape = weather_point_settings$shape,
labels = weather_point_settings$label,
size = 0.25,
title = 'Weather') +
tm_add_legend(type = 'line',
col = carto_pal(7, 'Teal')[7],
labels = 'Channels',
size = 0.8,
title = 'Water') +
tm_add_legend(type = 'symbol',
col = gauge_point_settings$color,
shape = gauge_point_settings$shape,
labels = gauge_point_settings$label,
size = 0.4) +
tm_add_legend(type = 'fill',
col = carto_pal(7, 'Teal')[4],
border.col = carto_pal(7, 'Teal')[7],
size = 0.4,
labels = 'Reservoir, Pond') +
tm_add_legend(type = 'symbol',
col = 'red2',
# border.col = 'yellow2',
shape = 4,
size = 0.4,
labels = 'Point source',
title = 'Point sources') +
tm_layout(title = 'a)')
# Plot of the DEM and hillshade
dem_map <- tm_shape(hill) +
tm_raster(palette = gray(0:100 / 100), n = 100, legend.show = FALSE) +
tm_shape(dem_agg) +
tm_raster(alpha = 0.3, palette = terrain.colors(25),
style = 'cont', title = 'Elevation (m a.s.l)',
legend.show = T, legend.is.portrait = T, legend.reverse = T) +
tm_shape(basin) +
tm_borders() +
tm_layout(title = 'b)')
# Plot of the land use layer
land_map <- tm_shape(basin) +
tm_borders() +
tm_shape(land) +
tm_polygons(col = 'type_group', border.col = carto_pal(7, 'Teal')[7],
palette = land_colors$color, legend.show = T, lwd = 0.05,
title = 'Land use / cover') +
tm_layout(title = 'c)')
# Plot of the soil map
n_soil <- nrow(unique(soil_agg))
soil_col <- colorRampPalette(carto_pal(12, "Safe"))
soil_map <- tm_shape(soil_agg) +
tm_raster(legend.show = F, palette = soil_col(n_soil), style = 'cat') +
tm_shape(basin) +
tm_borders() +
tm_compass(size = 1.5, position = 'RIGHT') +
tm_scale_bar(breaks = scale_bar_intervals, text.size = 0.5,
position = 'RIGHT') +
tm_layout(title = 'd)')
# Combine all plot panels to the overview map
overview_map <- tmap_arrange(stat_map, dem_map, land_map, soil_map, ncol = 2)
# Save the plot
tmap_save(overview_map, paste0(output_path, '/overview_map.png'), width = 9, height = 7)## stars object downsampled to 797 by 1254 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 797 by 1254 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 797 by 1254 cells. See tm_shape manual (argument raster.downsample)## Map saved to D:\git\swat-cs10\model_data\results\overview_map.png## Resolution: 2700 by 2100 pixels## Size: 9 by 7 inches (300 dpi)## stars object downsampled to 797 by 1254 cells. See tm_shape manual (argument raster.downsample)## stars object downsampled to 797 by 1254 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 797 by 1254 cells. See tm_shape manual (argument raster.downsample)