In hanel/hmkzLAPV: Vyuka HMKZ

knitr::opts_chunk$set(prompt = TRUE, comment = NA)

hmkzLAPV

R balík pro výuku předmětu Hydrologické dopady klimatické změny (FŽP, ČZU) v roce 2016.

Instalace

library(devtools)
install_github("hanel/hmkzLAPV")

nebo

library(devtools)
install_git("https://github.com/hanel/hmkzLAPV.git")

Úvod

Informace o jednotlivých LAPV jsou dostupné v tabulce lapv_tab:

library(hmkzLAPV)
data(lapv_tab)
head(lapv_tab)

Data pro jednotlivá LAPV načtete pomocí funkce get_lapv_data, jejímž argumentem je ID, tj. identifikátor povodí:

dta = get_lapv_data("AMERIKA")
head(dta)

Data obsahují:

tab = data.table(names(dta), c('datum', 'ID povodí', 'simulované srážky', 'simulovaná teplota', 'globální klimatický model', 'ID běhu modelu', 'ID simulace', 'pozorované srážky', 'pozorovaná teplota', 'odtok pro pozorované podmínky (simulace modelu Bilan) '))
knitr::kable(tab, col.names = c('název', 'význam'), format = 'markdown')

Zadání:

1. CVIČENÍ

• V rychlosti prostudujte článek Hurst, H. E. (1956). The Problem of long-term storage in reservoirs. International Association of Scientific Hydrology. Bulletin, 1:3, 13-27. DOI: 10.1080/02626665609493644

(A).

• Zagregujte data do ročního časového kroku (využijte funkci mon2yr)
• Pro různá N z intervalu (10, 110) vypočtěte K1 = R / (sigma N^0.5) a K2 = log (R / sigma) / log (N/2), kde R je rozpětí hodnot a sigma je směrodatná odchylka
• Pro různá N z intervalu (10, 110) vykreslete log(R/sigma) proti log(N/2)
• Vytvořte lineární model lm(log(R/sigma) ~ log(N/2)), zjistěte hodnotu regresního koeficientu (h) a dokreslete přímku do grafu
• opakujte pro bílý šum, AR(1) [případně AR(2), AR(3), atd.] proces a FGN proces - pro generování veličin použijte funkci toyGen
• výsledky porovnejte
• vytvořte funkci, která spočítá pro zadanou veličinu koeficient h
• !! vše uložte !!

(B)onus

devtools::install_github("jbkunst/highcharter")
ydta = mon2yr(dta)
library(highcharter)
library(magrittr)
highchart(type = "stock") %>% hc_tooltip(valueDecimals = 2) %>% 
  hc_add_series_times_values(ydta$DTM, ydta$obs_P, name = "Pozorovaná srážka", showInLegend = TRUE) %>% 
  hc_add_series_times_values(ydta$DTM, ydta$sim_P, name = "Simulovaná srážka") %>%
  hc_add_theme(hc_theme_darkunica())

viz http://jkunst.com/highcharter

(*)

• S využitím funkce agg2 postupně agregujte data po dvojicích
• Vykreslete graf log(k) proti log(sd(x) ^ 0.5), kde k je agregace a sd(x) je veličina agregovaná na k-té úrovni

---

2. CVIČENÍ

(A).

• funkce hurst je nyní součástí balíku, v čem spočívají zásadní rozdíly mezi touto funkcí a (chybnou) funkcí, kterou jsme používali minule? Tj.

hurst = function(x, N){
  res = data.frame(N = 10:100  )
  for (i in 1:nrow(res) ) {
    sx = sample(x, res[i, "N"])
    R = diff(range(sx))
    s = sd(sx)
    N = res[i, "N"]
    res[i, "R"] = R
    res[i, "sigma"] = s
    res[i, "k1"] = R / (s * N^.5)
    res[i, "k2"] = log(R/s) / log(N/2)
  }
  return(res)
}

• na základě vybrané veličiny x (např. pozorovaný odtok) simulujte pomocí funkce toyGen veličiny s využitím bílého šumu (WN), autoregresního procesu prvního řádu (AR1), a fractional gaussian noise (FGN). Vykreslete tyto řady a řady kumulativních odchylek od průměru (použijte funkci plotGen)
• jaké jsou typické hodnoty Hurstova koeficientu pro WN, AR1 a FGN? (lze zjistit pomocí attr(hurst(x), 'h'))?

• pomocí opakovaného samplování určete 90% interval spolehlivosti pro odhady Hurstova koeficientu pro jednotlivé procesy

• vykreslete časovou řadu (ročních) pozorovaných a simulovaných srážek nebo teploty a jejich kumulativní odchylky od průměru

• jak se liší odhady Hurstova koeficientu? Jsou rozdíly statisticky významné?

(B)onus

Někdy by bylo výhodné, aby grafy v RStudiu umožňovaly alespoň omezenou interaktivitu. Jednou z možností je využítí balíku manipulate, který umožňuje velmi jednoduchou implementaci ovládacích prvků (slider, picker, checkbox a button)

Příklady s využitím dat z balíku hmkzLAPV:

př1:

library(manipulate)

dta = get_lapv_data('AMERIKA')
mdta = mon2yr(dta)

manipulate(
  plot(mdta$DTM, y, type = type),
  y = picker("sim_P" = mdta$sim_P, "obs_P" = mdta$obs_P),
  type = picker("Points" = "p", "Line" = "l", "Step" = "s")
)

př2:

manipulate(
  {plot(mdta$DTM, y1, type = 'l', col = col1, ylim = range(y1, y2))
   lines(mdta$DTM, y2, col = col2)  
    },
    y1 = picker("sim_P" = mdta$sim_P, "obs_P" = mdta$obs_P),
    y2 = picker("sim_P" = mdta$sim_P, "obs_P" = mdta$obs_P),
    col1 = slider(1, 650, initial = 300),
    col2 = slider(1, 650, initial = 200)
)

př3:

manipulate(
  plot(mdta$DTM, y, type = type, xlim = c(mdta$DTM[from], mdta$DTM[to])),
  y = picker("sim_P" = mdta$sim_P, "obs_P" = mdta$obs_P),
  type = picker("Points" = "p", "Line" = "l", "Step" = "s"), 
  from = slider(1, nrow(mdta)),
  to = slider(1, nrow(mdta), initial = nrow(mdta))
)

(*)

• S využitím funkce agg2 postupně agregujte data po dvojicích
• Vykreslete graf log(k) proti log(sd(x) ^ 0.5), kde k je agregace a sd(x) je veličina agregovaná na k-té úrovni

---

3. CVIČENÍ - KOREKCE SYSTEMATICKÝCH CHYB

(A).

1. Zjistěte, jak se liší průměrné srážky a teplota v simulaci od pozorování
2. Vykreslete empirické distribuční funkce pozorovaných a simulovaných srážek/teploty pro všechny data a jednotlivé měsíce
3. Zjistěte, jaká je chyba v distribuční funkci simulovaných srážek a teploty pro měsíční (všechny vs jednotlivé měsíce) a roční data
4. Vyhlaďte rozdíly pomocí filtru loess
5. Opravte simulované veličiny
6. Vytvořte funkci umožňující korekci distribuční funkce pro jednotlivé měsíce

---

4. CVIČENÍ - HYDROLOGICKÝ MODEL A SYNTÉZA

Nové funkce v balíku:

• correct - provede korekci systematických chyb - viz ?correct
• bil.lapv - vytvoří instanci modelu Bilan s parametry pro zadanou LAPV - viz ?bil.lapv

Stručný úvod do modelu Bilan

• Bilan je model hydrologické bilance vyvíjený ve Výzkumném ústavu vodohospodářském T. G. Masaryka, v.v.i.
• 2 verze (stejné jádro) - uživatelské rozhraní x balík pro R
• nainstalujte pomocí zipového souboru bilan_2015-06-18.zip ve složce materialy, verze pro linux a OS X viz bilan_2015-06-18.tar.gz
• základní postup práce relevantní pro naše účely je následující:

library(bilan)

# načti data
dta = get_lapv_data("AMERIKA")

# vytvoř model 
b = bil.lapv("AMERIKA")

# model existuje, má příslušné parametry:
bil.get.params(b)

# ale žádná data
bil.get.data(b)

# Proto nahrajeme do modelu příslušná data z data.tablu dta
bil.set.values(b, dta[, .(DTM, P = obs_P, T = obs_T)])

# poté stačí spočítat potenciální evapotranspiraci 
bil.pet(b)

# a model spustit
res = bil.run(b)

res

Zadání:

• nahrajte data pro zvolenou LAPV
• pomocí funkce correct zkorigujte simulované srážky a teploty
• nainstalujte balík bilan
• nahrajte model pro zvolenou LAPV
• vytvořte data.frame, který bude obsahovat veličiny
  • DTM - datum
  • obs_RM - odtok modelovaný na základě pozorovaných srážek a teploty
  • sim_RM - odtok modelovaný na základě simulovaných srážek a teploty z klimatického modelu
  • cor_RM - odtok modelovaný na základě korigovaných srážek a teploty
• vykreslete odchylky sim_RM/obs_RM a cor_RM/obs_RM pro jednotlivé kvantily distribuční funkce a měsíce
• spočítejte hurstův koeficient pro tyto řady a vykreslete řady kumulativních odchylek od průměru

PROTOKOL

• použijte Rmarkdown
• struktura viz protokol.Rmd (zkopírujte a doplňte)
• odevzdání - mailem pdf nebo odkaz na html uložené na rpubs.com

hanel/hmkzLAPV documentation built on May 17, 2019, 2:28 p.m.