From cc2c21b1a6aaa32511ee29cdbfaddd7ea95355e6 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 18 Oct 2023 13:23:26 +0200 Subject: [PATCH 01/24] added shiny lakes and tb indepth, these will be combined in a future commit --- vignettes/tblakes_shiny.Rmd | 47 +++++++++++++++++++++++++++++++++++++ vignettes/tibetan-lakes.Rmd | 0 2 files changed, 47 insertions(+) create mode 100644 vignettes/tblakes_shiny.Rmd create mode 100644 vignettes/tibetan-lakes.Rmd diff --git a/vignettes/tblakes_shiny.Rmd b/vignettes/tblakes_shiny.Rmd new file mode 100644 index 00000000..bddf9da6 --- /dev/null +++ b/vignettes/tblakes_shiny.Rmd @@ -0,0 +1,47 @@ +--- +title: "Tibetan Lakes Shiny Example" +output: + bookdown::html_document2: + base_format: rmarkdown::html_vignette + fig_caption: yes +bibliography: "covid.bib" +link-citations: yes +vignette: > + %\VignetteIndexEntry{Tibetan Lakes Shiny Example} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +Here we will show how to use the shiny apps with the tibetan lakes dataset. + +```{r error=TRUE} +devtools::load_all() +library(INLA) +library(inlabru) +data <- read.csv('~/repos/tblakes/data_for_modeling.csv') +summary(data) +``` + + +now we will use the mesh builder app: + +```{r error=TRUE} +fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") +``` + +using the meshbuilder we have generated the following code: + +```{r error=TRUE} +mesh <- inla.mesh.2d(loc = data[c('LONG', 'LAT')], + max.edge = c(1.1, 8.1), + cutoff = 0.81, + offset=c(1.2, 4.6)) +plot(mesh) +``` + +next we will use the model tester app: + +```{r error=TRUE} +sp::coordinates(data) <- c('LONG', 'LAT') +fdmr::interactive_priors(data, measurement_data = data, time_variable = data$year) +``` \ No newline at end of file diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd new file mode 100644 index 00000000..e69de29b From 3dc3c004d548516a5a3b720e5a38bf75a843007b Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 18 Oct 2023 13:33:23 +0200 Subject: [PATCH 02/24] added some discussion --- vignettes/covid.bib | 52 +++++++++++++++++++++++++++++++++++++ vignettes/tibetan-lakes.Rmd | 15 +++++++++++ 2 files changed, 67 insertions(+) diff --git a/vignettes/covid.bib b/vignettes/covid.bib index 7f2bb2ee..2e47136f 100755 --- a/vignettes/covid.bib +++ b/vignettes/covid.bib @@ -1,3 +1,55 @@ +@article{brun2020limited, + title={Limited contribution of glacier mass loss to the recent increase in Tibetan Plateau lake volume}, + author={Brun, Fanny and Treichler, D{\'e}sir{\'e}e and Shean, David and Immerzeel, Walter W}, + journal={Frontiers in Earth Science}, + volume={8}, + pages={582060}, + year={2020}, + publisher={Frontiers Media SA}, + doi={10.3389/feart.2020.582060} +} + +@article{twentythreeproblems, +author = {Günter Blöschl, Marc F.P. Bierkens, Antonio Chambel, Christophe Cudennec, Georgia Destouni, Aldo Fiori, James W. Kirchner, Jeffrey J. McDonnell, Hubert H.G. Savenije, Murugesu Sivapalan, Christine Stumpp, Elena Toth, Elena Volpi, Gemma Carr, Claire Lupton, Josè Salinas, Borbála Széles, Alberto Viglione, Hafzullah Aksoy, Scott T. Allen, Anam Amin, Vazken Andréassian, Berit Arheimer, Santosh K. Aryal, Victor Baker, Earl Bardsley, Marlies H. Barendrecht, Alena Bartosova, Okke Batelaan, Wouter R. Berghuijs, Keith Beven, Theresa Blume, Thom Bogaard, Pablo Borges de Amorim, Michael E. Böttcher, Gilles Boulet, Korbinian Breinl, Mitja Brilly, Luca Brocca, Wouter Buytaert, Attilio Castellarin, Andrea Castelletti, Xiaohong Chen, Yangbo Chen, Yuanfang Chen, Peter Chifflard, Pierluigi Claps, Martyn P. Clark, Adrian L. Collins, Barry Croke, Annette Dathe, Paula C. David, Felipe P. J. de Barros, Gerrit de Rooij, Giuliano Di Baldassarre, Jessica M. Driscoll, Doris Duethmann, Ravindra Dwivedi, Ebru Eris, William H. Farmer, James Feiccabrino, Grant Ferguson, Ennio Ferrari, Stefano Ferraris, Benjamin Fersch, David Finger, Laura Foglia, Keirnan Fowler, Boris Gartsman, Simon Gascoin, Eric Gaume, Alexander Gelfan, Josie Geris, Shervan Gharari, Tom Gleeson, Miriam Glendell, Alena Gonzalez Bevacqua, María P. González-Dugo, Salvatore Grimaldi, A. B. Gupta, Björn Guse, Dawei Han, David Hannah, Adrian Harpold, Stefan Haun, Kate Heal, Kay Helfricht, Mathew Herrnegger, Matthew Hipsey, Hana Hlaváčiková, Clara Hohmann, Ladislav Holko, Christopher Hopkinson, Markus Hrachowitz, Tissa H. Illangasekare, Azhar Inam, Camyla Innocente, Erkan Istanbulluoglu, Ben Jarihani, Zahra Kalantari, Andis Kalvans, Sonu Khanal, Sina Khatami, Jens Kiesel, Mike Kirkby, Wouter Knoben, Krzysztof Kochanek, Silvia Kohnová, Alla Kolechkina, Stefan Krause, David Kreamer, Heidi Kreibich, Harald Kunstmann, Holger Lange, Margarida L. R. Liberato, Eric Lindquist, Timothy Link, Junguo Liu, Daniel Peter Loucks, Charles Luce, Gil Mahé, Olga Makarieva, Julien Malard, Shamshagul Mashtayeva, Shreedhar Maskey, Josep Mas-Pla, Maria Mavrova-Guirguinova, Maurizio Mazzoleni, Sebastian Mernild, Bruce Dudley Misstear, Alberto Montanari, Hannes Müller-Thomy, Alireza Nabizadeh, Fernando Nardi, Christopher Neale, Nataliia Nesterova, Bakhram Nurtaev, Vincent O. Odongo, Subhabrata Panda, Saket Pande, Zhonghe Pang, Georgia Papacharalampous, Charles Perrin, Laurent Pfister, Rafael Pimentel, María J. Polo, David Post, Cristina Prieto Sierra, Maria-Helena Ramos, Maik Renner, José Eduardo Reynolds, Elena Ridolfi, Riccardo Rigon, Monica Riva, David E. Robertson, Renzo Rosso, Tirthankar Roy, João H.M. Sá, Gianfausto Salvadori, Mel Sandells, Bettina Schaefli, Andreas Schumann, Anna Scolobig, Jan Seibert, Eric Servat, Mojtaba Shafiei, Ashish Sharma, Moussa Sidibe, Roy C. Sidle, Thomas Skaugen, Hugh Smith, Sabine M. Spiessl, Lina Stein, Ingelin Steinsland, Ulrich Strasser, Bob Su, Jan Szolgay, David Tarboton, Flavia Tauro, Guillaume Thirel, Fuqiang Tian, Rui Tong, Kamshat Tussupova, Hristos Tyralis, Remko Uijlenhoet, Rens van Beek, Ruud J. van der Ent, Martine van der Ploeg, Anne F. Van Loon, Ilja van Meerveld, Ronald van Nooijen, Pieter R. van Oel, Jean-Philippe Vidal, Jana von Freyberg, Sergiy Vorogushyn, Przemyslaw Wachniew, Andrew J. Wade, Philip Ward, Ida K. Westerberg, Christopher White, Eric F. Wood, Ross Woods, Zongxue Xu, Koray K. Yilmaz and Yongqiang Zhang}, +title = {Twenty-three unsolved problems in hydrology (UPH) – a community perspective}, +journal = {Hydrological Sciences Journal}, +volume = {64}, +number = {10}, +pages = {1141-1158}, +year = {2019}, +publisher = {Taylor & Francis}, +doi = {10.1080/02626667.2019.1620507}, +} + + +@article{immerzeel2020importance, + title={Importance and vulnerability of the world’s water towers}, + author={Immerzeel, Walter W and Lutz, Arthur F and Andrade, Marcos and Bahl, A and Biemans, Hester and Bolch, Tobias and Hyde, Sam and Brumby, S and Davies, BJ and Elmore, AC and others}, + journal={Nature}, + volume={577}, + number={7790}, + pages={364--369}, + year={2020}, + publisher={Nature Publishing Group UK London}, + doi={10.1038/s41586-019-1822-y} + } + +@article{zhang2017grl, +author = {Zhang, Guoqing and Yao, Tandong and Shum, C. K. and Yi, Shuang and Yang, Kun and Xie, Hongjie and Feng, Wei and Bolch, Tobias and Wang, Lei and Behrangi, Ali and Zhang, Hongbo and Wang, Weicai and Xiang, Yang and Yu, Jinyuan}, +title = {Lake volume and groundwater storage variations in Tibetan Plateau's endorheic basin}, +journal = {Geophysical Research Letters}, +volume = {44}, +number = {11}, +pages = {5550-5560}, +keywords = {lake volume, mass balance, groundwater storage, Tibetan Plateau}, +doi = {https://doi.org/10.1002/2017GL073773}, +url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL073773}, +eprint = {https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017GL073773}, +abstract = {Abstract The Tibetan Plateau (TP), the highest and largest plateau in the world, with complex and competing cryospheric-hydrologic-geodynamic processes, is particularly sensitive to anthropogenic warming. The quantitative water mass budget in the TP is poorly known. Here we examine annual changes in lake area, level, and volume during 1970s–2015. We find that a complex pattern of lake volume changes during 1970s–2015: a slight decrease of −2.78 Gt yr−1 during 1970s–1995, followed by a rapid increase of 12.53 Gt yr−1 during 1996–2010, and then a recent deceleration (1.46 Gt yr−1) during 2011–2015. We then estimated the recent water mass budget for the Inner TP, 2003–2009, including changes in terrestrial water storage, lake volume, glacier mass, snow water equivalent (SWE), soil moisture, and permafrost. The dominant components of water mass budget, namely, changes in lake volume (7.72 ± 0.63 Gt yr−1) and groundwater storage (5.01 ± 1.59 Gt yr−1), increased at similar rates. We find that increased net precipitation contributes the majority of water supply (74\%) for the lake volume increase, followed by glacier mass loss (13\%), and ground ice melt due to permafrost degradation (12\%). Other term such as SWE (1\%) makes a relatively small contribution. These results suggest that the hydrologic cycle in the TP has intensified remarkably during recent decades.}, +year = {2017} +} + + @misc{who2020world, title={World {H}ealth {O}rganization: {C}oronavirus {D}isease 2019 ({COVID}-19) {S}ituation {R}eport}, diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index e69de29b..d22136b7 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -0,0 +1,15 @@ +--- +title: "Tibetan Lakes Example" +output: + bookdown::html_document2: + base_format: rmarkdown::html_vignette + fig_caption: yes +bibliography: "covid.bib" +link-citations: yes +vignette: > + %\VignetteIndexEntry{Tibetan Lakes Example} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +The Tibetan Plateau and the High Mountain Asia region provide water resources to approximately 1 billion people and provide for their health, economies, and agriculture. These natural “water towers” are highly sensitive to climate change and therefore changes in regional water resources must be understood (@immerzeel2020importance). Indeed, understanding how climate change will alter cold region water resources has been described as a “grand challenge” of hydrology (@twentythreeproblems). Recent changes in lakes in the Tibetan Plateau endorheic basin have led to increases in lake areas. Conventionally attributed to increases in glacier melt due to climate change, a recent investigation has demonstrated that glacier mass loss contributes only a small fraction (19%) to the lake volume increase (@brun2020limited). Thus, an unresolved question in the Tibetan Plateau endorheic basin lakes is, what physical processes drive lake growth in the region and at what fraction? From ef0371e6d45bf3c489cf3314a92fc6cd524e7ce0 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 18 Oct 2023 14:16:08 +0200 Subject: [PATCH 03/24] added refs for datasets --- vignettes/covid.bib | 36 ++++++++++++++++++++++++++++++++++++ vignettes/tibetan-lakes.Rmd | 12 ++++++++++++ 2 files changed, 48 insertions(+) diff --git a/vignettes/covid.bib b/vignettes/covid.bib index 2e47136f..28a9464e 100755 --- a/vignettes/covid.bib +++ b/vignettes/covid.bib @@ -1,3 +1,39 @@ +@article{hersbach2020era5, + title={The ERA5 global reanalysis}, + author={Hersbach, Hans and Bell, Bill and Berrisford, Paul and Hirahara, Shoji and Hor{\'a}nyi, Andr{\'a}s and Mu{\~n}oz-Sabater, Joaqu{\'\i}n and Nicolas, Julien and Peubey, Carole and Radu, Raluca and Schepers, Dinand and others}, + journal={Quarterly Journal of the Royal Meteorological Society}, + volume={146}, + number={730}, + pages={1999--2049}, + year={2020}, + publisher={Wiley Online Library}, + doi={10.1002/qj.3803} +} + +@article{lehner2013global, + title={Global river hydrography and network routing: baseline data and new approaches to study the world's large river systems}, + author={Lehner, Bernhard and Grill, G{\"u}nther}, + journal={Hydrological Processes}, + volume={27}, + number={15}, + pages={2171--2186}, + year={2013}, + publisher={Wiley Online Library}, + doi={doi.org/10.1002/hyp.9740} +} + +@article{pekel2016high, + title={High-resolution mapping of global surface water and its long-term changes}, + author={Pekel, Jean-Fran{\c{c}}ois and Cottam, Andrew and Gorelick, Noel and Belward, Alan S}, + journal={Nature}, + volume={540}, + number={7633}, + pages={418--422}, + year={2016}, + publisher={Nature Publishing Group UK London}, + doi={10.1038/nature20584} +} + @article{brun2020limited, title={Limited contribution of glacier mass loss to the recent increase in Tibetan Plateau lake volume}, author={Brun, Fanny and Treichler, D{\'e}sir{\'e}e and Shean, David and Immerzeel, Walter W}, diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index d22136b7..8ada5830 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -13,3 +13,15 @@ vignette: > --- The Tibetan Plateau and the High Mountain Asia region provide water resources to approximately 1 billion people and provide for their health, economies, and agriculture. These natural “water towers” are highly sensitive to climate change and therefore changes in regional water resources must be understood (@immerzeel2020importance). Indeed, understanding how climate change will alter cold region water resources has been described as a “grand challenge” of hydrology (@twentythreeproblems). Recent changes in lakes in the Tibetan Plateau endorheic basin have led to increases in lake areas. Conventionally attributed to increases in glacier melt due to climate change, a recent investigation has demonstrated that glacier mass loss contributes only a small fraction (19%) to the lake volume increase (@brun2020limited). Thus, an unresolved question in the Tibetan Plateau endorheic basin lakes is, what physical processes drive lake growth in the region and at what fraction? + +```{r error=TRUE} +library(devtools) +devtools::load_all() +# fdmr::retrieve_tutorial_data(dataset = "tibetan lakes") + +# TODO : move this data into the Rdataset repo +data <- read.csv('~/repos/tblakes/data_for_modeling.csv') +summary(data) +``` + +We have compiled a dataset for 20 years of lake growth using the global surface water survey (@pekel2016high), the hydrobasins catchement dataset (@lehner2013global), and ERA5-land precipitation data (@hersbach2020era5). From 365766823617c3180491555281a2d71b60bddd29 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 18 Oct 2023 16:26:17 +0200 Subject: [PATCH 04/24] added plot for precip vs lake growth --- vignettes/tibetan-lakes.Rmd | 45 +++++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 8ada5830..4bf9b28a 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -20,8 +20,49 @@ devtools::load_all() # fdmr::retrieve_tutorial_data(dataset = "tibetan lakes") # TODO : move this data into the Rdataset repo -data <- read.csv('~/repos/tblakes/data_for_modeling.csv') +data <- read.csv('~/repos/tblakes/data_for_4dm.csv') summary(data) ``` -We have compiled a dataset for 20 years of lake growth using the global surface water survey (@pekel2016high), the hydrobasins catchement dataset (@lehner2013global), and ERA5-land precipitation data (@hersbach2020era5). +We have compiled a dataset for 20 years of lake growth using the global surface water survey (@pekel2016high), the hydrobasins catchment dataset (@lehner2013global), and ERA5-land precipitation data (@hersbach2020era5). The global surface water survey is a data set using LANDSAT data to classify pixels where lakes are annually. We have converted this to a time series by looking at how lakes have changed between years. This is done in Google Earth Engine [using the following codes](https://code.earthengine.google.com/3304560943f796c82f935403dcc54909). We then aggregate data to changes that have occured within a catchment using the hydrobasins catchment definitions (level 6). We do the same for the ERA5-land precipitation data. This is compiled into a CSV that is presented here. + +Let's look at this data and see what is going on. First lets look at the overall changes in lake areas. + + +```{r error=TRUE, fig.width=8, fig.height=4, fig.align = "center"} +library(tidyverse) + +# Split data into positive and negative, then summarize +summary_data <- data %>% + mutate(sign = ifelse(water_balance_m3 > 0, "Lake Growth", "Lake Decline")) %>% + group_by(year, sign) %>% + summarise(total = sum(water_balance_m3, na.rm = TRUE)) %>% + ungroup() + +# Summarize the data by year +sums <- data %>% + group_by(year) %>% + summarise( + precip_sum = sum(precip, na.rm = TRUE), + growth_ratio_sum = sum(growth_ratio, na.rm = TRUE), + decline_ratio_sum = sum(decline_ratio, na.rm = TRUE) + ) + +# Create the combined plot +ggplot() + + geom_col(data = summary_data, aes(x = year, y = total, fill = sign), position = "stack") + + geom_line(data = sums, aes(x = year, y = precip_sum * 100), color = 'cyan') + # Scale to match bar plot magnitude + scale_y_continuous( + name = "Water Balance (m3) and Precipitation (x100)", + sec.axis = sec_axis(~./10, name = "Precipitation") + ) + + theme_minimal(base_size = 15) + + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + + labs(fill = "Sign") +``` +We can see here that while lakes both grow and decline, the overall trend is that lakes are growing more than they are shrinking in the Tibetan plateau. Why is this? perhaps it is because of the rain. Our ERA5-land data set gives us precipitation per year (after aggregation). + + + + + From c1f70d250d7337c750c1d325f7ef959c8a5f7d21 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Thu, 19 Oct 2023 11:12:30 +0200 Subject: [PATCH 05/24] added a map --- vignettes/tibetan-lakes.Rmd | 57 ++++++++++++++++++++++++++++--------- 1 file changed, 43 insertions(+), 14 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 4bf9b28a..22ff88e3 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -12,7 +12,7 @@ vignette: > %\VignetteEncoding{UTF-8} --- -The Tibetan Plateau and the High Mountain Asia region provide water resources to approximately 1 billion people and provide for their health, economies, and agriculture. These natural “water towers” are highly sensitive to climate change and therefore changes in regional water resources must be understood (@immerzeel2020importance). Indeed, understanding how climate change will alter cold region water resources has been described as a “grand challenge” of hydrology (@twentythreeproblems). Recent changes in lakes in the Tibetan Plateau endorheic basin have led to increases in lake areas. Conventionally attributed to increases in glacier melt due to climate change, a recent investigation has demonstrated that glacier mass loss contributes only a small fraction (19%) to the lake volume increase (@brun2020limited). Thus, an unresolved question in the Tibetan Plateau endorheic basin lakes is, what physical processes drive lake growth in the region and at what fraction? +The Tibetan Plateau and the High Mountain Asia region provide water resources to approximately 1 billion people and provide for their health, economies, and agriculture. These natural "water towers" are highly sensitive to climate change and therefore changes in regional water resources must be understood (@immerzeel2020importance). Indeed, understanding how climate change will alter cold region water resources has been described as a "grand challenge" of hydrology (@twentythreeproblems). Recent changes in lakes in the Tibetan Plateau endorheic basin have led to increases in lake areas. Conventionally attributed to increases in glacier melt due to climate change, a recent investigation has demonstrated that glacier mass loss contributes only a small fraction (19%) to the lake volume increase (@brun2020limited). Thus, an unresolved question in the Tibetan Plateau endorheic basin lakes is, what physical processes drive lake growth in the region and at what fraction? ```{r error=TRUE} library(devtools) @@ -28,7 +28,6 @@ We have compiled a dataset for 20 years of lake growth using the global surface Let's look at this data and see what is going on. First lets look at the overall changes in lake areas. - ```{r error=TRUE, fig.width=8, fig.height=4, fig.align = "center"} library(tidyverse) @@ -48,21 +47,51 @@ sums <- data %>% decline_ratio_sum = sum(decline_ratio, na.rm = TRUE) ) -# Create the combined plot +# TODO : combine these plots +ggplot() + + geom_col(data = summary_data, aes(x = year, y = total, fill = sign), position = "stack") + ggplot() + - geom_col(data = summary_data, aes(x = year, y = total, fill = sign), position = "stack") + - geom_line(data = sums, aes(x = year, y = precip_sum * 100), color = 'cyan') + # Scale to match bar plot magnitude - scale_y_continuous( - name = "Water Balance (m3) and Precipitation (x100)", - sec.axis = sec_axis(~./10, name = "Precipitation") - ) + - theme_minimal(base_size = 15) + - theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + - labs(fill = "Sign") + geom_line(data = sums, aes(x = year, y = precip_sum * 100), color = 'cyan') # Scale +``` + +We can see here that while lakes both grow and decline, the overall trend is that lakes are growing more than they are shrinking in the Tibetan plateau. Why is this? perhaps it is because of the rain. Our ERA5-land data set gives us precipitation per year (after aggregation). But we see that the largest year that lakes grow, is not the largest amount of rain. What is going on here? + +```{r error=TRUE, fig.width=8, fig.height=8, fig.align = "center"} +ggplot(data, aes(x = precip, y = water_balance_m3, color = as.factor(year))) + + geom_point(alpha = 0.5) + + scale_color_viridis_d() + + theme_minimal() +``` + +what if we make a map? + +```{r error=TRUE, fig.width=8, fig.height=8, fig.align='center'} +library(tidyverse) +library(leaflet) + +# Summarize the data +summary_data <- data %>% + group_by(centroid_lon, centroid_lat) %>% + summarise( + total_balance = sum(abs(water_balance_m3), na.rm = TRUE) + ) + +# Create a map +leaflet(summary_data) %>% + addTiles() %>% + addCircles( + lng = ~centroid_lon, + lat = ~centroid_lat, + radius = ~50*total_balance, + color = "blue", + stroke = FALSE, + fillOpacity = 0.6 + ) ``` -We can see here that while lakes both grow and decline, the overall trend is that lakes are growing more than they are shrinking in the Tibetan plateau. Why is this? perhaps it is because of the rain. Our ERA5-land data set gives us precipitation per year (after aggregation). - + + From 41ed543033b6a92f88d2a10bb0a2cd8434f34559 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Thu, 19 Oct 2023 11:13:59 +0200 Subject: [PATCH 06/24] fixed bug in calculation for map --- vignettes/tibetan-lakes.Rmd | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 22ff88e3..c374d215 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -64,8 +64,7 @@ ggplot(data, aes(x = precip, y = water_balance_m3, color = as.factor(year))) + theme_minimal() ``` -what if we make a map? - +Here we can see a clear spatial distinction for regions that have had large lake growth compared to small lake growth. ```{r error=TRUE, fig.width=8, fig.height=8, fig.align='center'} library(tidyverse) library(leaflet) @@ -74,7 +73,7 @@ library(leaflet) summary_data <- data %>% group_by(centroid_lon, centroid_lat) %>% summarise( - total_balance = sum(abs(water_balance_m3), na.rm = TRUE) + total_balance = sum(water_balance_m3, na.rm = TRUE) ) # Create a map From 880b44b61453f1fd3290ab0090a36c46a5266b13 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Thu, 19 Oct 2023 11:19:01 +0200 Subject: [PATCH 07/24] added TODOs --- vignettes/tibetan-lakes.Rmd | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index c374d215..1538e073 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -66,6 +66,8 @@ ggplot(data, aes(x = precip, y = water_balance_m3, color = as.factor(year))) + Here we can see a clear spatial distinction for regions that have had large lake growth compared to small lake growth. ```{r error=TRUE, fig.width=8, fig.height=8, fig.align='center'} + +# TODO : this should include plotting the polygons of the catchments library(tidyverse) library(leaflet) @@ -89,6 +91,8 @@ leaflet(summary_data) %>% ) ``` + + From 2c2e2055b31aa5098a1cb215e2200c5c0b2d8da1 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Thu, 19 Oct 2023 11:29:41 +0200 Subject: [PATCH 08/24] added mesh --- vignettes/tibetan-lakes.Rmd | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 1538e073..e1d04354 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -64,7 +64,7 @@ ggplot(data, aes(x = precip, y = water_balance_m3, color = as.factor(year))) + theme_minimal() ``` -Here we can see a clear spatial distinction for regions that have had large lake growth compared to small lake growth. +Here we can see across the 20 year observation period a clear spatial distinction in lakegrowth however it is not as clear cut for precipitation. ```{r error=TRUE, fig.width=8, fig.height=8, fig.align='center'} # TODO : this should include plotting the polygons of the catchments @@ -75,7 +75,8 @@ library(leaflet) summary_data <- data %>% group_by(centroid_lon, centroid_lat) %>% summarise( - total_balance = sum(water_balance_m3, na.rm = TRUE) + total_balance = sum(water_balance_m3, na.rm = TRUE), + total_precip = sum(precip, na.rm = TRUE) ) # Create a map @@ -88,13 +89,32 @@ leaflet(summary_data) %>% color = "blue", stroke = FALSE, fillOpacity = 0.6 + ) %>% + addCircles( + lng = ~centroid_lon, + lat = ~centroid_lat, + radius = ~500*total_precip, + color = "red", + stroke = FALSE, + fillOpacity = 0.6 ) ``` +Let's try and model this. First we will use the meshbuilder to build a reasonable mesh. - - +```{r error=TRUE} - +fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", longitude_column = "centroid_lon", latitude_column = "centroid_lat") + +``` + +Using the mesh_builder we have the code to build the mesh for the inlabru model: + +```{r error=TRUE} +mesh <- fmesher::fm_mesh_2d_inla(loc = location_data, + max.edge = c(3, 7.9), + cutoff = 0.9, + offset=c(0.2, 4)) +``` \ No newline at end of file From a8460a7d22e56c3cf7bb0730f51c48d639c52fea Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Thu, 19 Oct 2023 12:02:00 +0200 Subject: [PATCH 09/24] added model eval, not working --- vignettes/tibetan-lakes.Rmd | 1 + 1 file changed, 1 insertion(+) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index e1d04354..b8be22a8 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -117,4 +117,5 @@ mesh <- fmesher::fm_mesh_2d_inla(loc = location_data, max.edge = c(3, 7.9), cutoff = 0.9, offset=c(0.2, 4)) +plot(mesh) ``` \ No newline at end of file From 6f98f5b5eb0c5c8330f575299e336858c693e1cb Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Sun, 22 Oct 2023 09:59:40 +0200 Subject: [PATCH 10/24] added text about 4d-modeller and inlabru --- vignettes/covid.bib | 3 ++- vignettes/tibetan-lakes.Rmd | 15 +++++++++++---- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/vignettes/covid.bib b/vignettes/covid.bib index 28a9464e..0088a870 100755 --- a/vignettes/covid.bib +++ b/vignettes/covid.bib @@ -1,4 +1,5 @@ -@article{hersbach2020era5, +@article{ +, title={The ERA5 global reanalysis}, author={Hersbach, Hans and Bell, Bill and Berrisford, Paul and Hirahara, Shoji and Hor{\'a}nyi, Andr{\'a}s and Mu{\~n}oz-Sabater, Joaqu{\'\i}n and Nicolas, Julien and Peubey, Carole and Radu, Raluca and Schepers, Dinand and others}, journal={Quarterly Journal of the Royal Meteorological Society}, diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index b8be22a8..aebf4ff0 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -100,9 +100,13 @@ leaflet(summary_data) %>% ) ``` -Let's try and model this. First we will use the meshbuilder to build a reasonable mesh. +Let's try and model this. We will use 4D-Modeller and inlabru. +[inlabru]() is designed to calculate fixed and random effects as well as continuous spatially and temporally distributed processes using SPDEs (). In order to make these computationally solveable, the continous processes are discretized on a finite element mesh. This mesh represents the spatial distribution of the process under study. At each node of the mesh, the model is used to calculate the outcome variable (in this case lake area changes). The mesh provides the spatial awareness to the model. Different processes can have different length scales, i.e., how far away does a process need to occur before it has no effect on the occurence of the process at my current location. These are controlled by priors (see below). +4D-Modeller is designed to make these model design decisions visual, interpretable, and accessible to users who may not have a strong background in R programming or bayesian spatiotemporal modeling. 4D-modeller has a set of shiny apps that help the user write the codes necessary to implement the full model and evaluate it's performance. + +First we will use the mesh_builder to build a reasonable mesh. ```{r error=TRUE} @@ -113,9 +117,12 @@ fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ell Using the mesh_builder we have the code to build the mesh for the inlabru model: ```{r error=TRUE} -mesh <- fmesher::fm_mesh_2d_inla(loc = location_data, +mesh <- fmesher::fm_mesh_2d_inla(loc = data['centroid_lat', 'centroid_lon'], max.edge = c(3, 7.9), cutoff = 0.9, offset=c(0.2, 4)) -plot(mesh) -``` \ No newline at end of file +# plot(mesh) +mesh +``` + + From 6e2814aa1d3fe559a6b43c0b17e938461807bf71 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Sun, 22 Oct 2023 10:44:56 +0200 Subject: [PATCH 11/24] added meshbuilder stuff --- vignettes/tibetan-lakes.Rmd | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index aebf4ff0..9d56e7ae 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -114,15 +114,20 @@ fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ell ``` -Using the mesh_builder we have the code to build the mesh for the inlabru model: +Using the mesh_builder we have the code to build the mesh for the inlabru model. Note, we needed to replace `location data` with `data[, c('centroid_lon', 'centroid_lat')]`. ```{r error=TRUE} -mesh <- fmesher::fm_mesh_2d_inla(loc = data['centroid_lat', 'centroid_lon'], +mesh <- fmesher::fm_mesh_2d_inla(loc = data[, c("centroid_lon", "centroid_lat")], max.edge = c(3, 7.9), cutoff = 0.9, offset=c(0.2, 4)) -# plot(mesh) +plot(mesh) mesh ``` +This mesh now represents where the model will calculate the lake changes for the tibetan lakes region. Now that we have the mesh we can use the priors app to specify the model. + +```{r error=TRUE} +fdmr::model_builder(spatial_data = sp_data, measurement_data = covid19_data, mesh = mesh, time_variable = "week") +``` From da6b4b45465ac3e129a5c2d970307a22b108cbed Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Tue, 24 Oct 2023 13:19:15 +0200 Subject: [PATCH 12/24] model builder --- vignettes/tibetan-lakes.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 9d56e7ae..7b529047 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -128,6 +128,6 @@ mesh This mesh now represents where the model will calculate the lake changes for the tibetan lakes region. Now that we have the mesh we can use the priors app to specify the model. ```{r error=TRUE} -fdmr::model_builder(spatial_data = sp_data, measurement_data = covid19_data, mesh = mesh, time_variable = "week") +fdmr::model_builder(spatial_data = data, measurement_data = data, mesh = mesh, time_variable = "week") ``` From f110a0c5fca45c7eb3e765e20f7ec6ac643214bb Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Tue, 24 Oct 2023 16:01:34 +0200 Subject: [PATCH 13/24] added more for model builder --- vignettes/tibetan-lakes.Rmd | 33 +++++++++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 7b529047..15fdf24e 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -21,6 +21,7 @@ devtools::load_all() # TODO : move this data into the Rdataset repo data <- read.csv('~/repos/tblakes/data_for_4dm.csv') +# sp::coordinates(data) <- c("centroid_lon", "centroid_lat") summary(data) ``` @@ -122,12 +123,40 @@ mesh <- fmesher::fm_mesh_2d_inla(loc = data[, c("centroid_lon", "centroid_lat")] cutoff = 0.9, offset=c(0.2, 4)) plot(mesh) -mesh ``` This mesh now represents where the model will calculate the lake changes for the tibetan lakes region. Now that we have the mesh we can use the priors app to specify the model. ```{r error=TRUE} -fdmr::model_builder(spatial_data = data, measurement_data = data, mesh = mesh, time_variable = "week") +# codes for making stuff work + +library(dplyr) +# Define a function to calculate z-scores, except for specified columns +calculate_zscore <- function(df, group_col, ignore_cols) { + df %>% + group_by(across(all_of(group_col))) %>% + mutate(across( + where(is.numeric) & !all_of(ignore_cols), + list(zscore = ~ ( . - mean(., na.rm = TRUE)) / sd(., na.rm = TRUE)), + .names = "zscore_{col}" + )) %>% + ungroup() +} + +# Specify columns to ignore +ignore_cols <- c("centroid_lon", "centroid_lat", "year") + +# Apply the function to the data +data_with_zscores <- calculate_zscore(data, "HYBAS_ID", ignore_cols) +``` +moving on +```{r error=TRUE} +# needs spatial data +sp::coordinates(data) <- c("centroid_lon", "centroid_lat") +``` +moving on + +```{r error=TRUE} +fdmr::model_builder(spatial_data = data, measurement_data = data, mesh = mesh, time_variable = "year") ``` From 10197de3bc17b91177ca3682bb6a747c7105f9e2 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Tue, 24 Oct 2023 17:36:57 +0200 Subject: [PATCH 14/24] model working --- vignettes/tibetan-lakes.Rmd | 44 ++++++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 15fdf24e..bb7b7c0f 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -152,11 +152,49 @@ data_with_zscores <- calculate_zscore(data, "HYBAS_ID", ignore_cols) moving on ```{r error=TRUE} # needs spatial data -sp::coordinates(data) <- c("centroid_lon", "centroid_lat") +# sp::coordinates(data) <- c("centroid_lon", "centroid_lat") +sp::coordinates(data_with_zscores) <- c("centroid_lon", "centroid_lat") ``` -moving on +moving ongi ```{r error=TRUE} -fdmr::model_builder(spatial_data = data, measurement_data = data, mesh = mesh, time_variable = "year") +fdmr::model_builder(spatial_data = data_with_zscores, measurement_data = data_with_zscores, mesh = mesh, time_variable = "year") ``` +now we use the code from the model builder app to make the model: +```{r error=TRUE} +group_index <- data_with_zscores$year +n_groups <- length(unique(data_with_zscores$year)) +spde <- INLA::inla.spde2.pcmatern( + mesh = mesh, + prior.range = c(0.05,0.1), + prior.sigma = c(0.05,0.2) + ) + +alphaprior <- list(theta = list( + prior = 'pccor1', + param = c(-0.2,0.8) + )) + +formula <- zscore_water_balance_m3 ~ 0 + Intercept(1) + zscore_precip + f( + main = coordinates, + model = spde, + group = group_index, + ngroup = n_groups, + control.group = list( + model = 'ar1', + hyper = alphaprior) + ) + +model_output <- inlabru::bru(formula, + # data = measurement_data, + data = data_with_zscores, + family = 'gaussian', + E = NULL, + control.family = NULL, + options = list( + verbose = FALSE + ) + ) +summary(model_output) +``` From 433eab5bdc3304bcf8372b0adf8fcdca1da27309 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Tue, 24 Oct 2023 17:39:36 +0200 Subject: [PATCH 15/24] added model viewer --- vignettes/tibetan-lakes.Rmd | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index bb7b7c0f..ea001b75 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -198,3 +198,12 @@ model_output <- inlabru::bru(formula, ) summary(model_output) ``` + + +now we can evaluate the model... + +```{r error=TRUE} +fdmr::model_viewer(model_output = model_output, mesh = mesh, measurement_data = data_with_zscores, data_distribution = "Gaussian") +``` + + From 5cb04f6b9e6f3bab1d03e12a3be46d1d373cdafe Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Sun, 29 Oct 2023 14:08:57 +0100 Subject: [PATCH 16/24] histograms --- vignettes/tibetan-lakes.Rmd | 3 +++ 1 file changed, 3 insertions(+) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index ea001b75..b4498a4a 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -207,3 +207,6 @@ fdmr::model_viewer(model_output = model_output, mesh = mesh, measurement_data = ``` +```{r error=TRUE} +hist(model_output$summary.fitted.values$mean, breaks=seq(-0.3, 0.3, l=100)) +``` \ No newline at end of file From 049a8e681bec8f1843aea28904aac05b0fd559c0 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Tue, 31 Oct 2023 11:11:11 +0100 Subject: [PATCH 17/24] model eval --- vignettes/tibetan-lakes.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index b4498a4a..a5650f72 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -103,7 +103,7 @@ leaflet(summary_data) %>% Let's try and model this. We will use 4D-Modeller and inlabru. -[inlabru]() is designed to calculate fixed and random effects as well as continuous spatially and temporally distributed processes using SPDEs (). In order to make these computationally solveable, the continous processes are discretized on a finite element mesh. This mesh represents the spatial distribution of the process under study. At each node of the mesh, the model is used to calculate the outcome variable (in this case lake area changes). The mesh provides the spatial awareness to the model. Different processes can have different length scales, i.e., how far away does a process need to occur before it has no effect on the occurence of the process at my current location. These are controlled by priors (see below). +[inlabru](https://inlabru-org.github.io/inlabru/index.html) is designed to calculate fixed and random effects as well as continuous spatially and temporally distributed processes using SPDEs. In order to make these computationally solveable, the continous processes are discretized on a finite element mesh. This mesh represents the spatial distribution of the process under study. At each node of the mesh, the model is used to calculate the outcome variable (in this case lake area changes). The mesh provides the spatial awareness to the model. Different processes can have different length scales, i.e., how far away does a process need to occur before it has no effect on the occurence of the process at my current location. These are controlled by priors (see below). One major assumption in this method is that the process being modeled is continuous. This has the benefit that the spatial field being modeled can be calculated for regions where there is no data. 4D-Modeller is designed to make these model design decisions visual, interpretable, and accessible to users who may not have a strong background in R programming or bayesian spatiotemporal modeling. 4D-modeller has a set of shiny apps that help the user write the codes necessary to implement the full model and evaluate it's performance. From 3dde82448e226ccd7784def2a4358dd23f7d9891 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Tue, 31 Oct 2023 15:07:01 +0100 Subject: [PATCH 18/24] nothing works anymore --- vignettes/hydro.Rmd | 2 ++ vignettes/tibetan-lakes.Rmd | 46 ++++++++++++++++++++++++++++--------- 2 files changed, 37 insertions(+), 11 deletions(-) diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 4669cf78..2084e993 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -13,6 +13,8 @@ If you do not have `INLA` or `inlabru` installed please check our [installation First we load INLA and then download the data for this tutorial from our [data repository](https://github.com/4DModeller/fdmr_data) using `retrieve_tutorial_data`. ```{r} +library(devtools) +devtools::load_all() library(INLA) library(magrittr) ``` diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index a5650f72..99e5b1bd 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -176,17 +176,41 @@ alphaprior <- list(theta = list( param = c(-0.2,0.8) )) -formula <- zscore_water_balance_m3 ~ 0 + Intercept(1) + zscore_precip + f( - main = coordinates, - model = spde, - group = group_index, - ngroup = n_groups, - control.group = list( - model = 'ar1', - hyper = alphaprior) - ) - -model_output <- inlabru::bru(formula, +# formula <- zscore_water_balance_m3 ~ 0 + Intercept(1) + zscore_precip + f( +# main = coordinates, +# model = spde, +# group = group_index, +# ngroup = n_groups, +# control.group = list( +# model = 'ar1', +# hyper = alphaprior) +# ) + +# formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model='linear') +formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) + +# + zscore_precip + f( +# main = coordinates, +# model = spde, +# group = group_index, +# ngroup = n_groups, +# control.group = list( +# model = 'ar1', +# hyper = alphaprior) +# ) + +# model_output <- inlabru::bru(formula, +# # data = measurement_data, +# data = data_with_zscores, +# family = 'gaussian', +# E = NULL, +# control.family = NULL, +# options = list( +# verbose = FALSE +# ) +# ) + +model_output <- inlabru::bru(formula2, # data = measurement_data, data = data_with_zscores, family = 'gaussian', From e81cff6ff6ebad5b235bb6cfadc45d0981d68a38 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 1 Nov 2023 15:57:05 +0100 Subject: [PATCH 19/24] deleted old --- vignettes/tblakes_shiny.Rmd | 47 ------------------------------------- 1 file changed, 47 deletions(-) delete mode 100644 vignettes/tblakes_shiny.Rmd diff --git a/vignettes/tblakes_shiny.Rmd b/vignettes/tblakes_shiny.Rmd deleted file mode 100644 index bddf9da6..00000000 --- a/vignettes/tblakes_shiny.Rmd +++ /dev/null @@ -1,47 +0,0 @@ ---- -title: "Tibetan Lakes Shiny Example" -output: - bookdown::html_document2: - base_format: rmarkdown::html_vignette - fig_caption: yes -bibliography: "covid.bib" -link-citations: yes -vignette: > - %\VignetteIndexEntry{Tibetan Lakes Shiny Example} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -Here we will show how to use the shiny apps with the tibetan lakes dataset. - -```{r error=TRUE} -devtools::load_all() -library(INLA) -library(inlabru) -data <- read.csv('~/repos/tblakes/data_for_modeling.csv') -summary(data) -``` - - -now we will use the mesh builder app: - -```{r error=TRUE} -fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") -``` - -using the meshbuilder we have generated the following code: - -```{r error=TRUE} -mesh <- inla.mesh.2d(loc = data[c('LONG', 'LAT')], - max.edge = c(1.1, 8.1), - cutoff = 0.81, - offset=c(1.2, 4.6)) -plot(mesh) -``` - -next we will use the model tester app: - -```{r error=TRUE} -sp::coordinates(data) <- c('LONG', 'LAT') -fdmr::interactive_priors(data, measurement_data = data, time_variable = data$year) -``` \ No newline at end of file From 17bf3583a2730f52af2d6197e1abf94da114e84a Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 1 Nov 2023 16:23:06 +0100 Subject: [PATCH 20/24] added temp data --- vignettes/tibetan-lakes.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 99e5b1bd..228585d7 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -20,7 +20,7 @@ devtools::load_all() # fdmr::retrieve_tutorial_data(dataset = "tibetan lakes") # TODO : move this data into the Rdataset repo -data <- read.csv('~/repos/tblakes/data_for_4dm.csv') +data <- read.csv('~/repos/tblakes/data_for_4dm_w_temp.csv') # sp::coordinates(data) <- c("centroid_lon", "centroid_lat") summary(data) ``` From c51ea4b33f8225f4fb4838ccef17a97bb0022c23 Mon Sep 17 00:00:00 2001 From: mnky9800n Date: Wed, 1 Nov 2023 16:29:27 +0100 Subject: [PATCH 21/24] functioning with temp --- vignettes/tibetan-lakes.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 228585d7..35cebff2 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -186,8 +186,8 @@ alphaprior <- list(theta = list( # hyper = alphaprior) # ) -# formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model='linear') -formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) +formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model='linear') + temperature(data_with_zscores$zscore_t2m) + residuals(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) +# formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) # + zscore_precip + f( # main = coordinates, From 75eff7860d59d854dfa8f137b0eb5919efa8f1af Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Wed, 1 Nov 2023 19:17:02 +0000 Subject: [PATCH 22/24] Small changes to text and removal of imports --- vignettes/hydro.Rmd | 2 -- vignettes/tibetan-lakes.Rmd | 71 ++++++++----------------------------- 2 files changed, 14 insertions(+), 59 deletions(-) diff --git a/vignettes/hydro.Rmd b/vignettes/hydro.Rmd index 3e2c613e..f4b55f26 100644 --- a/vignettes/hydro.Rmd +++ b/vignettes/hydro.Rmd @@ -13,8 +13,6 @@ If you do not have `INLA` or `inlabru` installed please check our [installation First we load INLA and then download the data for this tutorial from our [data repository](https://github.com/4DModeller/fdmr_data) using `retrieve_tutorial_data`. ```{r} -library(devtools) -devtools::load_all() library(INLA) library(magrittr) ``` diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index 35cebff2..abf25646 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -15,12 +15,8 @@ vignette: > The Tibetan Plateau and the High Mountain Asia region provide water resources to approximately 1 billion people and provide for their health, economies, and agriculture. These natural "water towers" are highly sensitive to climate change and therefore changes in regional water resources must be understood (@immerzeel2020importance). Indeed, understanding how climate change will alter cold region water resources has been described as a "grand challenge" of hydrology (@twentythreeproblems). Recent changes in lakes in the Tibetan Plateau endorheic basin have led to increases in lake areas. Conventionally attributed to increases in glacier melt due to climate change, a recent investigation has demonstrated that glacier mass loss contributes only a small fraction (19%) to the lake volume increase (@brun2020limited). Thus, an unresolved question in the Tibetan Plateau endorheic basin lakes is, what physical processes drive lake growth in the region and at what fraction? ```{r error=TRUE} -library(devtools) -devtools::load_all() -# fdmr::retrieve_tutorial_data(dataset = "tibetan lakes") +fdmr::retrieve_tutorial_data(dataset = "tibetan_lakes") -# TODO : move this data into the Rdataset repo -data <- read.csv('~/repos/tblakes/data_for_4dm_w_temp.csv') # sp::coordinates(data) <- c("centroid_lon", "centroid_lat") summary(data) ``` @@ -101,35 +97,33 @@ leaflet(summary_data) %>% ) ``` -Let's try and model this. We will use 4D-Modeller and inlabru. +Let's try and model this. We will use 4DModeller and inlabru. [inlabru](https://inlabru-org.github.io/inlabru/index.html) is designed to calculate fixed and random effects as well as continuous spatially and temporally distributed processes using SPDEs. In order to make these computationally solveable, the continous processes are discretized on a finite element mesh. This mesh represents the spatial distribution of the process under study. At each node of the mesh, the model is used to calculate the outcome variable (in this case lake area changes). The mesh provides the spatial awareness to the model. Different processes can have different length scales, i.e., how far away does a process need to occur before it has no effect on the occurence of the process at my current location. These are controlled by priors (see below). One major assumption in this method is that the process being modeled is continuous. This has the benefit that the spatial field being modeled can be calculated for regions where there is no data. -4D-Modeller is designed to make these model design decisions visual, interpretable, and accessible to users who may not have a strong background in R programming or bayesian spatiotemporal modeling. 4D-modeller has a set of shiny apps that help the user write the codes necessary to implement the full model and evaluate it's performance. +4DModeller is designed to make these model design decisions visual, interpretable, and accessible to users who may not have a strong background in R programming or bayesian spatiotemporal modeling. 4DModeller has a set of shiny apps that help the user write the codes necessary to implement the full model and evaluate it's performance. First we will use the mesh_builder to build a reasonable mesh. -```{r error=TRUE} - +```{r error=TRUE, eval=FALSE} fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", longitude_column = "centroid_lon", latitude_column = "centroid_lat") ``` -Using the mesh_builder we have the code to build the mesh for the inlabru model. Note, we needed to replace `location data` with `data[, c('centroid_lon', 'centroid_lat')]`. +Using the mesh_builder we have the code to build the mesh for the inlabru model. ```{r error=TRUE} mesh <- fmesher::fm_mesh_2d_inla(loc = data[, c("centroid_lon", "centroid_lat")], max.edge = c(3, 7.9), cutoff = 0.9, offset=c(0.2, 4)) -plot(mesh) + +fdmr::plot_mesh(mesh) ``` This mesh now represents where the model will calculate the lake changes for the tibetan lakes region. Now that we have the mesh we can use the priors app to specify the model. ```{r error=TRUE} -# codes for making stuff work - library(dplyr) # Define a function to calculate z-scores, except for specified columns calculate_zscore <- function(df, group_col, ignore_cols) { @@ -148,19 +142,16 @@ ignore_cols <- c("centroid_lon", "centroid_lat", "year") # Apply the function to the data data_with_zscores <- calculate_zscore(data, "HYBAS_ID", ignore_cols) -``` -moving on -```{r error=TRUE} -# needs spatial data -# sp::coordinates(data) <- c("centroid_lon", "centroid_lat") sp::coordinates(data_with_zscores) <- c("centroid_lon", "centroid_lat") ``` -moving ongi -```{r error=TRUE} +Now we're ready to use `model_builder` to investigate how changing priors affects model predictions. + +```{r error=TRUE, eval=FALSE} fdmr::model_builder(spatial_data = data_with_zscores, measurement_data = data_with_zscores, mesh = mesh, time_variable = "year") ``` -now we use the code from the model builder app to make the model: + +Now we use the code from the model builder app to make the model: ```{r error=TRUE} group_index <- data_with_zscores$year @@ -176,39 +167,7 @@ alphaprior <- list(theta = list( param = c(-0.2,0.8) )) -# formula <- zscore_water_balance_m3 ~ 0 + Intercept(1) + zscore_precip + f( -# main = coordinates, -# model = spde, -# group = group_index, -# ngroup = n_groups, -# control.group = list( -# model = 'ar1', -# hyper = alphaprior) -# ) - formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model='linear') + temperature(data_with_zscores$zscore_t2m) + residuals(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) -# formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) - -# + zscore_precip + f( -# main = coordinates, -# model = spde, -# group = group_index, -# ngroup = n_groups, -# control.group = list( -# model = 'ar1', -# hyper = alphaprior) -# ) - -# model_output <- inlabru::bru(formula, -# # data = measurement_data, -# data = data_with_zscores, -# family = 'gaussian', -# E = NULL, -# control.family = NULL, -# options = list( -# verbose = FALSE -# ) -# ) model_output <- inlabru::bru(formula2, # data = measurement_data, @@ -223,14 +182,12 @@ model_output <- inlabru::bru(formula2, summary(model_output) ``` +now we can use the model viewer to view the output of the model and create some plots -now we can evaluate the model... - -```{r error=TRUE} +```{r error=TRUE, eval=FALSE} fdmr::model_viewer(model_output = model_output, mesh = mesh, measurement_data = data_with_zscores, data_distribution = "Gaussian") ``` - ```{r error=TRUE} hist(model_output$summary.fitted.values$mean, breaks=seq(-0.3, 0.3, l=100)) ``` \ No newline at end of file From 0113d72223af0ed38f0706d519c8e3f2cb282b9e Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Wed, 1 Nov 2023 19:17:28 +0000 Subject: [PATCH 23/24] Ran styler over tut --- vignettes/tibetan-lakes.Rmd | 62 ++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan-lakes.Rmd index abf25646..39f68369 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan-lakes.Rmd @@ -46,10 +46,10 @@ sums <- data %>% # TODO : combine these plots ggplot() + - geom_col(data = summary_data, aes(x = year, y = total, fill = sign), position = "stack") + geom_col(data = summary_data, aes(x = year, y = total, fill = sign), position = "stack") ggplot() + - geom_line(data = sums, aes(x = year, y = precip_sum * 100), color = 'cyan') # Scale + geom_line(data = sums, aes(x = year, y = precip_sum * 100), color = "cyan") # Scale ``` We can see here that while lakes both grow and decline, the overall trend is that lakes are growing more than they are shrinking in the Tibetan plateau. Why is this? perhaps it is because of the rain. Our ERA5-land data set gives us precipitation per year (after aggregation). But we see that the largest year that lakes grow, is not the largest amount of rain. What is going on here? @@ -63,7 +63,6 @@ ggplot(data, aes(x = precip, y = water_balance_m3, color = as.factor(year))) + Here we can see across the 20 year observation period a clear spatial distinction in lakegrowth however it is not as clear cut for precipitation. ```{r error=TRUE, fig.width=8, fig.height=8, fig.align='center'} - # TODO : this should include plotting the polygons of the catchments library(tidyverse) library(leaflet) @@ -82,7 +81,7 @@ leaflet(summary_data) %>% addCircles( lng = ~centroid_lon, lat = ~centroid_lat, - radius = ~50*total_balance, + radius = ~ 50 * total_balance, color = "blue", stroke = FALSE, fillOpacity = 0.6 @@ -90,7 +89,7 @@ leaflet(summary_data) %>% addCircles( lng = ~centroid_lon, lat = ~centroid_lat, - radius = ~500*total_precip, + radius = ~ 500 * total_precip, color = "red", stroke = FALSE, fillOpacity = 0.6 @@ -106,17 +105,18 @@ Let's try and model this. We will use 4DModeller and inlabru. First we will use the mesh_builder to build a reasonable mesh. ```{r error=TRUE, eval=FALSE} -fdmr::mesh_builder(spatial_data = data, obs_data = data, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", longitude_column = "centroid_lon", latitude_column = "centroid_lat") - +fdmr::mesh_builder(spatial_data = data, obs_data = data, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", longitude_column = "centroid_lon", latitude_column = "centroid_lat") ``` Using the mesh_builder we have the code to build the mesh for the inlabru model. ```{r error=TRUE} -mesh <- fmesher::fm_mesh_2d_inla(loc = data[, c("centroid_lon", "centroid_lat")], - max.edge = c(3, 7.9), - cutoff = 0.9, - offset=c(0.2, 4)) +mesh <- fmesher::fm_mesh_2d_inla( + loc = data[, c("centroid_lon", "centroid_lat")], + max.edge = c(3, 7.9), + cutoff = 0.9, + offset = c(0.2, 4) +) fdmr::plot_mesh(mesh) ``` @@ -131,7 +131,7 @@ calculate_zscore <- function(df, group_col, ignore_cols) { group_by(across(all_of(group_col))) %>% mutate(across( where(is.numeric) & !all_of(ignore_cols), - list(zscore = ~ ( . - mean(., na.rm = TRUE)) / sd(., na.rm = TRUE)), + list(zscore = ~ (. - mean(., na.rm = TRUE)) / sd(., na.rm = TRUE)), .names = "zscore_{col}" )) %>% ungroup() @@ -157,28 +157,28 @@ Now we use the code from the model builder app to make the model: group_index <- data_with_zscores$year n_groups <- length(unique(data_with_zscores$year)) spde <- INLA::inla.spde2.pcmatern( - mesh = mesh, - prior.range = c(0.05,0.1), - prior.sigma = c(0.05,0.2) - ) + mesh = mesh, + prior.range = c(0.05, 0.1), + prior.sigma = c(0.05, 0.2) +) alphaprior <- list(theta = list( - prior = 'pccor1', - param = c(-0.2,0.8) - )) + prior = "pccor1", + param = c(-0.2, 0.8) +)) -formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model='linear') + temperature(data_with_zscores$zscore_t2m) + residuals(data_with_zscores[, c('centroid_lon', 'centroid_lat')], model='spde', group=group_index, ngroup=n_groups, control.group=list(model='ar1', hyper=alphaprior)) +formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model = "linear") + temperature(data_with_zscores$zscore_t2m) + residuals(data_with_zscores[, c("centroid_lon", "centroid_lat")], model = "spde", group = group_index, ngroup = n_groups, control.group = list(model = "ar1", hyper = alphaprior)) model_output <- inlabru::bru(formula2, - # data = measurement_data, - data = data_with_zscores, - family = 'gaussian', - E = NULL, - control.family = NULL, - options = list( - verbose = FALSE - ) - ) + # data = measurement_data, + data = data_with_zscores, + family = "gaussian", + E = NULL, + control.family = NULL, + options = list( + verbose = FALSE + ) +) summary(model_output) ``` @@ -189,5 +189,5 @@ fdmr::model_viewer(model_output = model_output, mesh = mesh, measurement_data = ``` ```{r error=TRUE} -hist(model_output$summary.fitted.values$mean, breaks=seq(-0.3, 0.3, l=100)) -``` \ No newline at end of file +hist(model_output$summary.fitted.values$mean, breaks = seq(-0.3, 0.3, l = 100)) +``` From 94869350cfa41948851b40349428da2ba4dbd4fc Mon Sep 17 00:00:00 2001 From: Gareth Jones Date: Wed, 1 Nov 2023 21:01:13 +0000 Subject: [PATCH 24/24] Add retrieval of data from fdmr_data repo, add vignette to pkgdown list --- _pkgdown.yml | 2 ++ vignettes/{tibetan-lakes.Rmd => tibetan_lakes.Rmd} | 9 ++++++--- 2 files changed, 8 insertions(+), 3 deletions(-) rename vignettes/{tibetan-lakes.Rmd => tibetan_lakes.Rmd} (97%) diff --git a/_pkgdown.yml b/_pkgdown.yml index cfebb6a8..fce7dace 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -45,6 +45,8 @@ navbar: href: articles/hydro.html - text: "Priors exploration" href: articles/priors.html + - text: Tibetan Lakes + href: articles/tibetan_lakes.html - text: --- - text: "FAQ" - text: INLA Crash FAQ diff --git a/vignettes/tibetan-lakes.Rmd b/vignettes/tibetan_lakes.Rmd similarity index 97% rename from vignettes/tibetan-lakes.Rmd rename to vignettes/tibetan_lakes.Rmd index 39f68369..6b691d26 100644 --- a/vignettes/tibetan-lakes.Rmd +++ b/vignettes/tibetan_lakes.Rmd @@ -17,7 +17,9 @@ The Tibetan Plateau and the High Mountain Asia region provide water resources to ```{r error=TRUE} fdmr::retrieve_tutorial_data(dataset = "tibetan_lakes") -# sp::coordinates(data) <- c("centroid_lon", "centroid_lat") +data_filepath <- fdmr::get_tutorial_datapath(dataset = "tibetan_lakes", file = "tibetan_lakes_data.csv") +data <- read.csv(data_filepath) + summary(data) ``` @@ -64,7 +66,6 @@ ggplot(data, aes(x = precip, y = water_balance_m3, color = as.factor(year))) + Here we can see across the 20 year observation period a clear spatial distinction in lakegrowth however it is not as clear cut for precipitation. ```{r error=TRUE, fig.width=8, fig.height=8, fig.align='center'} # TODO : this should include plotting the polygons of the catchments -library(tidyverse) library(leaflet) # Summarize the data @@ -104,6 +105,8 @@ Let's try and model this. We will use 4DModeller and inlabru. First we will use the mesh_builder to build a reasonable mesh. +> **_NOTE:_** Due to the size of the mesh being created the initial load of the app may take some time + ```{r error=TRUE, eval=FALSE} fdmr::mesh_builder(spatial_data = data, obs_data = data, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", longitude_column = "centroid_lon", latitude_column = "centroid_lat") ``` @@ -170,7 +173,6 @@ alphaprior <- list(theta = list( formula2 <- zscore_water_balance_m3 ~ 0 + Intercept(1) + precipitation(data_with_zscores$zscore_precip, model = "linear") + temperature(data_with_zscores$zscore_t2m) + residuals(data_with_zscores[, c("centroid_lon", "centroid_lat")], model = "spde", group = group_index, ngroup = n_groups, control.group = list(model = "ar1", hyper = alphaprior)) model_output <- inlabru::bru(formula2, - # data = measurement_data, data = data_with_zscores, family = "gaussian", E = NULL, @@ -179,6 +181,7 @@ model_output <- inlabru::bru(formula2, verbose = FALSE ) ) + summary(model_output) ```