Merge branch 'main' of github.com:geco-bern/les

analysis/box_modelling.R

@@ -333,3 +333,18 @@ df_dist |>
   labs(x = "Time", y = "C pool")
 
 ggsave(here("book/images/disturbance_recovery.png"), width = 4, height = 3)
+
+# for exam
+ggplot() +
+  geom_line(aes(year, co2),
+            data = df_const |> 
+              filter(year >= 1850),
+            color = "black"
+  ) +
+  geom_vline(xintercept = 2023, linetype = "dotted") +
+  theme_classic() +
+  labs(x = "Year", 
+       y = expression(paste("CO"[2], " (ppm)")),
+       subtitle = expression(paste("Atmospheric CO"[2])))
+
+ggsave(here::here("fig/co2_stabilisation.png"), width = 4, height = 3)

analysis/plot_map_aridityindex.R

@@ -0,0 +1,145 @@
+library(sf)
+library(terra)
+library(dplyr)
+library(ggplot2)
+library(readr)
+library(rgeco)
+library(sf)
+library(khroma)
+
+source(here::here("R/plot_discrete_cbar.R"))
+
+## Global map of AI ----------------
+# breaks for color bar
+breaks <- c(0, 0.03, 0.1, 0.2, 0.3, 0.4, 0.5, 0.65, 0.8, 1, 1.2, 1.6, 2, 2.5, 3, Inf)
+# breaks <- c(0, 0.03, 0.2, 0.5, 0.65, Inf)
+
+# raster_ai <- rast("~/data/greve/ep_over_p_cru_ncep.nc")
+
+# data from https://doi.org/10.6084/m9.figshare.7504448.v4
+rasta <- rast("~/data/aridityindex_zomer_2022/Global-AI_ET0_v3_annual/ai_v3_yr.tif")
+
+# values are provided as integers, multiplied by 1e4
+values(rasta) <- values(rasta) * 1e-4
+
+# aggregate to 0.1 deg
+rasta_agg <- aggregate(
+  rasta,
+  fact = res(rasta)[1]^-1*0.1,
+  fun = "mean"
+)
+
+# bin values to get a discrete color scale (personal preference)
+rasta_bin <- rasta_agg
+values(rasta_bin) <- cut(
+  values(rasta_agg),
+  breaks = breaks, 
+  right = FALSE
+)
+
+# get coast outlines
+layer_coast <- rnaturalearth::ne_coastline(
+  scale = 110, 
+  returnclass = "sf"
+  )
+
+# get ocean layer
+layer_ocean <- rnaturalearth::ne_download(
+  scale = 110,
+  type = "ocean",
+  category = "physical",
+  returnclass = "sf",
+  destdir = here::here("data/")
+)
+
+# construct map
+ggmap <- ggplot() +
+
+  # aridity index raster
+  tidyterra::geom_spatraster(
+    data = rasta_bin,
+    show.legend = FALSE
+    ) +
+
+  # coastline
+  geom_sf(
+    data = coast,
+    colour = 'black',
+    linewidth = 0.1
+  ) +
+
+  # ocean to mask zeros
+  geom_sf(
+    data = layer_ocean,
+    color = NA,
+    fill = "azure3"
+    ) +
+
+  # color palette from the khroma package
+  scale_fill_roma(discrete = TRUE, name = "")  +
+  coord_sf(
+    ylim = c(-60, 85),
+    expand = FALSE   # to draw map strictly bounded by the specified extent
+  ) +
+  xlab('') +
+  ylab('') +
+  theme_bw() +
+  theme(axis.ticks.y.right = element_line(),
+        axis.ticks.x.top = element_line(),
+        panel.grid = element_blank(),
+        plot.background = element_rect(fill = "white")
+  )
+
+gglegend <- plot_discrete_cbar(
+  breaks = breaks,
+  colors = c(khroma::color("roma")(length(breaks)-1)),
+  legend_title = "",
+  legend_direction = "vertical",
+  width = 0.03,
+  font_size = 3,
+  expand_size_y = 0.5,
+  spacing = "constant"
+)
+
+cowplot::plot_grid(ggmap, gglegend, ncol = 2, rel_widths = c(1, 0.10))
+
+ggsave(here::here("book/images/map_aridityindex.png"), width = 12*0.8, height = 7*0.8)
+
+## "Donohue" plot --------------
+
+# get fAPAR at 0.5 deg
+rasta_fpar <- terra::rast("~/data/MODIS-C006_MOD15A2_LAI_FPAR_zmaw/MODIS-C006_MOD15A2__LAI_FPAR__LPDAAC__GLOBAL_0.5degree__UHAM-ICDC__2000_2018__MON__fv0.02_ANNMAX_MEAN.nc")["fpar"]$fpar
+
+# resample AI data to fAPAR grid
+rasta_ai_agg <- resample(rasta_agg, rasta_fpar, method = "bilinear")
+
+df <- as.data.frame(rasta_fpar, xy = TRUE) |> 
+  left_join(
+    as.data.frame(rasta_ai_agg, xy = TRUE),
+    by = c("x", "y")
+  ) |> 
+  as_tibble()
+
+df |> 
+  ggplot(aes(awi_pm_sr_yr, fpar)) +
+  geom_hex(
+    bins = 50,
+    # show.legend = FALSE
+    ) +
+  xlim(0, 1.5) +
+  theme_classic() +
+  scale_fill_gradientn(
+    colours = colorRampPalette( c("gray65", "navy", "red", "yellow"))(5),
+    # trans = "log",
+    ) +
+  labs(
+    x = expression(paste(italic(P)/PET)),
+    y = "fAPAR"
+  )
+
+ggsave(here::here("book/images/donohue.png"), width = 6, height = 4)
+
+
+
+
+

analysis/plot_radiation_sites_biomes.R

@@ -305,3 +305,87 @@ ddf_sub |>
   ) +
   facet_wrap(~sitename)
 
+## Extra plot --------------------
+ddf_sub |> 
+  filter(sitename == "DE-Hai") |> 
+  pivot_longer(cols = c(netrad, heat, le)) |>  # , res
+  ggplot(aes(date, value, color = name)) +
+  geom_hline(yintercept = 0, linetype = "dotted") +
+  geom_line() +
+  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
+  scale_color_manual(
+    name = "",
+    labels = c(
+      heat = expression(paste(italic("H"))),
+      le = expression(paste(lambda, italic("E"))),
+      netrad = expression(paste(italic("R")[net]))
+      # res = "Residual"
+    ),
+    values = c(
+      heat = "#E69F00",
+      le = "#56B4E9",
+      netrad = "black"
+      # res = "grey70"
+    )
+  ) +
+  theme_classic() +
+  labs(
+    x = "",
+    y = expression(paste("Energy flux (W m"^-2, ")"))
+  )
+
+ggsave(
+  here::here("fig/seasonal_cycle_de-hai.png"),
+  width = 5, 
+  height = 3
+)
+
+hdf_sub |> 
+  filter(sitename == "DE-Hai") |> 
+  pivot_longer(cols = c(netrad, heat, le)) |>  # , res
+  ggplot(aes(hm, value, color = name)) +
+  geom_hline(yintercept = 0, linetype = "dotted") +
+  geom_line() +
+  scale_x_time() +
+  # scale_x_datetime(date_breaks= "2 hours", date_labels = "%H:%M") +
+  scale_color_manual(
+    name = "",
+    labels = c(
+      heat = expression(paste(italic("H"))),
+      le = expression(paste(lambda, italic("E"))),
+      netrad = expression(paste(italic("R")[net]))
+      # res = "Residual"
+    ),
+    values = c(
+      heat = "#E69F00",
+      le = "#56B4E9",
+      netrad = "black"
+      # res = "grey70"
+    )
+  ) +
+  # scale_linetype_manual(
+  #   name = "",
+  #   labels = c(
+  #     heat = expression(paste(italic("H"))),
+  #     le = expression(paste(lambda, italic("E"))),
+  #     netrad = expression(paste(italic("R")[net])),
+  #     res = "Residual"
+  #   ),
+  #   values = c(
+  #     heat = "solid",
+  #     le = "solid",
+  #     netrad = "solid",
+  #     res = "dashed"
+  #   )
+  # ) +
+  theme_classic() +
+  labs(
+    x = "Time of day",
+    y = expression(paste("Energy flux (W m"^-2, ")"))
+  )
+
+ggsave(
+  here::here("fig/diurnal_cycle_de-hai.png"),
+  width = 6, 
+  height = 3
+)

book/ecohydrology.qmd

@@ -270,6 +270,24 @@ gg1 <- df_swc |>
   xlim(0, 0.5) + 
   theme(legend.position="none")
 
+# # extra plot
+# df_swc |> 
+#   filter(soilclass %in% c("Sand", "Loam", "Clay")) |> 
+#   ggplot(aes(swc, psi, color = soilclass)) +
+#   geom_line() +
+#   geom_hline(yintercept = c(-1000, -150000), linetype = "dotted") +
+#   khroma::scale_color_okabeito() +
+#   labs(
+#     x = "Volumetric soil water content",
+#     y = expression(paste("Soil matric potential (mm)"))
+#   ) +
+#   theme_classic() +
+#   ylim(-0.5e6, 0) + 
+#   xlim(0, 0.5) + 
+#   theme(legend.position="none")
+# 
+# ggsave(here::here("fig/soil_matric_potential_vs_swc.png"), width = 4, height = 3)
+
 gg2 <- df_swc |> 
   ggplot(aes(swc, -psi, color = soilclass)) +
   geom_line() +
@@ -1016,11 +1034,21 @@ Related patterns emerge from the annual sums of ecosystem water fluxes (AET, PET
 knitr::include_graphics("images/budyko_fluxnet.png")
 ```
 
-<!-- - Start by: constant EF under un-limited conditions, decline as water becomes limiting -->
-
-<!-- ### (Donohue graph) -->
-
-<!-- ### Budyko {#sec-budyko} -->
+<!-- XXX Add: -->
+<!-- Budyko framework -->
+
+<!-- - Chatchment characteristics -->
+<!-- - Soil -->
+<!-- - Vegetation -->
+<!-- - Topography -->
+<!-- - Groundwater -->
+<!-- - Hydroclimate -->
+<!-- - Precipitation magnitude -->
+<!-- - Precipitation distribution -->
+<!-- - Seasonality of precipitation
+and PET -->
+
+<!-- - Fig. 3 in Caracciolo et al. 2018 https://doi.org/10.1007/s11269-018-1984-7 -->
 
 ## Landscape-scale heterogeneity
 

book/globalcarbonpatterns.qmd

@@ -29,6 +29,8 @@ $$
 I_\mathrm{0} = (1-\alpha_\mathrm{sw})\; \tau(z, f_c) \; I_\mathrm{TOA}
 $$ {#eq-sin}
 
+<!-- XXX give values and units of all variables -->
+
 ```{r echo=FALSE}
 #| label: fig-solar_radiation
 #| fig-cap: "Map of solar radiation incident at the Earth surface, averaged over 1981-2010, based on @chelsa_bioclim. Figure from [https://en.wikipedia.org/wiki/Solar_irradiance#/media/File:RSDS_wiki.png](https://en.wikipedia.org/wiki/Solar_irradiance#/media/File:RSDS_wiki.png)."
@@ -47,6 +49,8 @@ I_\mathrm{TOA} = I_S r_E^{-2} \cos \theta_z
 $$ {#eq-s-toa}
 $\theta_z$ is the solar zenith angle. The term $\cos \theta_z$ accounts for the dependence of the solar radiation on the angle at which the sun's rays reach the Earth surface (no terrain considered). It varies cyclically over the course of a day (with *hour-of-day*) and over the course of a year (with *day-of-year*). The zenith angle is zero when the sun is directly above the observer - in the zenith (@fig-solar_geometry). At this point, the intensity of the solar radiation is at its maximum.
 
+<!-- XXX explain why ~1/r^2 -->
+
 Note that $r_E$ is not constant over the course of a year as the Earth rotates around the sun not following a perfect circle but an ellipse.
 
 ```{r echo=FALSE}
@@ -489,6 +493,7 @@ knitr::include_graphics("images/seasonal_cycle.png")
 7. Consider mature forest ecosystems (*sensu* @odum69sci, @fig-odum) in a moist tropical, temperate and boreal region. Where do you expect the largest annual NEP? Where do you expect the largest daily NEP during the peak growing season?
 :::
 
+<!-- XXX remove exercise 5: not trivial explanation -->
 
 #### The breathing of the Earth