Adapating code from https://github.com/milos-agathon/precipitation-maps to compare spain precipitation in 2023 and previous years

pluja.R

## Preparations

pkgs <- c("tidyverse", "pRecipe", "giscoR", "terra", "rayshader", "patchwork")
rlang::check_installed(pkgs)
suppressPackageStartupMessages({
    invisible(lapply(pkgs, library, character.only = TRUE, quietly = TRUE))
})

## Data

pRecipe::download_data(dataset = "mswep", path = getwd(), domain = "raw", timestep = "yearly")
list.files()

# CH
# country_sf <- giscoR::gisco_get_countries(country = "CH", resolution = "1")
# mswep_data <- terra::rast("mswep_tp_mm_global_197902_202301_025_yearly.nc") %>% 
#     terra::crop(country_sf)
# terra::plot(mswep_data[[1]])
# plot(sf::st_geometry(country_sf), add = TRUE)

# ES
country_sf <- giscoR::gisco_get_countries(country = "ES", resolution = "1")
mswep_data <- terra::rast("mswep_tp_mm_global_197902_202301_025_yearly.nc") %>%
    terra::crop(country_sf)
terra::plot(mswep_data[[1]])
plot(sf::st_geometry(country_sf), add = TRUE)

names(mswep_data) <- 1979:2023
mswep_df <- mswep_data %>%
    as.data.frame(xy = TRUE) %>% 
    pivot_longer(!c("x", "y"), names_to = "year", values_to = "precipitation") %>% 
    filter(year != 2023)

mswep_last_df <- mswep_data %>%
    as.data.frame(xy = TRUE) %>% 
    pivot_longer(!c("x", "y"), names_to = "year", values_to = "precipitation") %>% 
    filter(year == 2023)

# Ploting the map

theme_for_the_win <- function(){
    theme_minimal() +
        theme(
            axis.line = element_blank(),
            axis.title.x = element_blank(),
            axis.title.y = element_blank(),
            axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            legend.position = "right",
            legend.title = element_text(
                size = 11, color = "grey10"
            ),
            legend.text = element_text(
                size = 10, color = "grey10"
            ),
            panel.grid.major = element_line(
                color = NA
            ),
            panel.grid.minor = element_line(
                color = NA
            ),
            plot.background = element_rect(
                fill = NA, color = NA
            ),
            legend.background = element_rect(
                fill = "white", color = NA
            ),
            panel.border = element_rect(
                fill = NA, color = NA
            ),
            plot.margin = unit(
                c(
                    t = 0, r = 0,
                    b = 0, l = 0
                ), "lines"
            )
        )
}

breaks <- classInt::classIntervals(
    mswep_last_df$precipitation,
    n = 5,
    style = "equal"
)$brks

colors <- hcl.colors(
    n = length(breaks),
    palette = "Temps",
    rev = TRUE
)

map1 <- ggplot(
    data = mswep_last_df
) +
    geom_raster(
        aes(
            x = x,
            y = y,
            fill = precipitation
        )
    ) +
    geom_contour(
        aes(
            x = x,
            y = y,
            z = precipitation 
        ), color = "white" # add this line
    ) +
    geom_sf(
        data = country_sf,
        fill = "transparent",
        color = "grey10",
        size = .5
    ) +
    scale_fill_gradientn(
        name = "mm",
        colors = colors,
        breaks = breaks,
        labels = round(breaks, 0), # use round(breaks, 0)
        limits = c(
            min(mswep_last_df$precipitation),
            max(mswep_last_df$precipitation)+50
        )
    ) +
    # facet_wrap(~year) +
    guides(
        fill = guide_colourbar(
            direction = "vertical",
            barheight = unit(50, "mm"),
            barwidth = unit(5, "mm"),
            title.position = "top",
            label.position = "right",
            title.hjust = .5,
            label.hjust = .5,
            ncol = 1,
            byrow = FALSE
        )
    ) +
    theme_for_the_win() +
    labs(title = "Precipitació total de 2023-01 a 2023-07")

map1

# AVERAGE PRECIPATION

mswep_average_df <- mswep_df |>
    dplyr::group_by(
        x, y, .drop = FALSE
    ) |>
    dplyr::summarise(
        mean = mean(precipitation)
    )

head(mswep_average_df)

breaks <- classInt::classIntervals(
    mswep_average_df$mean,
    n = 5,
    style = "equal"
)$brks

colors <- hcl.colors(
    n = length(breaks),
    palette = "Temps",
    rev = TRUE
)

map2 <- ggplot(
    data = mswep_average_df
) +
    geom_raster(
        aes(
            x = x,
            y = y,
            fill = mean
        )
    ) +
    geom_contour(
        aes(
            x = x,
            y = y,
            z = mean 
        ), color = "white" # add this line
    ) +
    geom_sf(
        data = country_sf,
        fill = "transparent",
        color = "grey10",
        size = .5
    ) +
    scale_fill_gradientn(
        name = "mm",
        colors = colors,
        breaks = breaks, 
        labels = round(breaks, 0), # use round(breaks, 0)
        limits = c(
            min(mswep_average_df$mean),
            max(mswep_average_df$mean)
        )
    ) +
    guides(
        fill = guide_colourbar(
            direction = "vertical",
            barheight = unit(50, "mm"),
            barwidth = unit(5, "mm"),
            title.position = "top",
            label.position = "right",
            title.hjust = .5,
            label.hjust = .5,
            ncol = 1,
            byrow = FALSE
        )
    ) +
    theme_for_the_win() +
    labs(title = "Precipitació anual mitjana de 1979 a 2022")

map2

(map2 + map1) +
    labs(caption = "Autor @jrosell | Font: mswep")