# Data Visualisation in R / Adatábrázolás R-ben
# Presentation of research sites on an overview map / Kutatási helyszínek bemutatása átnézeti térképen
# Ákos Bede-Fazekas, 13 Nov 2025

# Please change the path below to the path of the unzipped "files" folder 
setwd("s:/Oktatás, tanulás/Oktatás, tanulás 2025_26_1_ősz/Adatábrázolás R-ben/20251113 Bede-Fazekas Ákos - térképkészítés/előkészítés/files")

# Install (if needed) and load packages
needed_packages <- c("sf", "raster", "prettymapr", "rworldmap")
invisible(lapply(X = needed_packages, FUN = \(package) {
	if (!require(package = package, character.only = TRUE)) {
		install.packages(pkgs = package, verbose = FALSE, quiet = TRUE, INSTALL_opts = c("--no-lock"), repos = "http://cran.us.r-project.org")
		library(package = package, character.only = TRUE)
	}
}))

# Load the research sites
sites <- read.csv(file = "input/sites.csv", sep = ";", dec = ".", encoding = "UTF-8")
print(sites)
sites <- na.omit(sites)

# First attempt
plot(x = sites$longitude, y = sites$latitude)
# But what about lines and polygons?...

# This is where package "sf" is needed
# Convert to a spatial dataset ("simple features")
sites <- st_as_sf(
	x = sites,
	coords = c("longitude", "latitude"), # always the west-east coordinate first, the south-north coordinate second
	crs = 4326 # WGS-84
)
class(sites) # Simple Features (sf) and data.frame -> can be used as a data.frame

print(sites)

# Plot colored by a column
plot(sites[, "year"])

# Plot only the geometry (i.e., the location)
plot(st_geometry(sites))

# With axes
plot(st_geometry(sites), axes = TRUE)

# And graticule (grid lines at pretty coordinates)
plot(st_geometry(sites), axes = TRUE, graticule = TRUE)

# We want to add the year as labels right to the points
# Coordinates are needed
coordinates <- st_coordinates(sites)
head(coordinates)

# Now we can add the labels to the specified locations (with some rightward shift)
plot(st_geometry(sites), axes = TRUE, graticule = TRUE)
text(labels = as.character(sites$year), x = coordinates[, 1] + 0.25, y = coordinates[, 2])

# There are sites  visited in multiple years -> let's keep one record per locations ("name" column)
sites$name
sites_unique <- sites[!duplicated(sites$name), ]
coordinates <- st_coordinates(sites_unique)
plot(st_geometry(sites_unique), axes = TRUE, graticule = TRUE)
text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.25, y = coordinates[, 2])

# Let's save the plot to a PNG file
file_name <- "map.png"
png(file = file_name, width = 1400, height = 1000)
plot(st_geometry(sites_unique), axes = TRUE, graticule = TRUE)
text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.25, y = coordinates[, 2])
dev.off()

# Different canvas size -> point & text sizes, horizontal shift, margins and axis label size should be adjusted
png(file = file_name, width = 1400, height = 1000)
par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
plot(st_geometry(sites_unique), pch = 16, cex = 1.5, cex.axis = 1.5, axes = TRUE, graticule = TRUE)
text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5)
dev.off()

# I hate opening the the PNG file every time, so I will force the OS to open it after created
source("open_map.r")
open_map

# Map drawing involves many potential errors, but whether there are errors or not, we must close the PNG device to see the result.
# Let's use error-proof drawing!
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = { # R tries to run the following commands
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5)
	},
	finally = { # whether it failed or not, it will close the device and open the file
		dev.off()
		open_map(file_name)
	}
)

# Let's add some layers (country border, hidrology etc.) to help the orientation of the map reader
# Load the files
load("input/hungary.RData")
load("input/lakes.RData")
load("input/rivers.RData")

# Let's visualize
plot(hungary)

# The dataset "hungary" contains one line ("LINESTRING")
print(hungary)

# While dataset "rivers" contains several lines
print(rivers)

# The datasets "lakes" contains several polygons ("POLYGON") for a single entry ("MULTIPOLYGON")
print(lakes)

# Let's visualize them together
plot(hungary)

# Draw the three datasets as three layers on the same map
plot(hungary)
plot(lakes, add = TRUE)
plot(rivers, add = TRUE)

# Add the three layers to the map of our sites
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5)
		plot(hungary, add = TRUE)
		plot(lakes, add = TRUE)
		plot(rivers, add = TRUE)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# Where are the new layers??
# The coordinates are super strange (not latitudes and longitudes)
plot(hungary, axes = TRUE)
head(st_coordinates(hungary))

# The coordinate reference system (CRS, projection) of the layers differ
st_crs(sites)$epsg # 4326 = WGS-84 in degrees
st_crs(hungary)$epsg # 23700 = Egységes Országos Vetület (EOV) in meters

# Transform (reproject) the datasets to WGS-84
hungary <- st_transform(x = hungary, crs = 4326)
lakes <- st_transform(x = lakes, crs = 4326)
rivers <- st_transform(x = rivers, crs = 4326)

# Now try again
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5)
		plot(hungary, add = TRUE)
		plot(lakes, add = TRUE)
		plot(rivers, add = TRUE)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# Always the first layer defines the extent of the map
# If we want to show the entire country, we should
#    (1) change the drawing order, or
#    (2) use xlim and ylim parameters of the first plot()
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5)
		plot(lakes, add = TRUE)
		plot(rivers, add = TRUE)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# Load a digital elevation model (DEM) to visualize the topography
dem <- raster("input/dem.tif")
class(dem) # it is a single-layer raster ("RasterLayer")

# How does it look like?
plot(dem)

# Use a more common color palette
colors <- colorRampPalette(c("#336600", "#F3CA89", "#D9A627", "#C2B0B0", "#E5E5E5"))(255)
plot(dem, col = colors)

# When adding to an existing plot as a new layer, the legend should be removed
plot(dem, legend = FALSE, col = colors)

# Add the raster and then the border of Hungary once again
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		plot(dem, legend = FALSE, col = colors, add = TRUE)
		plot(hungary, add = TRUE) # this was covered by the raster, so we draw it again (i.e., above the raster)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5)
		plot(lakes, add = TRUE)
		plot(rivers, add = TRUE)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# Let's prettify it
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, lwd = 5, col = "gray35", cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		plot(dem, legend = FALSE, col = colors, add = TRUE)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, col = "purple4", add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5, col = "purple4")
		plot(hungary, lwd = 5, col = "gray35", add = TRUE)
		plot(lakes, col = "darkblue", border = "transparent", add = TRUE)
		plot(rivers, lwd = 2, col = "darkblue", add = TRUE)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# Add north arrow and scale bar to the bottom right corner and a thick black border to the map
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, lwd = 5, col = "gray35", cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		plot(dem, legend = FALSE, col = colors, add = TRUE)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, col = "purple4", add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5, col = "purple4")
		plot(hungary, lwd = 5, col = "gray35", add = TRUE)
		plot(lakes, col = "darkblue", border = "transparent", add = TRUE)
		plot(rivers, lwd = 2, col = "darkblue", add = TRUE)
		suppressMessages(addscalebar(plotunit = "km", plotepsg = 4326, pos = "bottomright", padin = c(0.5, 0.5), htin = 0.2, labelpadin = 0.2, label.cex = 1.7, label.col = "black"))
		addnortharrow(pos = "bottomright", scale = 1.2, padin = c(0.5, 1.0))
		box(lwd = 3) # thick black border
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# We are now planning to add an inset figure showing Europe and the position of Hungary therein.
# Let's load the world's countries (from package "rworldmap") and keep the geometry only
countries <- st_geometry(st_as_sf(getMap(resolution = "low")))
plot(countries)

# Load the main rivers from an ESRI shapefile (shp)
main_rivers <- st_read(dsn = "input/main_rivers.shp", quiet = TRUE)
str(main_rivers)

# Remove the unnecessary data columns
main_rivers <- st_geometry(main_rivers)
plot(main_rivers)

# Add the inset figure
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, lwd = 5, col = "gray35", cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		plot(dem, legend = FALSE, col = colors, add = TRUE)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, col = "purple4", add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5, col = "purple4")
		plot(hungary, lwd = 5, col = "gray35", add = TRUE)
		plot(lakes, col = "darkblue", border = "transparent", add = TRUE)
		plot(rivers, lwd = 2, col = "darkblue", add = TRUE)
		suppressMessages(addscalebar(plotunit = "km", plotepsg = 4326, pos = "bottomright", padin = c(0.5, 0.5), htin = 0.2, labelpadin = 0.2, label.cex = 1.7, label.col = "black"))
		addnortharrow(pos = "bottomright", scale = 1.2, padin = c(0.5, 1.0))
		box(lwd = 3)
		par(fig = c(0.06, 0.4, 0.77, 0.955), new = TRUE, mar = c(0, 0, 0, 0)) # top left corner, without margin
		plot(countries, lwd = 1.5, border = "black", col = "white", axes = FALSE, xlab = "", ylab = "", xlim = c(5, 40), ylim = c(35, 52))
		plot(main_rivers, lwd = 1.8, col = "gray75", add = TRUE)
		box(lwd = 2) # thick black border
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# Add gray background (i.e., a gray rectangle larger than the map's extent)
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, lwd = 5, col = "gray35", cex.axis = 1.5, axes = TRUE, graticule = TRUE)
		plot(dem, legend = FALSE, col = colors, add = TRUE)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, col = "purple4", add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5, col = "purple4")
		plot(hungary, lwd = 5, col = "gray35", add = TRUE)
		plot(lakes, col = "darkblue", border = "transparent", add = TRUE)
		plot(rivers, lwd = 2, col = "darkblue", add = TRUE)
		suppressMessages(addscalebar(plotunit = "km", plotepsg = 4326, pos = "bottomright", padin = c(0.5, 0.5), htin = 0.2, labelpadin = 0.2, label.cex = 1.7, label.col = "black"))
		addnortharrow(pos = "bottomright", scale = 1.2, padin = c(0.5, 1.0))
		box(lwd = 3)
		par(fig = c(0.06, 0.4, 0.77, 0.955), new = TRUE, mar = c(0, 0, 0, 0))
		plot(countries, lwd = 1.5, border = "black", col = "white", axes = FALSE, xlab = "", ylab = "", xlim = c(5, 40), ylim = c(35, 52))
		rect(xleft = -10, ybottom = 30, xright = 60, ytop = 60, col = "gray90") # rectangle filled with gray color
		plot(countries, lwd = 1.5, border = "black", col = "white", add = TRUE) # draw the countries again (above the gray rectangle)
		plot(main_rivers, lwd = 1.8, col = "gray75", add = TRUE)
		box(lwd = 2)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

# We want to show the extent of the main map at the inset figure (as a red rectangle)
# Let's get the bounding box of Hungary
bounding_box <- st_bbox(hungary)
print(bounding_box)

# Add some extra space along the edges
bounding_box <- bounding_box + c(-0.1, -0.1, 0.1, 0.1)

# Use this as the xlim and ylim parameters of the main map and add a red rectangle to the inset map
# (This could be done with rect(), now I use extent())
file_name <- "map.png"
png(file = file_name, width = 1400, height = 1000)
tryCatch(
	expr = {
		par(mar = c(3, 4, 1, 1), mgp = c(1, 1.5, 0), tck = -0.01)
		plot(hungary, lwd = 5, col = "gray35", axes = TRUE, graticule = TRUE, cex.axis = 1.5, xlim = bounding_box[c("xmin", "xmax")], ylim = bounding_box[c("ymin", "ymax")]) # xlim and ylim are specified
		plot(dem, legend = FALSE, col = colors, add = TRUE)
		plot(st_geometry(sites_unique), pch = 16, cex = 1.5, col = "purple4", add = TRUE)
		text(labels = as.character(sites_unique$year), x = coordinates[, 1] + 0.15, y = coordinates[, 2], cex = 1.5, col = "purple4")
		plot(hungary, lwd = 5, col = "gray35", add = TRUE)
		plot(lakes, col = "darkblue", border = "transparent", add = TRUE)
		plot(rivers, lwd = 2, col = "darkblue", add = TRUE)
		suppressMessages(addscalebar(plotunit = "km", plotepsg = 4326, pos = "bottomright", padin = c(0.5, 0.5), htin = 0.2, labelpadin = 0.2, label.cex = 1.7, label.col = "black"))
		addnortharrow(pos = "bottomright", scale = 1.2, padin = c(0.5, 1.0))
		box(lwd = 3)
		par(fig = c(0.06, 0.4, 0.77, 0.955), new = TRUE, mar = c(0, 0, 0, 0))
		plot(countries, lwd = 1.5, border = "black", col = "white", axes = FALSE, xlab = "", ylab = "", xlim = c(5, 40), ylim = c(35, 52))
		rect(xleft = -10, ybottom = 30, xright = 60, ytop = 60, col = "gray90")
		plot(countries, lwd = 1.5, border = "black", col = "white", add = TRUE)
		plot(main_rivers, lwd = 1.8, col = "gray75", add = TRUE)
		plot(extent(as.numeric(bounding_box[c("xmin", "xmax", "ymin", "ymax")])), col = "red", lwd = 3, add = TRUE) # red rectangle
		box(lwd = 2)
	},
	finally = {
		dev.off()
		open_map(file_name)
	}
)

