# INBO CODING CLUB 13 December 2022 Welcome! ## Share your code snippet If you want to share your code snippet, copy paste your snippet within a section of three backticks (```): As an **example**: ``` library(tidyverse) ``` (*you can copy paste this example and add your code further down*) ## Yellow sticky notes No yellow sticky notes online. Put your name + " | " and add a "*" each time you solve a challenge (see below). ## Participants Name | Challenges --- | --- Damiano Oldoni | *** Raïsa Carmen |** Lynn Pallemaerts |** Ward Langeraert |* Luc De Bruyn | * Emma Cartuyvels |** Dirk Maes | Hans Van Calster | ** Pieter Huybrechts |** Lynn Delgat |** Matthieu Chastel | ## Challenge 1 ### Damiano solution ### Ward solution ```` --- title: "Read and preprocess geese data" author: "Ward Langeraert" date: "`r Sys.Date()`" output: html_document: df_print: paged toc: true toc_float: true number_sections: true editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # CHALLENGE 1 Convert the code below to a Rmd document called `1_geese_read_data.Rmd` and make an html version of it. title : "Read and preprocess geese data" ## Setup Load libraries: ```{r, message=FALSE, warning=FALSE} library(tidyverse) # to do datascience library(geepack) # to do modelling library(INBOtheme) # to apply INBO style to graphs library(sf) # to work with geospatial vector data library(leaflet) # to make dynamic maps library(htmltools) # to make nice html labels for dynamic maps ``` ## Introduction In this document we will: 1. read geese data 2. explore data 3. preprocess data ## Read data Read catches and counts of geese in Flanders: ```{r} # credits to Hans for use here::here() catch_fl <- read_csv(here::here("data/20221213/20221213_geese_counts.txt"), na = "NA", col_types = cols( province = col_character(), location = col_character(), year = col_double(), latinName = col_character(), commonName = col_character(), counts = col_double(), catched = col_double(), adult = col_double(), pulli = col_double(), not_catched = col_double() )) ``` Number of catch data: ```{r} nrow(catch_fl) ``` Preview: ```{r} head(catch_fl, n = 30) ``` ## Explore data ### Taxonomic information Species present: ```{r} catch_fl %>% distinct(latinName, commonName) ``` ### Geographic information Data are geographically grouped by province and municipality (`location`): ```{r} catch_fl %>% distinct(province, location) ``` ### Temporal information The data are temporally defined at year level: ```{r} years <- catch_fl %>% distinct(year) %>% pull() ``` from: ```{r} min(years) ``` to: ```{r} max(years) ``` ## Preprocess data Data not linked to any `province` or `location` (`NAs`) will be removed. Number of rows removed: ```{r, echo=FALSE} catch_fl %>% filter(is.na(province) | is.na(location)) %>% nrow ``` Final dataset: ```{r, echo=FALSE} catch_fl <- catch_fl %>% filter(!is.na(province) & !is.na(location)) catch_fl ``` ```` ### Hans solution ```` --- title: "Read and preprocess geese data" author: "Hans Van Calster" date: "`r Sys.Date()`" output: html_document: toc: true toc_float: true number_sections: true df_print: paged --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) knitr::opts_knit$set(root.dir = here::here()) # Load libraries: library(tidyverse) # to do datascience library(INBOtheme) # to apply INBO style to graphs ``` # Introduction In this document we will: 1. read geese data 2. explore data 3. preprocess data # Read data Read catches and counts of geese in Flanders: ```{r read-geese-counts} catch_fl <- read_csv("data/20221213/20221213_geese_counts.txt", na = "NA", col_types = cols( province = col_character(), location = col_character(), year = col_double(), latinName = col_character(), commonName = col_character(), counts = col_double(), catched = col_double(), adult = col_double(), pulli = col_double(), not_catched = col_double() )) ``` Number of catch data: ```{r} nrow(catch_fl) ``` Preview: ```{r} head(catch_fl, n = 30) ``` # Explore data ## Taxonomic information Species present: ```{r taxa} catch_fl %>% distinct(latinName, commonName) ``` ## Geographic information Data are geographically grouped by province and municipality (`location`): ```{r locations} catch_fl %>% distinct(province, location) ``` ## Temporal information The data are temporally defined at year level: ```{r years} years <- catch_fl %>% distinct(year) %>% pull() ``` from: ```{r echo=FALSE} min(years) ``` to: ```{r echo=FALSE} max(years) ``` ## Preprocess data Data not linked to any `province` or `location` (`NAs`) will be removed. Number of rows removed: ```{r pre-process} catch_fl %>% filter(is.na(province) | is.na(location)) %>% nrow catch_fl <- catch_fl %>% filter(!is.na(province) & !is.na(location)) ``` Final dataset: ```{r filtered-data} catch_fl ``` ```` ### Pieter's solution ```` --- title: "Read and preprocess geese data" author: "PieterH" date: "`r format(Sys.time(), '%d %B, %Y')`" output: html_document: df_print: paged toc: true toc_depth: 4 toc_float: true number_sections: true collapsed: false smooth_scroll: true --- set `html_document:` to `df_print: paged` for pagination in dataframe outputs You can execute R code in the YAML header! You can set general options for all chunks here, in this case we want to echo by default! ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # Read and preprocess geese data ## Setup This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>. When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this: ```{r, message = FALSE} library(tidyverse) # to do datascience library(geepack) # to do modelling library(INBOtheme) # to apply INBO style to graphs library(sf) # to work with geospatial vector data library(leaflet) # to make dynamic maps library(htmltools) # to make nice html labels for dynamic maps ``` ## Introduction In this document we will: 1. read geese data 2. explore data 3. preprocess data You can also embed plots, for example: ## Read data Press `Ctrl+Alt+I` to insert a code chunk! ### Read catches and counts of geese in Flanders: ```{r} catch_fl <- read_csv(here::here("./data/20221213/20221213_geese_counts.txt"), na = "NA", col_types = cols( province = col_character(), location = col_character(), year = col_double(), latinName = col_character(), commonName = col_character(), counts = col_double(), catched = col_double(), adult = col_double(), pulli = col_double(), not_catched = col_double() )) ``` ### Number of catch data: ```{r} nrow(catch_fl) ``` #### Preview: ```{r} head(catch_fl, n = 30) ``` ## Explore data ### Taxonomic information #### Species present: ```{r} catch_fl %>% distinct(latinName, commonName) ``` ## Geographic information Data are geographically grouped by province and municipality (`location`): ```{r} catch_fl %>% distinct(province, location) ``` ## Temporal information The data are temporally defined at year level: ```{r} years <- catch_fl %>% distinct(year) %>% pull() ``` from: ```{r, echo = FALSE} min(years) ``` `echo = FALSE` results in the code not being shown, but it is evaluated to: ```{r} max(years) ``` ## Preprocess data Data not linked to any `province` or `location` (`NAs`) will be removed. Number of rows removed: ```{r, echo = FALSE} catch_fl %>% filter(is.na(province) | is.na(location)) %>% nrow ``` ### Final dataset: ```{r, echo = FALSE} catch_fl <- catch_fl %>% filter(!is.na(province) & !is.na(location)) catch_fl ``` ```` ## Challenge 2 ### Hans solution ```` --- title: "Data visualization" author: "Hans Van Calster" date: "`r Sys.Date()`" output: html_document: toc: true toc_float: true number_sections: true df_print: paged code_folding: hide --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE, tidy = "styler") knitr::opts_knit$set(root.dir = here::here()) library(tidyverse) # to do datascience library(geepack) # to do modelling library(INBOtheme) # to apply INBO style to graphs library(sf) # to work with geospatial vector data library(leaflet) # to make dynamic maps library(htmltools) # to make nice html labels for dynamic maps ``` In this section we will show how number of catches varies by year, province and species. Both static plots and dynamic maps are generated. # Read data Read catches and counts of geese in Flanders: ```{r read-geese-counts} catch_fl <- read_csv("data/20221213/20221213_geese_counts.txt", na = "NA", col_types = cols( province = col_character(), location = col_character(), year = col_double(), latinName = col_character(), commonName = col_character(), counts = col_double(), catched = col_double(), adult = col_double(), pulli = col_double(), not_catched = col_double() )) ``` Number of catch data: ```{r} nrow(catch_fl) ``` Preview: ```{r} head(catch_fl, n = 30) ``` # Explore data ## Taxonomic information Species present: ```{r taxa} catch_fl %>% distinct(latinName, commonName) ``` ## Geographic information Data are geographically grouped by province and municipality (`location`): ```{r locations} catch_fl %>% distinct(province, location) ``` ## Temporal information The data are temporally defined at year level: ```{r years} years <- catch_fl %>% distinct(year) %>% pull() ``` from: ```{r echo=FALSE} min(years) ``` to: ```{r echo=FALSE} max(years) ``` ## Preprocess data Data not linked to any `province` or `location` (`NAs`) will be removed. Number of rows removed: ```{r pre-process} catch_fl %>% filter(is.na(province) | is.na(location)) %>% nrow catch_fl <- catch_fl %>% filter(!is.na(province) & !is.na(location)) ``` Final dataset: ```{r filtered-data} catch_fl ``` # Static plots ## Show catches {.tabset} ### per province ```{r} catch_per_province <- catch_fl %>% group_by(province) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched_total)) ggplot(catch_per_province, aes(x = province, y = catched_total)) + geom_bar(stat = 'identity') + scale_x_discrete(breaks = 2009:2018) ``` ### per year ```{r} catch_per_year <- catch_fl %>% group_by(year) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched_total)) ggplot(catch_per_year, aes(x = year, y = catched_total)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ### per year and province ```{r} catch_per_year_province <- catch_fl %>% group_by(year, province) %>% summarize(catched = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched)) ggplot(catch_per_year_province, aes(x = year, y = catched, fill = province)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ## Catch analysis at species level ### Species selection Before we proceed to analyse the catches at species level, specify the species we are interested to by `commonName`: ```{r echo=TRUE} species <- c( "Brandgans", "Canadese gans", "Grauwe gans", "Soepgans", "Nijlgans" ) ``` ### Catches per year and species ```{r} catch_species <- catch_fl %>% filter(commonName %in% species) %>% group_by(year, commonName) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% arrange(commonName) ggplot(catch_species, aes(x = year, y = catched_total, fill = commonName)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ### Data modelling We apply a GEE (generalized estimating equations) model to data from 2010. ```{r} model_per_species <- map( species, function(s) { dfs <- catch_fl %>% filter(commonName == s & year >= 2010) %>% arrange(location, year) %>% mutate(year = as_factor(as.character(year)), location = as_factor(location)) geeglm(counts ~ 0 + year, family = poisson, data = dfs, waves = year, id = location) }) names(model_per_species) <- species overview_model <- map(model_per_species, ~summary(.)) overview_gee <- map2_dfr( overview_model, names(overview_model), function(model, name) { coefficients(model)[,1:2] %>% rownames_to_column(var = "year") %>% as_tibble() %>% mutate(species = name, year = str_sub(year, start = 5)) %>% select(species, everything()) }) overview_gee <- overview_gee %>% mutate( lwr = exp(Estimate - Std.err), upr = exp(Estimate + Std.err), Estimate = exp(Estimate) ) ``` ```{r} ggplot(overview_gee, aes(x = year, y = Estimate, ymin = lwr, ymax = upr)) + geom_errorbar(colour = "cyan3") + geom_point(colour = "cyan4") + facet_grid(.~species) + xlab("year") + ylab("Estimated number of geese per location") + theme_inbo(14) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 9)) + labs(title = "All provinces") ``` # Dynamic maps We make dynamic _leaflet_ maps of total number of catches per province. ```{r} pr_fl <- st_read("./data/20221213/20221213_flemish_provinces.gpkg", quiet = TRUE) pr_fl <- pr_fl %>% dplyr::left_join(catch_per_province, by = c("TX_PROV_DESCR_NL" = "province")) ``` ```{r} bins <- seq(0, 14000, by = 2000) pal <- colorBin("YlOrRd", domain = pr_fl$catched_total, bins = bins) labels <- sprintf( "<strong>%s</strong><br/>%g catches", pr_fl$TX_PROV_DESCR_NL, pr_fl$catched_total ) %>% lapply(HTML) map_catch_pr <- leaflet(pr_fl) %>% addTiles() %>% addPolygons( fillColor = ~pal(catched_total), weight = 2, opacity = 1, color = "white", dashArray = "3", fillOpacity = 0.7, highlightOptions = highlightOptions( weight = 5, color = "#666", dashArray = "", fillOpacity = 0.7, bringToFront = TRUE), label = labels, labelOptions = labelOptions( style = list("font-weight" = "normal", padding = "3px 8px"), textsize = "15px", direction = "auto")) map_catch_pr ``` ```` ### Pieter's solution ```` --- title: "Data visualization" author: "PieterH" date: "2022-12-13" output: html_document: code_folding: hide toc: true toc_depth: 4 toc_float: collapsed: false --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, cache.globals = TRUE) ``` We are dependent on 1_geese_read_data: ```{r} catch_fl <- data.table::fread("catch_fl.csv", encoding = "UTF-8") ``` Load libraries: ```{r class.source = 'fold-show', message = FALSE} library(tidyverse) # to do datascience library(geepack) # to do modelling library(INBOtheme) # to apply INBO style to graphs library(sf) # to work with geospatial vector data library(leaflet) # to make dynamic maps library(htmltools) # to make nice html labels for dynamic maps ``` # Data visualization In this section we will show how number of catches varies by year, province and species. Both static plots and dynamic maps are generated. ## Static plots ### Show catches {.tabset} #### per province ```{r} catch_per_province <- catch_fl %>% group_by(province) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched_total)) ggplot(catch_per_province, aes(x = province, y = catched_total)) + geom_bar(stat = 'identity') + scale_x_discrete(breaks = 2009:2018) ``` #### per year ```{r} catch_per_year <- catch_fl %>% group_by(year) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched_total)) ggplot(catch_per_year, aes(x = year, y = catched_total)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` #### per year and province ```{r} catch_per_year_province <- catch_fl %>% group_by(year, province) %>% summarize(catched = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched)) ggplot(catch_per_year_province, aes(x = year, y = catched, fill = province)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ### Catch analysis at species level #### Species selection Before we proceed to analyse the catches at species level, specify the species we are interested to by `commonName`: SHOW THIS CHUNK CODE! ```{r class.source = 'fold-show'} species <- c( "Brandgans", "Canadese gans", "Grauwe gans", "Soepgans", "Nijlgans" ) ``` #### Catches per year and species ```{r message = FALSE} catch_species <- catch_fl %>% filter(commonName %in% species) %>% group_by(year, commonName) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% arrange(commonName) ggplot(catch_species, aes(x = year, y = catched_total, fill = commonName)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` #### Data modelling We apply a GEE (generalized estimating equations) model to data from 2010. ```{r} model_per_species <- map( species, function(s) { dfs <- catch_fl %>% filter(commonName == s & year >= 2010) %>% arrange(location, year) %>% mutate(year = as_factor(as.character(year)), location = as_factor(location)) geeglm(counts ~ 0 + year, family = poisson, data = dfs, waves = year, id = location) }) names(model_per_species) <- species overview_model <- map(model_per_species, ~summary(.)) overview_gee <- map2_dfr( overview_model, names(overview_model), function(model, name) { coefficients(model)[,1:2] %>% rownames_to_column(var = "year") %>% as_tibble() %>% mutate(species = name, year = str_sub(year, start = 5)) %>% select(species, everything()) }) overview_gee <- overview_gee %>% mutate( lwr = exp(Estimate - Std.err), upr = exp(Estimate + Std.err), Estimate = exp(Estimate) ) ggplot(overview_gee, aes(x = year, y = Estimate, ymin = lwr, ymax = upr)) + geom_errorbar(colour = "cyan3") + geom_point(colour = "cyan4") + facet_grid(.~species) + xlab("year") + ylab("Estimated number of geese per location") + theme_inbo(14) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 9)) + labs(title = "All provinces") ``` ## Dynamic maps We make dynamic _leaflet_ maps of total number of catches per province. ```{r} pr_fl <- st_read(here::here("./data/20221213/20221213_flemish_provinces.gpkg")) pr_fl <- pr_fl %>% dplyr::left_join(catch_per_province, by = c("TX_PROV_DESCR_NL" = "province")) bins <- seq(0, 14000, by = 2000) pal <- colorBin("YlOrRd", domain = pr_fl$catched_total, bins = bins) labels <- sprintf( "<strong>%s</strong><br/>%g catches", pr_fl$TX_PROV_DESCR_NL, pr_fl$catched_total ) %>% lapply(HTML) map_catch_pr <- leaflet(pr_fl) %>% addTiles() %>% addPolygons( fillColor = ~pal(catched_total), weight = 2, opacity = 1, color = "white", dashArray = "3", fillOpacity = 0.7, highlightOptions = highlightOptions( weight = 5, color = "#666", dashArray = "", fillOpacity = 0.7, bringToFront = TRUE), label = labels, labelOptions = labelOptions( style = list("font-weight" = "normal", padding = "3px 8px"), textsize = "15px", direction = "auto")) map_catch_pr ``` ```` ### Solution Ward ```` --- title: "Data visualization" author: "Ward Langeraert" date: "2022-12-13" output: html_document: df_print: paged toc: true toc_float: collapsed: FALSE number_sections: true code_folding: show editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # CHALLENGE 2 Title: Data visualization In this section we will show how number of catches varies by year, province and species. Both static plots and dynamic maps are generated. ## Setup Load libraries: ```{r, message=FALSE, warning=FALSE} library(tidyverse) # to do datascience library(geepack) # to do modelling library(INBOtheme) # to apply INBO style to graphs library(sf) # to work with geospatial vector data library(leaflet) # to make dynamic maps library(htmltools) # to make nice html labels for dynamic maps ``` ## Read data Read catches and counts of geese in Flanders: ```{r} catch_fl <- read.csv(here::here("data/20221213/20221213_geese_counts_processed.txt")) ``` ## Static plots {.tabset .tabset-fade .tabset-pills} Show catches (make a tabbed section out of it) ### per province (1st tab) ```{r} catch_per_province <- catch_fl %>% group_by(province) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched_total)) ggplot(catch_per_province, aes(x = province, y = catched_total)) + geom_bar(stat = 'identity') + scale_x_discrete(breaks = 2009:2018) ``` ### per year (2nd tab) ```{r} catch_per_year <- catch_fl %>% group_by(year) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched_total)) ggplot(catch_per_year, aes(x = year, y = catched_total)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ### per year and province (3rd tab) ```{r} catch_per_year_province <- catch_fl %>% group_by(year, province) %>% summarize(catched = sum(catched, na.rm = TRUE)) %>% ungroup() %>% arrange(desc(catched)) ggplot(catch_per_year_province, aes(x = year, y = catched, fill = province)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ## Catch analysis at species level ### Species selection Before we proceed to analyse the catches at species level, specify the species we are interested to by `commonName`: ```{r} species <- c( "Brandgans", "Canadese gans", "Grauwe gans", "Soepgans", "Nijlgans" ) ``` ### Catches per year and species ```{r} catch_species <- catch_fl %>% filter(commonName %in% species) %>% group_by(year, commonName) %>% summarize(catched_total = sum(catched, na.rm = TRUE)) %>% arrange(commonName) ggplot(catch_species, aes(x = year, y = catched_total, fill = commonName)) + geom_bar(stat = 'identity') + scale_x_continuous(breaks = 2009:2018) ``` ### Data modelling We apply a GEE (generalized estimating equations) model to data from 2010. ```{r, class.source = 'fold-hide'} model_per_species <- map( species, function(s) { dfs <- catch_fl %>% filter(commonName == s & year >= 2010) %>% arrange(location, year) %>% mutate(year = as_factor(as.character(year)), location = as_factor(location)) geeglm(counts ~ 0 + year, family = poisson, data = dfs, waves = year, id = location) }) names(model_per_species) <- species overview_model <- map(model_per_species, ~summary(.)) overview_gee <- map2_dfr( overview_model, names(overview_model), function(model, name) { coefficients(model)[,1:2] %>% rownames_to_column(var = "year") %>% as_tibble() %>% mutate(species = name, year = str_sub(year, start = 5)) %>% select(species, everything()) }) overview_gee <- overview_gee %>% mutate( lwr = exp(Estimate - Std.err), upr = exp(Estimate + Std.err), Estimate = exp(Estimate)) ggplot(overview_gee, aes(x = year, y = Estimate, ymin = lwr, ymax = upr)) + geom_errorbar(colour = "cyan3") + geom_point(colour = "cyan4") + facet_grid(.~species) + xlab("year") + ylab("Estimated number of geese per location") + theme_inbo(14) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 9)) + labs(title = "All provinces") ``` ## Dynamic maps We make dynamic _leaflet_ maps of total number of catches per province. ```{r, class.source = 'fold-hide'} pr_fl <- st_read(here::here("data/20221213/20221213_flemish_provinces.gpkg")) pr_fl <- pr_fl %>% dplyr::left_join(catch_per_province, by = c("TX_PROV_DESCR_NL" = "province")) bins <- seq(0, 14000, by = 2000) pal <- colorBin("YlOrRd", domain = pr_fl$catched_total, bins = bins) labels <- sprintf( "<strong>%s</strong><br/>%g catches", pr_fl$TX_PROV_DESCR_NL, pr_fl$catched_total ) %>% lapply(HTML) map_catch_pr <- leaflet(pr_fl) %>% addTiles() %>% addPolygons( fillColor = ~pal(catched_total), weight = 2, opacity = 1, color = "white", dashArray = "3", fillOpacity = 0.7, highlightOptions = highlightOptions( weight = 5, color = "#666", dashArray = "", fillOpacity = 0.7, bringToFront = TRUE), label = labels, labelOptions = labelOptions( style = list("font-weight" = "normal", padding = "3px 8px"), textsize = "15px", direction = "auto")) map_catch_pr ``` ```` ## Challenge 3