diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..4a9d8e8
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,2 @@
+*.png filter=lfs diff=lfs merge=lfs -text
+*.gif filter=lfs diff=lfs merge=lfs -text
diff --git a/case_study_3/.gitignore b/case_study_3/.gitignore
new file mode 100644
index 0000000..5870e0b
--- /dev/null
+++ b/case_study_3/.gitignore
@@ -0,0 +1,2 @@
+case_study_3_reproducible_cache
+case_study_3_additional_cache
diff --git a/case_study_3/case_study_3_additional.Rmd b/case_study_3/case_study_3_additional.Rmd
new file mode 100644
index 0000000..40c7444
--- /dev/null
+++ b/case_study_3/case_study_3_additional.Rmd
@@ -0,0 +1,161 @@
+---
+title: "Additional exploration"
+output: github_document
+---
+
+<!---# Source - https://stackoverflow.com/a
+# Posted by Thomas
+# Retrieved 2025-12-09, License - CC BY-SA 3.0 -->
+
+
+First some code from the reproducible example is repeated to use the same time frame
+
+```{r echo=FALSE}
+knitr::opts_chunk$set(cache = TRUE)
+invisible(knitr::purl("case_study_3_reproducable.Rmd", output="temp.R", quiet=TRUE))
+knitr::read_chunk("temp.R")
+```
+
+```{r ref.label='requiredPkg'}
+```
+```{r ref.label='datetime'}
+```
+
+```{r echo=FALSE}
+unlink("temp.R")
+```
+Lets take a bit wider interval to see the context
+
+
+```{r interv}
+small_scale_radar_interval<-lubridate::interval(time-hours(10), time+hours(10))
+```
+## Robin Radars located at sea
+
+To investigate the high densities at sea we analyse two robin radars located there. 
+
+```{r uvaRR}
+require(uvaRR)
+trks <- purrr::map(
+  c("rws01", "borssele"),
+  ~ rr_track(
+    database = .x,
+    ts =c(int_start(small_scale_radar_interval), int_end(small_scale_radar_interval)) ,
+    interv = 60 * 60,
+    geom = T,
+    classification = c("BIRD"),
+    movement_metrics = c("dir", "dist", "edist", "spd", "str", "z")
+  ) |>
+    tibble::add_column(radar = .x)
+) |>
+  bind_rows() |>
+  mutate(trajectory = sf::st_as_sfc(trajectory, crs = 4326)) |>
+  st_as_sf() |>
+  mutate(
+    radarid_unique = trajectory_radarid |>
+      purrr::map(
+        ~ unlist(purrr::map(
+          purrr::map(strsplit(gsub("\\{|\\}", "", .x), ","), as.numeric),
+          unique
+        ))
+      ),
+    radarid_unique_chr = purrr::map_chr(
+      purrr::map(radarid_unique, as.character),
+      paste,
+      collapse = " "
+    ),
+    radar_type = case_match(
+      radarid_unique_chr,
+      c("11 14", "5 8") ~ "mixed",
+      c("14", "8") ~ "vertical",
+      c("11", "5") ~ "horizontal"
+    )
+  )
+```
+There is an increased density of birds in Luchterduinen (`rws01`) while this is absent from Borssele.
+```{r rrplot}
+ggplot(trks) +
+  geom_histogram(aes(x = timestamp_start, fill = classification)) +
+  facet_grid(radar ~ .) +
+  geom_vline(xintercept = time)
+ggplot(trks) +
+  geom_histogram(aes(x = timestamp_start, fill = radar_type)) +
+  facet_grid(radar ~ .) +
+  geom_vline(xintercept = time)
+```
+A directional histogram seems to confirm that the movements are vey directional corresponding to migratory movements
+```{r dirhist}
+ggplot(
+  trks |>
+    dplyr::filter(
+      radar_type != "vertical",
+      timestamp_start <
+        lubridate::ceiling_date(time + lubridate::hours(2), "hours"),
+      timestamp_start >
+        lubridate::floor_date(time - lubridate::hours(2), "hours")
+    )
+) +
+  geom_histogram(aes(x = dir, fill = radar_type), binwidth = 30) +
+  facet_grid(radar ~ lubridate::floor_date(timestamp_start, "hour")) +
+  coord_polar() +
+  scale_x_continuous(limits = c(0, 360), breaks = 0:3*90)
+```
+
+
+## MR1 Radars
+
+
+```{r mr1}
+bs_radar_databases = c("sbrs_kid_bs291", "sbrs_ams_bs290")
+bs_data <- purrr::map(
+  bs_radar_databases, ~
+    DBI::dbGetQuery(
+      as.is = F,
+      uvaAuth::get_connection_postgresql(.x),
+      glue::glue(
+        "select * from collection left join rf_classification on rf_classification.echo = collection.row
+      left join rfclasses on rfclasses.id = rf_classification.class  where time_stamp between '{format(lubridate::int_start(small_scale_radar_interval))}' and '{format(lubridate::int_end(small_scale_radar_interval))}';
+"
+      )
+    ) |>
+      dplyr::mutate(
+        speed = feature3,
+        azimuth = feature2,
+        altitude_AGL = feature1
+      ) %>%
+      dplyr::select(-starts_with("feature")) %>%
+      tibble::add_column(database = .x)
+) |>
+  bind_rows()
+
+```
+
+```{r explore}
+bs_data_filtered<-bs_data |> dplyr::filter(!name%in%c('insect', 'undefined', 'nonbio','precipitation'))
+table(bs_data$name)
+
+ggplot(bs_data_filtered) +
+  geom_histogram(aes(x = time_stamp, fill = name)) +
+  facet_grid(database ~ .) +
+  geom_vline(xintercept = time)
+```
+
+Here we see that the radar in Kijkduin has a migratory peak in the morning while that is absent from Amsterdam.
+
+Lets explore the directions in Kijkduin:
+```{r dirhistBs}
+ggplot(
+  bs_data_filtered |>
+    dplyr::filter(
+     database=='sbrs_kid_bs291',
+      time_stamp <
+        lubridate::ceiling_date(time + lubridate::hours(2), "hours"),
+      time_stamp >
+        lubridate::floor_date(time - lubridate::hours(2), "hours")
+    )
+) +
+  geom_histogram(aes(x = azimuth, fill = name), binwidth = 15) +
+  facet_wrap(. ~ lubridate::floor_date(time_stamp, "hour")) +
+  coord_polar() +
+  scale_x_continuous(limits = c(0, 360), breaks = 0:3*90)
+```
diff --git a/case_study_3/case_study_3_additional.md b/case_study_3/case_study_3_additional.md
new file mode 100644
index 0000000..9ffa171
--- /dev/null
+++ b/case_study_3/case_study_3_additional.md
@@ -0,0 +1,235 @@
+Additional exploration
+================
+
+<!---# Source - https://stackoverflow.com/a
+# Posted by Thomas
+# Retrieved 2025-12-09, License - CC BY-SA 3.0 -->
+
+First some code from the reproducible example is repeated to use the
+same time frame
+
+``` r
+require(bioRad)
+require(getRad)
+require(lubridate)
+require(dplyr)
+require(sf)
+require(ggplot2)
+require(stars)
+```
+
+``` r
+time <- as_datetime("2022-10-19 8:00:00")
+radarCode <- c("nldhl", "nlhrw")
+radarLocation <- get_weather_radars() |>
+  dplyr::filter(radar %in% radarCode) |>
+  st_transform(crs = 28992)
+radarGrid <- raster::raster(st_buffer(radarLocation, units::set_units(105, "km")),
+  res = 500
+)
+```
+
+Lets take a bit wider interval to see the context
+
+``` r
+small_scale_radar_interval<-lubridate::interval(time-hours(10), time+hours(10))
+```
+
+## Robin Radars located at sea
+
+To investigate the high densities at sea we analyse two robin radars
+located there.
+
+``` r
+require(uvaRR)
+trks <- purrr::map(
+  c("rws01", "borssele"),
+  ~ rr_track(
+    database = .x,
+    ts =c(int_start(small_scale_radar_interval), int_end(small_scale_radar_interval)) ,
+    interv = 60 * 60,
+    geom = T,
+    classification = c("BIRD"),
+    movement_metrics = c("dir", "dist", "edist", "spd", "str", "z")
+  ) |>
+    tibble::add_column(radar = .x)
+) |>
+  bind_rows() |>
+  mutate(trajectory = sf::st_as_sfc(trajectory, crs = 4326)) |>
+  st_as_sf() |>
+  mutate(
+    radarid_unique = trajectory_radarid |>
+      purrr::map(
+        ~ unlist(purrr::map(
+          purrr::map(strsplit(gsub("\\{|\\}", "", .x), ","), as.numeric),
+          unique
+        ))
+      ),
+    radarid_unique_chr = purrr::map_chr(
+      purrr::map(radarid_unique, as.character),
+      paste,
+      collapse = " "
+    ),
+    radar_type = case_match(
+      radarid_unique_chr,
+      c("11 14", "5 8") ~ "mixed",
+      c("14", "8") ~ "vertical",
+      c("11", "5") ~ "horizontal"
+    )
+  )
+```
+
+    ## no assignable properties filter specified in WHERE-clause
+
+    ## classification filter included in WHERE-clause
+
+    ## no tracktype filter (RaAz, RaEl, RaAzEl) specified in WHERE-clause
+
+    ## nr_of_plots filter included in WHERE-clause is less than 1500
+
+    ## no polygon filter specified in WHERE-clause
+
+    ## start while loop
+
+    ## 0h 1 min 36 sec
+
+    ## no assignable properties filter specified in WHERE-clause
+
+    ## classification filter included in WHERE-clause
+
+    ## no tracktype filter (RaAz, RaEl, RaAzEl) specified in WHERE-clause
+
+    ## nr_of_plots filter included in WHERE-clause is less than 1500
+
+    ## no polygon filter specified in WHERE-clause
+
+    ## start while loop
+
+    ## 0h 0 min 33 sec
+
+There is an increased density of birds in Luchterduinen (`rws01`) while
+this is absent from Borssele.
+
+``` r
+ggplot(trks) +
+  geom_histogram(aes(x = timestamp_start, fill = classification)) +
+  facet_grid(radar ~ .) +
+  geom_vline(xintercept = time)
+```
+
+    ## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
+
+![](case_study_3_additional_files/figure-gfm/rrplot-1.png)<!-- -->
+
+``` r
+ggplot(trks) +
+  geom_histogram(aes(x = timestamp_start, fill = radar_type)) +
+  facet_grid(radar ~ .) +
+  geom_vline(xintercept = time)
+```
+
+    ## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
+
+![](case_study_3_additional_files/figure-gfm/rrplot-2.png)<!-- --> A
+directional histogram seems to confirm that the movements are vey
+directional corresponding to migratory movements
+
+``` r
+ggplot(
+  trks |>
+    dplyr::filter(
+      radar_type != "vertical",
+      timestamp_start <
+        lubridate::ceiling_date(time + lubridate::hours(2), "hours"),
+      timestamp_start >
+        lubridate::floor_date(time - lubridate::hours(2), "hours")
+    )
+) +
+  geom_histogram(aes(x = dir, fill = radar_type), binwidth = 30) +
+  facet_grid(radar ~ lubridate::floor_date(timestamp_start, "hour")) +
+  coord_polar() +
+  scale_x_continuous(limits = c(0, 360), breaks = 0:3*90)
+```
+
+    ## Warning: Removed 32 rows containing missing values or values outside the scale range
+    ## (`geom_bar()`).
+
+![](case_study_3_additional_files/figure-gfm/dirhist-1.png)<!-- -->
+
+## MR1 Radars
+
+``` r
+bs_radar_databases = c("sbrs_kid_bs291", "sbrs_ams_bs290")
+bs_data <- purrr::map(
+  bs_radar_databases, ~
+    DBI::dbGetQuery(
+      as.is = F,
+      uvaAuth::get_connection_postgresql(.x),
+      glue::glue(
+        "select * from collection left join rf_classification on rf_classification.echo = collection.row
+      left join rfclasses on rfclasses.id = rf_classification.class  where time_stamp between '{format(lubridate::int_start(small_scale_radar_interval))}' and '{format(lubridate::int_end(small_scale_radar_interval))}';
+"
+      )
+    ) |>
+      dplyr::mutate(
+        speed = feature3,
+        azimuth = feature2,
+        altitude_AGL = feature1
+      ) %>%
+      dplyr::select(-starts_with("feature")) %>%
+      tibble::add_column(database = .x)
+) |>
+  bind_rows()
+```
+
+``` r
+bs_data_filtered<-bs_data |> dplyr::filter(!name%in%c('insect', 'undefined', 'nonbio','precipitation'))
+table(bs_data$name)
+```
+
+    ## 
+    ##     bird_flock         insect     large_bird         nonbio passerine_type 
+    ##           1336           8781             99           1872           9903 
+    ##  precipitation     swift_type      undefined      unid_bird     wader_type 
+    ##            178            536              1           5295           1623
+
+``` r
+ggplot(bs_data_filtered) +
+  geom_histogram(aes(x = time_stamp, fill = name)) +
+  facet_grid(database ~ .) +
+  geom_vline(xintercept = time)
+```
+
+    ## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
+
+![](case_study_3_additional_files/figure-gfm/explore-1.png)<!-- -->
+
+Here we see that the radar in Kijkduin has a migratory peak in the
+morning while that is absent from Amsterdam.
+
+Lets explore the directions in Kijkduin:
+
+``` r
+ggplot(
+  bs_data_filtered |>
+    dplyr::filter(
+     database=='sbrs_kid_bs291',
+      time_stamp <
+        lubridate::ceiling_date(time + lubridate::hours(2), "hours"),
+      time_stamp >
+        lubridate::floor_date(time - lubridate::hours(2), "hours")
+    )
+) +
+  geom_histogram(aes(x = azimuth, fill = name), binwidth = 15) +
+  facet_wrap(. ~ lubridate::floor_date(time_stamp, "hour")) +
+  coord_polar() +
+  scale_x_continuous(limits = c(0, 360), breaks = 0:3*90)
+```
+
+    ## Warning: Removed 1487 rows containing non-finite outside the scale range
+    ## (`stat_bin()`).
+
+    ## Warning: Removed 44 rows containing missing values or values outside the scale range
+    ## (`geom_bar()`).
+
+![](case_study_3_additional_files/figure-gfm/dirhistBs-1.png)<!-- -->
diff --git a/case_study_3/case_study_3_additional_files/figure-gfm/dirhist-1.png b/case_study_3/case_study_3_additional_files/figure-gfm/dirhist-1.png
new file mode 100644
index 0000000..7cc32cb
--- /dev/null
+++ b/case_study_3/case_study_3_additional_files/figure-gfm/dirhist-1.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5d6867fb95b7954f8c534af37871a5bb92809a66f71faac3dfc8d13be22b22cc
+size 54072
diff --git a/case_study_3/case_study_3_additional_files/figure-gfm/dirhistBs-1.png b/case_study_3/case_study_3_additional_files/figure-gfm/dirhistBs-1.png
new file mode 100644
index 0000000..9fe4d79
--- /dev/null
+++ b/case_study_3/case_study_3_additional_files/figure-gfm/dirhistBs-1.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fc7f0fcf575807def86624eec31cbb04646674214bf946235a23d504f2981a8f
+size 60052
diff --git a/case_study_3/case_study_3_additional_files/figure-gfm/explore-1.png b/case_study_3/case_study_3_additional_files/figure-gfm/explore-1.png
new file mode 100644
index 0000000..6d747b8
--- /dev/null
+++ b/case_study_3/case_study_3_additional_files/figure-gfm/explore-1.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b9285de5e4811fb309a12d28dd75c1e779dd68d6288fdc12de12939ec091117c
+size 28539
diff --git a/case_study_3/case_study_3_additional_files/figure-gfm/rrplot-1.png b/case_study_3/case_study_3_additional_files/figure-gfm/rrplot-1.png
new file mode 100644
index 0000000..632e5e1
--- /dev/null
+++ b/case_study_3/case_study_3_additional_files/figure-gfm/rrplot-1.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0e312b2597739a853a98bdedde364e215f7857b2cbb33b4a574042e5f52dc731
+size 20656
diff --git a/case_study_3/case_study_3_additional_files/figure-gfm/rrplot-2.png b/case_study_3/case_study_3_additional_files/figure-gfm/rrplot-2.png
new file mode 100644
index 0000000..7d62000
--- /dev/null
+++ b/case_study_3/case_study_3_additional_files/figure-gfm/rrplot-2.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4e7a855502c9f8cb01f7fea68f30fcba36fcacb2e3d7632cd9c23b5b0611e8a5
+size 18789
diff --git a/case_study_3/case_study_3_reproducible.Rmd b/case_study_3/case_study_3_reproducible.Rmd
new file mode 100644
index 0000000..6db04fc
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible.Rmd
@@ -0,0 +1,271 @@
+---
+title: "Case Study weather radar"
+output: github_document
+---
+
+```{r setup, echo=FALSE}
+knitr::opts_chunk$set(cache = TRUE)
+```
+
+
+
+## Setup
+
+
+
+```{r requiredPkg, cache=FALSE, message=FALSE}
+require(bioRad)
+require(getRad)
+require(lubridate)
+require(dplyr)
+require(sf)
+require(ggplot2)
+require(stars)
+```
+
+In this case study we zoom in to the morning of the 19th of October 2022. We select both radar in the Netherlands. For a quick scan of the data the crow visualizer [aloft](https://crow.aloftdata.eu#/?radar=bejab&date=2022-10-19&interval=24&timedisplay=radarLocal&vpiMode=mtr&vpColorScheme=viridis&lang=en) repository can be consulted.
+
+```{r datetime}
+time <- as_datetime("2022-10-19 8:00:00")
+radarCode <- c("nldhl", "nlhrw")
+radarLocation <- get_weather_radars() |>
+  dplyr::filter(radar %in% radarCode) |>
+  st_transform(crs = 28992)
+radarGrid <- raster::raster(st_buffer(radarLocation, units::set_units(105, "km")),
+  res = 500
+)
+```
+
+On this night high densities of numbers were counted along the dutch coast, as seen in the following video on youtube.
+
+[![Watch the video](https://img.youtube.com/vi/8x4Dos7GbvI/default.jpg)](https://youtu.be/8x4Dos7GbvI)
+
+
+
+For visualizations we load a polygon map, in this case from the `cbsodataR` package. This is a slightly nicer map compared to `rnaturalearth` which globally provides reasonable maps.
+```{r map}
+require(cbsodataR)
+map <- cbs_get_sf("provincie", year = 2022)
+```
+
+
+
+# Mornings activity
+
+To explore the mornings activbity we can use the vpts. 
+
+```{r vpts}
+vpts <- get_vpts(radar = radarCode, as.Date(time))
+plot(vpts$nlhrw, main="Herwijnen")
+plot(vpts$nldhl, main="Den Helder")
+```
+
+
+# Static map
+
+The first step is to download radar data and calculate vertical distributions from the data
+```{r examplePlot}
+# Set names is to ensure the resulting list is named according to the radar code
+pvol_list <- setNames(
+  get_pvol(
+    radarCode, time,
+    elev_max = 89
+  ), radarCode
+)
+vp_list <- purrr::map(pvol_list, calculate_vp)
+```
+
+The next step is to create for each radar separately a vertically integrated estimate of density. Here we convert all cases where strong corrections have been applied to `NA`
+
+```{r createStars}
+ppi_strs <- purrr::map2(pvol_list, vp_list, integrate_to_ppi, raster = radarGrid) |>
+  purrr::map(purrr::pluck, "data") |>
+  purrr::map(stars::st_as_stars) |>
+  c(list(along = "radar")) |>
+  do.call(what = "c") |>
+  mutate(VID = dplyr::if_else(R > 10, NA, VID))
+```
+
+To combine both radars we create a dataset where we identify which radar on which location as been minimally corrected. 
+
+
+
+```{r exampleR}
+minR <- st_apply(select(ppi_strs, R), 1:2, function(x) {
+  if (all(is.na(x))) {
+    x
+  } else {
+    x == min(x, na.rm = T)
+  }
+}, keep = T) |>
+  setNames("minR") |>
+  aperm(c(2, 3, 1))
+```
+By  multiplying the VID with the minimal R value (0/1)and summing the values for each location we can get the best VID estimate for each location
+```{r mul}
+ppi_minR <- c(ppi_strs, minR) |>
+  mutate(VID = VID * minR) |>
+  st_apply(c("x", "y"), sum, na.rm = T)
+```
+
+```{r exampleRplot}
+#| fig-height: 9
+ggplot() +
+  geom_stars(data = ppi_minR, aes(fill = VID), color = NA, sf = T) +
+  scale_fill_viridis_c(
+    name = "Bird density [1/km^2]",
+    option = "turbo",
+    limit = c(2, 800),
+    transform = "log10"
+  ) +
+  geom_sf(data = map, fill = NA, color = "black", size = .3) +
+  geom_sf(data=radarLocation)+
+  ggspatial::annotation_scale() +
+  coord_sf(
+    expand = F,
+    xlim = st_bbox(radarGrid)[c("xmin", "xmax")],
+    ylim = st_bbox(radarGrid)[c("ymin", "ymax")]
+  )
+```
+
+
+In this map we see high density of birds along the coast this corresponds to [citizen science observations](https://trektellen.nl/maps/species/-2/0/0/20221019/20221019/1/0) along the coast. This is mainly due to an exceptional number of redwing passing by as noted by [SOVON](https://www.vogelbescherming.nl/actueel/bericht/koperwiekentrek).
+
+# Animation
+
+For the animation we first generate a vector with all timestamps we want to visualize
+```{r animationPPIlist}
+animation_interval <- interval(
+  time - lubridate::minutes(180),
+  time + lubridate::minutes(20)
+)
+animation_timestamps <- seq(
+  int_start(animation_interval),
+  int_end(animation_interval),
+  "5 mins"
+)
+```
+
+The next step is a function to translate each timestamps in a dataset for visualization. This function is based on the example before, however extra attention has been given to make the code memory efficient. For larger timeseries you might want to implement caching of intermediate results. 
+
+```{r visfun}
+vid_for_timestamp <- function(x) {
+  pvol_list <- getRad::get_pvol(
+    radarCode,
+    x,
+    param = c("DBZH", "VRADH", "RHOHV"),
+    elev_max = 89
+  ) |>
+    setNames(radarCode)
+  vp_list <- purrr::map(pvol_list, bioRad::calculate_vp, warnings = F)
+  ppi_list <- purrr::map2(
+    pvol_list,
+    vp_list,
+    bioRad::integrate_to_ppi,
+    raster = radarGrid
+  ) |>
+    purrr::map(purrr::pluck, "data") |>
+    purrr::map(stars::st_as_stars) |>
+    c(list(along = "radar")) |>
+    do.call(what = "c") |>
+    mutate(VID = dplyr::if_else(R > 10, NA, VID)) |>
+    select(R, VID)
+  minR <- st_apply(
+    select(ppi_list, R),
+    1:2,
+    function(x) {
+      if (all(is.na(x))) {
+        x
+      } else {
+        x == min(x, na.rm = T)
+      }
+    },
+    keep = T
+  ) |>
+    setNames("minR") |>
+    aperm(c(2, 3, 1))
+  ppi_minR <- c(ppi_list, minR) |>
+    mutate(VID = VID * minR) |>
+    st_apply(c("x", "y"), sum, na.rm = T) |>
+    select(VID) |>
+    sf::st_as_sf(as_points = T) |>
+    tibble::add_column(time = x) |>
+    filter(!is.na(VID)) |>
+    mutate(
+      X = sf::st_coordinates(geometry)[, "X"],
+      Y = sf::st_coordinates(geometry)[, "Y"]
+    ) |>
+    sf::st_drop_geometry()
+}
+```
+Using insistently we can ensure that temporal failures are tried, this makes the code slightly more robust for a failed download attempt
+```{r ins}
+vid_for_timestamp_insistently <- purrr::insistently(
+  vid_for_timestamp,
+  quiet = FALSE,
+  rate = purrr::rate_backoff(
+    max_times = 5,
+    pause_base = 5,
+    jitter = F,
+    pause_cap = 60
+  )
+)
+```
+
+```{r anim_list}
+ppi_df <- purrr::map(
+  .progress = FALSE,
+  animation_timestamps,
+  vid_for_timestamp_insistently
+)
+```
+```{r plotting_prep}
+ppi_df_summarized <- ppi_df |>
+  bind_rows()
+```
+
+
+ Using the `gganimate` package we can create an animation from all timestamps. 
+```{r plottingR}
+require(gganimate)
+animation <- ggplot() +
+  geom_tile(
+    data = ppi_df_summarized,
+    aes(fill = VID, x = X, y = Y),
+    color = NA
+  ) +
+  scale_fill_viridis_c(
+    name = "Bird density [1/km^2]",
+    option = "turbo",
+    limit = c(2, 800),
+    transform = "log10"
+  ) +
+  geom_sf(data = map, fill = NA, color = "black", size = .3) +
+  geom_sf(data = radarLocation) +
+  coord_sf(
+    expand = F,
+    xlim = st_bbox(radarGrid)[c("xmin", "xmax")],
+    ylim = st_bbox(radarGrid)[c("ymin", "ymax")]
+  ) +
+  labs(
+    title = "Bird density over the Netherlands",
+    subtitle = "Time: {current_frame}",
+  ) +
+  transition_manual(time)
+```
+
+This animation can then be controlled and printed in different ways.  For example, as a `gif` or video file. 
+
+```{r printAnimation1R}
+#| fig-height: 9
+animate(
+  animation,
+  nframes = length(animation_timestamps),
+  fps = 5
+)
+```
+
+In this animation we again see the high density of migratory birds along the coast. This corresponds to observations by a radar (MR1) located there. At the same time a similar radar located in Amsterdam did not observe a migration peek, in the animation we see that a high density of birds migrates south of Amterdam and over the IJselmeer but not throught Amsterdam.
+
+
+
diff --git a/case_study_3/case_study_3_reproducible.md b/case_study_3/case_study_3_reproducible.md
new file mode 100644
index 0000000..415e25b
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible.md
@@ -0,0 +1,326 @@
+Case Study weather radar
+================
+
+## Setup
+
+``` r
+require(bioRad)
+require(getRad)
+require(lubridate)
+require(dplyr)
+require(sf)
+require(ggplot2)
+require(stars)
+```
+
+In this case study we zoom in to the morning of the 19th of October
+2022. We select both radar in the Netherlands. For a quick scan of the
+data the crow visualizer
+[aloft](https://crow.aloftdata.eu#/?radar=bejab&date=2022-10-19&interval=24&timedisplay=radarLocal&vpiMode=mtr&vpColorScheme=viridis&lang=en)
+repository can be consulted.
+
+``` r
+time <- as_datetime("2022-10-19 8:00:00")
+radarCode <- c("nldhl", "nlhrw")
+radarLocation <- get_weather_radars() |>
+  dplyr::filter(radar %in% radarCode) |>
+  st_transform(crs = 28992)
+radarGrid <- raster::raster(st_buffer(radarLocation, units::set_units(105, "km")),
+  res = 500
+)
+```
+
+On this night high densities of numbers were counted along the dutch
+coast, as seen in the following video on youtube.
+
+[![Watch the
+video](https://img.youtube.com/vi/8x4Dos7GbvI/default.jpg)](https://youtu.be/8x4Dos7GbvI)
+
+For visualizations we load a polygon map, in this case from the
+`cbsodataR` package. This is a slightly nicer map compared to
+`rnaturalearth` which globally provides reasonable maps.
+
+``` r
+require(cbsodataR)
+```
+
+    ## Loading required package: cbsodataR
+
+``` r
+map <- cbs_get_sf("provincie", year = 2022)
+```
+
+# Mornings activity
+
+To explore the mornings activbity we can use the vpts.
+
+``` r
+vpts <- get_vpts(radar = radarCode, as.Date(time))
+plot(vpts$nlhrw, main="Herwijnen")
+```
+
+![](case_study_3_reproducible_files/figure-gfm/vpts-1.png)<!-- -->
+
+``` r
+plot(vpts$nldhl, main="Den Helder")
+```
+
+![](case_study_3_reproducible_files/figure-gfm/vpts-2.png)<!-- -->
+
+# Static map
+
+The first step is to download radar data and calculate vertical
+distributions from the data
+
+``` r
+# Set names is to ensure the resulting list is named according to the radar code
+pvol_list <- setNames(
+  get_pvol(
+    radarCode, time,
+    elev_max = 89
+  ), radarCode
+)
+vp_list <- purrr::map(pvol_list, calculate_vp)
+```
+
+    ## Running vol2birdSetUp
+
+    ## Warning: radial velocities will be dealiased...
+
+    ## Running vol2birdSetUp
+
+    ## Warning: radial velocities will be dealiased...
+
+The next step is to create for each radar separately a vertically
+integrated estimate of density. Here we convert all cases where strong
+corrections have been applied to `NA`
+
+``` r
+ppi_strs <- purrr::map2(pvol_list, vp_list, integrate_to_ppi, raster = radarGrid) |>
+  purrr::map(purrr::pluck, "data") |>
+  purrr::map(stars::st_as_stars) |>
+  c(list(along = "radar")) |>
+  do.call(what = "c") |>
+  mutate(VID = dplyr::if_else(R > 10, NA, VID))
+```
+
+    ## Loading required namespace: raster
+
+To combine both radars we create a dataset where we identify which radar
+on which location as been minimally corrected.
+
+``` r
+minR <- st_apply(select(ppi_strs, R), 1:2, function(x) {
+  if (all(is.na(x))) {
+    x
+  } else {
+    x == min(x, na.rm = T)
+  }
+}, keep = T) |>
+  setNames("minR") |>
+  aperm(c(2, 3, 1))
+```
+
+By multiplying the VID with the minimal R value (0/1)and summing the
+values for each location we can get the best VID estimate for each
+location
+
+``` r
+ppi_minR <- c(ppi_strs, minR) |>
+  mutate(VID = VID * minR) |>
+  st_apply(c("x", "y"), sum, na.rm = T)
+```
+
+``` r
+ggplot() +
+  geom_stars(data = ppi_minR, aes(fill = VID), color = NA, sf = T) +
+  scale_fill_viridis_c(
+    name = "Bird density [1/km^2]",
+    option = "turbo",
+    limit = c(2, 800),
+    transform = "log10"
+  ) +
+  geom_sf(data = map, fill = NA, color = "black", size = .3) +
+  geom_sf(data=radarLocation)+
+  ggspatial::annotation_scale() +
+  coord_sf(
+    expand = F,
+    xlim = st_bbox(radarGrid)[c("xmin", "xmax")],
+    ylim = st_bbox(radarGrid)[c("ymin", "ymax")]
+  )
+```
+
+    ## Loading required namespace: raster
+
+    ## Warning in scale_fill_viridis_c(name = "Bird density [1/km^2]", option =
+    ## "turbo", : log-10 transformation introduced infinite values.
+
+![](case_study_3_reproducible_files/figure-gfm/exampleRplot-1.png)<!-- -->
+
+In this map we see high density of birds along the coast this
+corresponds to [citizen science
+observations](https://trektellen.nl/maps/species/-2/0/0/20221019/20221019/1/0)
+along the coast. This is mainly due to an exceptional number of redwing
+passing by as noted by
+[SOVON](https://www.vogelbescherming.nl/actueel/bericht/koperwiekentrek).
+
+# Animation
+
+For the animation we first generate a vector with all timestamps we want
+to visualize
+
+``` r
+animation_interval <- interval(
+  time - lubridate::minutes(180),
+  time + lubridate::minutes(20)
+)
+animation_timestamps <- seq(
+  int_start(animation_interval),
+  int_end(animation_interval),
+  "5 mins"
+)
+```
+
+The next step is a function to translate each timestamps in a dataset
+for visualization. This function is based on the example before, however
+extra attention has been given to make the code memory efficient. For
+larger timeseries you might want to implement caching of intermediate
+results.
+
+``` r
+vid_for_timestamp <- function(x) {
+  pvol_list <- getRad::get_pvol(
+    radarCode,
+    x,
+    param = c("DBZH", "VRADH", "RHOHV"),
+    elev_max = 89
+  ) |>
+    setNames(radarCode)
+  vp_list <- purrr::map(pvol_list, bioRad::calculate_vp, warnings = F)
+  ppi_list <- purrr::map2(
+    pvol_list,
+    vp_list,
+    bioRad::integrate_to_ppi,
+    raster = radarGrid
+  ) |>
+    purrr::map(purrr::pluck, "data") |>
+    purrr::map(stars::st_as_stars) |>
+    c(list(along = "radar")) |>
+    do.call(what = "c") |>
+    mutate(VID = dplyr::if_else(R > 10, NA, VID)) |>
+    select(R, VID)
+  minR <- st_apply(
+    select(ppi_list, R),
+    1:2,
+    function(x) {
+      if (all(is.na(x))) {
+        x
+      } else {
+        x == min(x, na.rm = T)
+      }
+    },
+    keep = T
+  ) |>
+    setNames("minR") |>
+    aperm(c(2, 3, 1))
+  ppi_minR <- c(ppi_list, minR) |>
+    mutate(VID = VID * minR) |>
+    st_apply(c("x", "y"), sum, na.rm = T) |>
+    select(VID) |>
+    sf::st_as_sf(as_points = T) |>
+    tibble::add_column(time = x) |>
+    filter(!is.na(VID)) |>
+    mutate(
+      X = sf::st_coordinates(geometry)[, "X"],
+      Y = sf::st_coordinates(geometry)[, "Y"]
+    ) |>
+    sf::st_drop_geometry()
+}
+```
+
+Using insistently we can ensure that temporal failures are tried, this
+makes the code slightly more robust for a failed download attempt
+
+``` r
+vid_for_timestamp_insistently <- purrr::insistently(
+  vid_for_timestamp,
+  quiet = FALSE,
+  rate = purrr::rate_backoff(
+    max_times = 5,
+    pause_base = 5,
+    jitter = F,
+    pause_cap = 60
+  )
+)
+```
+
+``` r
+ppi_df <- purrr::map(
+  .progress = FALSE,
+  animation_timestamps,
+  vid_for_timestamp_insistently
+)
+```
+
+    ## Loading required namespace: raster
+
+``` r
+ppi_df_summarized <- ppi_df |>
+  bind_rows()
+```
+
+Using the `gganimate` package we can create an animation from all
+timestamps.
+
+``` r
+require(gganimate)
+animation <- ggplot() +
+  geom_tile(
+    data = ppi_df_summarized,
+    aes(fill = VID, x = X, y = Y),
+    color = NA
+  ) +
+  scale_fill_viridis_c(
+    name = "Bird density [1/km^2]",
+    option = "turbo",
+    limit = c(2, 800),
+    transform = "log10"
+  ) +
+  geom_sf(data = map, fill = NA, color = "black", size = .3) +
+  geom_sf(data = radarLocation) +
+  coord_sf(
+    expand = F,
+    xlim = st_bbox(radarGrid)[c("xmin", "xmax")],
+    ylim = st_bbox(radarGrid)[c("ymin", "ymax")]
+  ) +
+  labs(
+    title = "Bird density over the Netherlands",
+    subtitle = "Time: {current_frame}",
+  ) +
+  transition_manual(time)
+```
+
+    ## Loading required namespace: raster
+
+This animation can then be controlled and printed in different ways. For
+example, as a `gif` or video file.
+
+``` r
+animate(
+  animation,
+  nframes = length(animation_timestamps),
+  fps = 5
+)
+```
+
+    ## Warning in scale_fill_viridis_c(name = "Bird density [1/km^2]", option =
+    ## "turbo", : log-10 transformation introduced infinite values.
+
+![](case_study_3_reproducible_files/figure-gfm/printAnimation1R-1.gif)<!-- -->
+
+In this animation we again see the high density of migratory birds along
+the coast. This corresponds to observations by a radar (MR1) located
+there. At the same time a similar radar located in Amsterdam did not
+observe a migration peek, in the animation we see that a high density of
+birds migrates south of Amterdam and over the IJselmeer but not throught
+Amsterdam.
diff --git a/case_study_3/case_study_3_reproducible_files/figure-gfm/compressed.gif b/case_study_3/case_study_3_reproducible_files/figure-gfm/compressed.gif
new file mode 100644
index 0000000..2f5c750
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible_files/figure-gfm/compressed.gif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:06d444a2a0aa3fcb164dd9dc16a2f225da4ca6fe3a47cac98be2da21444d11d2
+size 15184199
diff --git a/case_study_3/case_study_3_reproducible_files/figure-gfm/exampleRplot-1.png b/case_study_3/case_study_3_reproducible_files/figure-gfm/exampleRplot-1.png
new file mode 100644
index 0000000..7330ec3
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible_files/figure-gfm/exampleRplot-1.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:62d79e15961cb89c7a561ac58cb57817f5998e9d76642f392d32d30f031c2651
+size 762488
diff --git a/case_study_3/case_study_3_reproducible_files/figure-gfm/printAnimation1R-1.gif b/case_study_3/case_study_3_reproducible_files/figure-gfm/printAnimation1R-1.gif
new file mode 100644
index 0000000..e0c9d43
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible_files/figure-gfm/printAnimation1R-1.gif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2506a7a4d6fec552bd38902718be77fe49e09999159c6a8e09ad5b708d945d60
+size 16802155
diff --git a/case_study_3/case_study_3_reproducible_files/figure-gfm/vpts-1.png b/case_study_3/case_study_3_reproducible_files/figure-gfm/vpts-1.png
new file mode 100644
index 0000000..ab8eaa9
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible_files/figure-gfm/vpts-1.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6c33483f408bbe74cd6fbb7ac223ae9d9b98e8a809660a6ef3df2a68ea5f7edc
+size 33033
diff --git a/case_study_3/case_study_3_reproducible_files/figure-gfm/vpts-2.png b/case_study_3/case_study_3_reproducible_files/figure-gfm/vpts-2.png
new file mode 100644
index 0000000..fc02745
--- /dev/null
+++ b/case_study_3/case_study_3_reproducible_files/figure-gfm/vpts-2.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e54f1bcdf04e14f5ffedf9746ea15744621445af879a65bd818dbbd9540d9b4b
+size 39611