-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
143 lines (99 loc) · 5.89 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# aemetools <a href="https://limnotrack.github.io/aemetools/"><img src="man/figures/logo.png" alt="aemetools website" align="right" height="120"/></a>
<!-- badges: start -->
[](https://lifecycle.r-lib.org/articles/stages.html#stable) [](https://github.com/limnotrack/aemetools/actions/workflows/R-CMD-check.yaml) [](https://app.codecov.io/gh/limnotrack/aemetools?branch=main)
<!-- badges: end -->
aemetools is designed to work with [AEME](https://github.com/limnotrack/AEME/). It contains a range of functions to assist in setting up simulations for a lake site.
## Development
This package was developed by [LimnoTrack](http://limnotrack.com/) as part of the Lake Ecosystem Restoration New Zealand Modelling Platform (LERNZmp) project. <a href="http://limnotrack.com/"><img src="man/figures/limnotrack_border.jpg" alt="LimnoTrack website" align="right" height="80"/></a>
## Overview
Currently, this package can be used to:
- Download meteorological data from [ERA5-Land](https://www.ecmwf.int/en/era5-land) for any point in New Zealand from 1999-2022, or download ERA5-Land netCDF files for any area in the world.
- Set up and run hydrological simulations using the suite of models from the [`airGR`](https://hydrogr.github.io/airGR/) package using catchment, reach and lake data.
- Conduct a sensitivity analysis on the parameters for the AEME models.
- Calibrate the AEME models using lake observational data.
```{r include=FALSE, echo=FALSE}
data("era5_ref_table", package = "aemetools")
era5_names <- paste0(tolower(era5_ref_table$variable[-1]), collapse = ", ")
```
## Installation
You can install the development version of aemetools from [GitHub](https://github.com/) with:
``` r
# install.packages("devtools")
devtools::install_github("limnotrack/aemetools")
```
```{r setup, warning=FALSE, message=FALSE}
library(aemetools)
```
### Download NZ point meteorological data
Currently, there is ERA5-Land data (\~9km grid spacing) archived for New Zealand (166.5/-46.6/178.6/-34.5) for the time period 1980-2023 with the main meteorological variables (`r era5_names`) required to drive hydrological and hydrodynamic models. This can be easily downloaded using the example below. There is a `parallel` switch which allows you to use multiple cores on your computer to speed up the download.
```{r download-era5}
lon <- 176.2717
lat <- -38.079
variables <- c("MET_tmpair", "MET_pprain")
met <- get_era5_point(lat = lat, lon = lon, years = 2020:2021,
variables = variables)
summary(met)
```
```{r plot-era5}
library(ggplot2)
library(tidyr)
met |>
pivot_longer(cols = c(MET_tmpair, MET_pprain)) |>
ggplot(aes(x = Date, y = value)) +
geom_line() +
facet_wrap(~name, scales = "free_y", ncol = 1) +
theme_bw()
```
### Calibrate AEME model
See the vignette [here](https://limnotrack.github.io/aemetools/articles/calibrate-aeme.html).
### Sensitivity analysis for AEME models
See the vignette [here](https://limnotrack.github.io/aemetools/articles/sensitivity-analysis.html).
### Hydrological modelling - Run GR4J model
Here is simple example set up for one of the lake inflows into Lake Rotorua. First, the input for the model are generated using the stream ID (nzsegment), and spatial features (sf objects) of the reaches, lake and catchment (including sub-catchments), observed discharge (if available) meteorological data (air temperature and precipitation). It recursively creates an upstream network using the nzsegment, then combines the subcatchments of all the upstream reaches (`sf::st_union()`) to calculate the area of the catchment.
```{r make-gr4j-inputs, results='asis', fig.keep='all', warning=FALSE}
data_dir <- system.file("extdata/hydro/", package = "aemetools")
lake <- readRDS(file.path(data_dir, "lake.rds"))
reaches <- readRDS(file.path(data_dir, "reaches.rds"))
catchments <- readRDS(file.path(data_dir, "catchments.rds"))
met <- readRDS(file.path(data_dir, "met.rds"))
obs_flow <- readRDS(file.path(data_dir, "obs_flow.rds"))
FUN_MOD <- airGR::RunModel_GR4J
id <- 4087861 # nzsegment
inputs <- make_GR_inputs(id = id, reaches = reaches, lake = lake,
catchments = catchments, obs_flow = obs_flow, met = met,
lat = lat, FUN_MOD = FUN_MOD,
plot = TRUE)
```
Within the `airGR` package, there are calibration algorithms which allows you to calibrate the hydrological model if discharge data for the reach is available. The calibrated parameters can be passed to the `run_GR` function to run the selected model.
```{r run-gr4j, warning=FALSE}
#' airGR uses indices to run the model, so first we split our observed data in
#' half (0.5) for calibration and validation periods based on when the
#' observation data starts (which is provided in `inputs$data$start`).
idx_spl <- floor(nrow(inputs$data[inputs$start:nrow(inputs$data), ])
* 0.5)
#' Use a model warmup period as everything before when the observations start.
warmup <- 1:(inputs$start - 1)
# Set the indices for the calibration period
cal_idx <- inputs$start:(idx_spl + inputs$start)
# Run the calibration and assign the output
calib <- calib_GR(inputs = inputs, warmup = warmup, run_index = cal_idx)
# Extract the calibrated parameters
param <- calib$ParamFinalR
# Run the model
output <- run_GR(inputs = inputs, param = param,
warmup = warmup, run_index = cal_idx)
# Plot the output
plot(output, Qobs = inputs$data$Qmm[cal_idx])
```