Add energy figure

vignettes/wind-graph.Rmd

@@ -22,10 +22,12 @@ library(ncdf4)
 library(leaflet)
 library(leaflet.extras)
 library(igraph)
+library(ggplot2)
+library(gridExtra)
 grl <- readRDS(system.file("extdata", "18LX_grl.rda", package = "GeoPressureR"))
 ```
 
-In this vignette, we will enhance the [basic graph](./articles/basic-graph) by adding the wind speed. This allows to defines the movement model on the airspeed (instead of groundspeed). In addition, we can extract windsupport information on the journey of the bird. 
+In this vignette, we will enhance the [basic graph](./articles/basic-graph) by adding the wind speed. This allows to defines the movement model on the airspeed (instead of groundspeed). In addition, we can extract windsupport information on the journey of the bird.
 
 ```{r}
 pam_data <- pam_read(
@@ -39,11 +41,10 @@ static_prob <- readRDS(system.file("extdata", "18LX_static_prob.rda", package =
 
 ## Download wind data
 
-Wind data is available at high resolution (1hr, 0.25°, 37 pressure level) on [ERA5 hourly data on pressure levels](https://doi.org/10.24381/cds.bd0915c6). And this data is easily accessible through the [`ecmfr`](https://bluegreen-labs.github.io/ecmwfr) package. 
-As the flight are of short duration, we suggest to download a file for each flight.
+Wind data is available at high resolution (1hr, 0.25°, 37 pressure level) on [ERA5 hourly data on pressure levels](https://doi.org/10.24381/cds.bd0915c6). And this data is easily accessible through the [`ecmfr`](https://bluegreen-labs.github.io/ecmwfr) package. As the flight are of short duration, we suggest to download a file for each flight.
 
 The first step is to setup-up your CDS access. You will need to create an account on [https://cds.climate.copernicus.eu/](https://cds.climate.copernicus.eu/user/register) to generate your API key and uid number. You can stored them in your `.Rprofile` with the commented code.
- 
+
 ```{r, eval=F}
 #Sys.setenv( cds.key="Insert_your_CDS_API_KEY_here")
 #Sys.setenv( cds.user="Insert_your_CDS_UID_here")
@@ -62,7 +63,7 @@ area = extent(static_prob[[1]])
 area = c(area@ymax, area@xmin, area@ymin, area@xmax)
 ```
 
-It will be faster to send all the requests (one per flight) with `wf_request()` using `transfer = F` without waiting to get the data from one before requesting the next. These query can be very quick or sometime take a couple of hours, usually not more. You can monitor easily at [https://cds.climate.copernicus.eu/cdsapp#!/yourrequests](https://cds.climate.copernicus.eu/cdsapp#!/yourrequests) the status of your requests, delete them or download them manually if needed. 
+It will be faster to send all the requests (one per flight) with `wf_request()` using `transfer = F` without waiting to get the data from one before requesting the next. These query can be very quick or sometime take a couple of hours, usually not more. You can monitor easily at <https://cds.climate.copernicus.eu/cdsapp#!/yourrequests> the status of your requests, delete them or download them manually if needed.
 
 ```{r, eval=F}
 req <- list()
@@ -101,7 +102,7 @@ Define the folder where to download the data. `wf_transfer()` will download the
 dir.save <- '~'
 ```
 
-The following code will return the status of the request if then have not yet been processed, or return the data otherwise. I suggest to check on [https://cds.climate.copernicus.eu/cdsapp#!/yourrequests] if the request have been completed before running this code. It can take a couple of minutes to a few hours.    
+The following code will return the status of the request if then have not yet been processed, or return the data otherwise. I suggest to check on [<https://cds.climate.copernicus.eu/cdsapp#!/yourrequests>] if the request have been completed before running this code. It can take a couple of minutes to a few hours.
 
 ```{r, eval=F}
 for (i_s in seq_len(nrow(pam_data$sta)-1)){
@@ -118,10 +119,9 @@ We first create the graph identically to in [basic graph](./articles/basic-graph
 grl <- graph_create(static_prob, thr_prob_percentile = .99, thr_gs = 150)  
 ```
 
-
 ## Add wind to graph
 
-We can compute the windspeed experienced by the bird if he had flew each possible transition (i.e. edge in the graph). Based on this windspeed and groundspeed, we also compute the airspeed. All of these are stored as complex value with the real part representing the E-W component and the imaginary part corresponding to the N-S. 
+We can compute the windspeed experienced by the bird if he had flew each possible transition (i.e. edge in the graph). Based on this windspeed and groundspeed, we also compute the airspeed. All of these are stored as complex value with the real part representing the E-W component and the imaginary part corresponding to the N-S.
 
 ```{r, eval=F}
 filename = paste0(dir.save,"/","18IC_")
@@ -130,20 +130,20 @@ grl <- graph_add_wind(grl, pressure=pam_data$pressure, filename, thr_as = 100)
 
 ## Compute the edges probability
 
-Now that the have computed the airspeed required for performing the transition of each edge, we can improve the computation of the probability by modeling the probability of airspeed rather than groundspeed. 
+Now that the have computed the airspeed required for performing the transition of each edge, we can improve the computation of the probability by modeling the probability of airspeed rather than groundspeed.
 
 We first search the morphological information of the Great Reed Warbler using the [AVONET database](https://doi.org/10.6084/m9.figshare.16586228.v5). You can also overwrite any of these value if you know them. See `flight_bird()` for more details
+
 ```{r}
 bird <- flight_bird("Acrocephalus arundinaceus")
 speed <- seq(0,80)
 prob <- flight_prob(speed, method = "power", bird = bird, low_speed_fix = 10, 
                     fun_power = function(power) { (1 / power)^3 })
 plot(speed, prob, type="l", xlab="Airspeed [km/h]", ylab="Probability")
-grl$p <- grl$ps * flight_prob(grl$gs, method = "power", bird = bird, low_speed_fix = 10)
+grl$p <- grl$ps * flight_prob(grl$as, method = "power", bird = bird, low_speed_fix = 10)
 ```
 
-
-## Output 1: Shortest path
+## Output 1: Shortest path with wind
 
 In graph theory, [the shortest path](https://en.wikipedia.org/wiki/Shortest_path_problem) correspond to the set of nodes whose sum of the edges weights are as small as possible. By weighting the edges with the minus of the log of the probability, this corresponds to finding the most likely trajectory of our bird. We use the [igraph package](https://igraph.org/r/) to compute the shortestpath
 
@@ -153,13 +153,22 @@ g <- graph_from_data_frame(data.frame(
   to = grl$t,
   weight = -log(grl$p)
 ))
-sp <- shortest_paths(g, from = paste(grl$equipement), to = paste(grl$retrival))
+# In case there are no retrival site, we select the position with the hight probability according to the marginal
+# retrival <- which.max(as.matrix(static_prob_marginal[[length(static_prob_marginal)]])) + grl$sz[1] * grl$sz[2] * (grl$sz[3] - 1)
+# stopifnot(retrival %in% grl$retrival)
+retrival <- grl$retrival
+sp <- shortest_paths(g, from = paste(grl$equipement), to = paste(retrival))
+
 # Convert igraph representation to lat-lon
 grl$shortest_path <- graph_path2lonlat(as.numeric(sp$vpath[[1]]$name), grl)
 ```
 
 We can visualize the shortest path with the windpseed direction (arrow) and magnitude (color) experienced during this particular flight
 
+<details>
+
+<summary>See `fun_marker_color()` and `fun_NSEW()`</summary>
+
 ```{r}
 fun_marker_color <- function(norm){
     if (norm < 20){
@@ -203,6 +212,8 @@ fun_NSEW <- function(angle){
 }
 ```
 
+</details>
+
 ```{r,}
 sta_duration <- unlist(lapply(static_prob,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))
 
@@ -236,28 +247,7 @@ for (i_s in seq_len(grl$sz[3]-1)){
 m
 ```
 
-
-
-```{r, eval = F}
-path <- as.data.frame(grl$shortest_path)
-path$sta_id <- seq_len(nrow(path))
-static_timeserie <- geopressure_ts_path(path, pam_data$pressure)
-geopressureviz = list(
-  pam_data=pam_data,
-  static_prob=static_prob,
-  pressure_prob=pressure_prob,
-  light_prob=light_prob,
-  pressure_timeserie = static_timeserie
-)
-```
-
-```{r, eval=F}
-appDir <- system.file("geopressureviz", package = "GeoPressureR")
-shiny::runApp(appDir, launch.browser = getOption("browser"))
-# shiny::runApp('./inst/geopressureviz', launch.browser = getOption("browser"))
-```
-
-## Output 2: Proability map of stationary period
+## Output 2: Proability map of stationary period in GeoPressureViz
 
 Estimating the position of the bird for each stationary period is generally the most sought-after output of tracking studies. Using the graph built, we can compute this exactly (i.e., without iterative approach such as MCMC). This problem is the same as computing the marginal distribution of a Markov process which can be solved mathematically.
 
@@ -278,33 +268,99 @@ for (i_r in seq_len(length(grl_marginal))) {
     addCircles(lng = grl$shortest_path$lon[i_s], lat = grl$shortest_path$lat[i_s], opacity = 1, color = "#000", weight = 10, group = info_str)
 }
 l %>%
+  addPolylines(lng = grl$shortest_path$lon, lat = grl$shortest_path$lat, opacity = .5, color = "#808080", weight = 0.5) %>%
+  addCircles(lng = grl$shortest_path$lon, lat = grl$shortest_path$lat, opacity = .5, color = "#000", weight = sta_duration^(0.3)*10) %>% 
   addLayersControl(
     overlayGroups = li_s,
     options = layersControlOptions(collapsed = FALSE)
   ) %>%
   hideGroup(tail(li_s, length(li_s) - 1))
 ```
 
-## Output 3: Simulate path
+We can also visualize the shortest path and marginal map in GeoPressureViz. There are few preparation step before.
+
+First, we need to query the pressure timeserie on the shortest path location. 
+
+```{r, eval = F}
+shortest_path <- as.data.frame(grl$shortest_path)
+shortest_path_timeserie <- geopressure_ts_path(shortest_path, pam_data$pressure)
+```
+
+Then, we need to load a lot of the previous data and aligned them (same stationary period). We finally write the file in your home directory.
+
+```{r, eval = F}
+grl <- readRDS(system.file("extdata", "18LX_grl.rda", package = "GeoPressureR"))
+grl_marginal <- readRDS(system.file("extdata", "18LX_grl_marginal.rda", package = "GeoPressureR"))
+shortest_path_timeserie <- readRDS(system.file("extdata", "18LX_shortest_path_timeserie.rda", package = "GeoPressureR"))
+static_prob <- readRDS(system.file("extdata", "18LX_static_prob.rda", package = "GeoPressureR"))
+light_prob <- readRDS(system.file("extdata", "18LX_light_prob.rda", package = "GeoPressureR")) 
+pressure_prob <- readRDS(system.file("extdata", "18LX_pressure_prob.rda", package = "GeoPressureR")) 
+
+sta_marginal <- unlist(lapply(grl_marginal, function(x) raster::metadata(x)$sta_id))
+sta_pres <- unlist(lapply(pressure_prob, function(x) raster::metadata(x)$sta_id))
+sta_light <- unlist(lapply(light_prob, function(x) raster::metadata(x)$sta_id))
+pressure_prob <- pressure_prob[sta_pres %in% sta_marginal]
+light_prob <- light_prob[sta_light %in% sta_marginal]
+
+geopressureviz = list(
+  pam_data=pam_data,
+  static_prob=static_prob,
+  static_prob_marginal = grl_marginal,
+  pressure_prob=pressure_prob,
+  light_prob=light_prob,
+  pressure_timeserie = shortest_path_timeserie
+)
+save(geopressureviz, file = "~/geopressureviz.RData")
+```
+
+And start the app!
+
+```{r, eval=F}
+appDir <- system.file("geopressureviz", package = "GeoPressureR")
+shiny::runApp(appDir, launch.browser = getOption("browser"))
+# shiny::runApp('./inst/geopressureviz', launch.browser = getOption("browser"))
+```
+
+## Output 3: Simulate path and compute energy
+
+We can compute a few simulation paths. Because they are independent and without error, so you don't need many, although the computation to request more is almost the same as a few. 
 
 ```{r, results='hide'}
 nj=10
 path <- graph_simulation(grl, nj)
 ```
+
+We can compute the energy [Joules] expense of each flight for the bird
+
 ```{r}
-col <- rep(RColorBrewer::brewer.pal(9, "Set1"), times = ceiling(grl$sz[3]/9))
-m <- leaflet(width = "100%") %>% addProviderTiles(providers$Stamen.TerrainBackground) %>%  addFullscreenControl()
-for (i in seq_len(nj)) {
-  m <- m %>%
-    addPolylines(lng = path$lon[i,], lat = path$lat[i,], opacity = 0.7, weight = 1, color = "#808080") 
-}
-for (i in seq_len(grl$sz[3])){
-   m <- m %>%
-    addCircles(lng = path$lon[,i], lat = path$lat[,i], opacity = .4, weight = 10, color = col[i])
-}
-m <- m %>% leaflet::addLegend("bottomright", colors = col[1:grl$sz[3]], labels = seq_len(grl$sz[3]),
-    title = "Stationay period",
-    opacity = 1
-  ) 
-  m
+edge <- graph_path2edge(path$id, grl)
+
+# Convert airpseed from km/hr to m/s and use the bird constructure earlier to compute the mechanical power [W=J/s]
+p_mech <- flight_power(abs(grl$as[edge]) * 1000 / 60 / 60, bird = bird)
+
+# Convert the power into energy [J]
+energy <- p_mech * rep(head(grl$flight_duration,-1)*60*60,nj)
+dim(energy) <- dim(edge)
+```
+
+And compare the histogram of ground, wind, airspeed, flight duration and energy for each of the simulation. 
+
+```{r, fig.height=5}
+energy_df <- data.frame(
+  energy = as.vector(energy),
+  as = abs(grl$as[edge]),
+  gs = abs(grl$gs[edge]),
+  ws = abs(grl$ws[edge]),
+  sta_id_s = rep(head(grl$sta_id,-1), nj),
+  sta_id_t = rep(tail(grl$sta_id,-1), nj),
+  flight_duration = rep(head(grl$flight_duration,-1), nj)
+)
+energy_df$name = paste(energy_df$sta_id_s,energy_df$sta_id_t, sep=" ->")
+
+plot1 <- ggplot(energy_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot()
+plot2 <- ggplot(energy_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() 
+plot3 <- ggplot(energy_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() 
+plot4 <- ggplot(energy_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point()
+plot5 <- ggplot(energy_df, aes(reorder(name, sta_id_s), energy)) + geom_boxplot() 
+grid.arrange(plot1, plot2, plot3, plot4, plot5, nrow=5)
 ```