Pressure Trajectories
library(GeoPressureR)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(plotly)
library(GeoLocTools)
setupGeolocation()
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_wind_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_grl.Rdata"))
col <- rep(RColorBrewer::brewer.pal(8, "Dark2"), times = ceiling(max(pam$sta$sta_id) / 8))
The large circles indicate the shortest path (most likely trajectory) estimated by the graph approach, with surrounding smaller points and lines representing other simulated paths. The size of each point is proportional to the duration of stay.
sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
}))
pal <- colorFactor(col, as.factor(seq_len(length(col))))
m <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl()
for (i in seq_len(nrow(path_sim$lon))) {
m <- m %>%
addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.5, weight = 1, color = "#808080") %>%
addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = .7, weight = 1, color = pal(factor(shortest_path$sta_id, levels = pam$sta$sta_id)))
}
m %>%
addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = pal(factor(shortest_path$sta_id, levels = pam$sta$sta_id)), weight = sta_duration^(0.3) * 10, popup = paste0("sta_id=", shortest_path$sta_id))
The marginal probability map estimates the overall probability of position at each stationary period. It is the most useful quantification of the uncertainty of the position of the bird. Select each stationary period (labeled by the first number, followed by the dates of duration) to view the progression of the full annual cycle.
li_s <- list()
l <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl()
for (i_r in seq_len(length(static_prob_marginal))) {
i_s <- metadata(static_prob_marginal[[i_r]])$sta_id
info <- metadata(static_prob_marginal[[i_r]])$temporal_extent
info_str <- paste0(i_s, " | ", info[1], "->", info[2])
li_s <- append(li_s, info_str)
l <- l %>%
addRasterImage(static_prob_marginal[[i_r]], colors = "OrRd", opacity = 0.8, group = info_str) #method = "ngb") #%>%
#addCircles(lng = shortest_path$lon[i_s], lat = shortest_path$lat[i_s], opacity = 1, color = "#000", weight = 10, group = info_str)
}
l %>%
addLayersControl(
overlayGroups = li_s,
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(tail(li_s, length(li_s) - 1))
Altitudes are computed based on pressure measurements, corrected based on the assumed location of the shortest path. This correction accounts for the natural variation of pressure as estimated by ERA-5.
Each stationary period is indicated by a separate color, while each migratory period is colored black. Click and drag a viewing area to zoom into different events for further details.
p <- ggplot() +
# geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
geom_line(data = do.call("rbind", shortest_path_timeserie), aes(x = date, y = altitude)) +
geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = altitude, col = factor(sta_id))) +
theme_bw() +
scale_colour_manual(values = col) +
scale_y_continuous(name = "Altitude (m)")
ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)
Wind data (ERA5) displayed during migratory periods along the shortest path simulation.
fun_marker_color <- function(norm){
if (norm < 20){
"darkpurple"
} else if (norm < 35){
"darkblue"
} else if (norm < 50){
"lightblue"
} else if (norm < 60){
"lightgreen"
} else if (norm < 80){
"yellow"
} else if (norm < 100){
"lightred"
} else {
"darkred"
}
}
fun_NSEW <- function(angle){
angle <- angle %% (pi* 2)
angle <- angle*180/pi
if (angle < 45/2){
"E"
} else if (angle < 45*3/2){
"NE"
} else if (angle < 45*5/2){
"N"
} else if (angle < 45*7/2){
"NW"
} else if (angle < 45*9/2){
"W"
} else if (angle < 45*11/2){
"SW"
} else if (angle < 45*13/2){
"S"
}else if (angle < 45*15/2){
"SE"
} else {
"E"
}
}
sta_duration <- unlist(lapply(static_prob_marginal,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))
m <-leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>% addFullscreenControl() %>%
addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3)*10)
for (i_s in seq_len(grl$sz[3]-1)){
if (grl$flight_duration[i_s]>5){
edge <- which(grl$s == shortest_path$id[i_s] & grl$t == shortest_path$id[i_s+1])
label = paste0( i_s,': ', grl$flight[[i_s]]$start, " - ", grl$flight[[i_s]]$end, "<br>",
"F. dur.: ", round(grl$flight_duration[i_s]), ' h <br>',
"GS: ", round(abs(grl$gs[edge])), ' km/h, ',fun_NSEW(Arg(grl$gs[edge])),'<br>',
"WS: ", round(abs(grl$ws[edge])), ' km/h, ',fun_NSEW(Arg(grl$ws[edge])),'<br>',
"AS: ", round(abs(grl$as[edge])), ' km/h, ',fun_NSEW(Arg(grl$as[edge])),'<br>')
iconArrow <- makeAwesomeIcon(icon = "arrow-up",
library = "fa",
iconColor = "#FFF",
iconRotate = (90 - Arg(grl$ws[edge])/pi*180) %% 360,
squareMarker = TRUE,
markerColor = fun_marker_color(abs(grl$ws[edge])))
m <- m %>% addAwesomeMarkers(lng = (shortest_path$lon[i_s] + shortest_path$lon[i_s+1])/2,
lat = (shortest_path$lat[i_s] + shortest_path$lat[i_s+1])/2,
icon = iconArrow, popup = label)
}
}
m
Comparison of each migratory flight duration (hours), ground speed (gs, km/h), wind speed (ws, km/h), air speed (as, km/h), maximum altitude (m), and mean altitude (m). Each flight is labeled as the starting stationary period to the next. Hover over each plot for details
edge <- t(graph_path2edge(path_sim$id, grl))
nj <- ncol(edge)
nsta <- ncol(path_sim$lon)
speed_df <- data.frame(
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),
dist = geosphere::distGeo(
cbind(as.vector(t(path_sim$lon[,1:nsta-1])), as.vector(t(path_sim$lat[,1:nsta-1]))),
cbind(as.vector(t(path_sim$lon[,2:nsta])), as.vector(t(path_sim$lat[,2:nsta])))
) / 1000
) %>% mutate(
name = paste(sta_id_s,sta_id_t, sep="-")
)
alt_df = do.call("rbind", shortest_path_timeserie) %>%
arrange(date) %>%
mutate(
sta_id_s = cummax(sta_id),
sta_id_t = sta_id_s+1
) %>%
filter(sta_id == 0 & sta_id_s > 0 ) %>%
group_by(sta_id_s, sta_id_t) %>%
mutate(diff_altitude = altitude-lag(altitude, default = altitude[1])) %>%
summarise(
alt_min = min(altitude),
alt_max = max(altitude),
alt_mean = mean(altitude),
alt_med = median(altitude),
alt_sumdabsdiff = sum(abs(diff_altitude)),
alt_sumposdiff = sum(ifelse(diff_altitude>0,diff_altitude,0)),
)
com_df <- speed_df %>%
left_join(alt_df)
plot1 <- ggplot(com_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot2 <- ggplot(com_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot3 <- ggplot(com_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot4 <- ggplot(com_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point() + theme_bw() +scale_x_discrete(name = "")
plot5 <- ggplot(com_df, aes(reorder(name, sta_id_s), alt_max)) + geom_point() + theme_bw() +scale_x_discrete(name = "")
plot6 <- ggplot(com_df, aes(reorder(name, sta_id_s), alt_mean)) + geom_point() + theme_bw() +scale_x_discrete(name = "")
subplot(ggplotly(plot4), ggplotly(plot1), ggplotly(plot2), ggplotly(plot3), ggplotly(plot5), ggplotly(plot6), nrows=6, titleY=T)
Histogram of error between recorded pressure and ERA5 for each stationary period
pressure_ts_bind <- do.call("rbind", shortest_path_timeserie) %>%
filter(!is.na(sta_id))
pam$pressure %>%
left_join(pressure_ts_bind %>% dplyr::select(c("date", "pressure0")), by = "date") %>%
mutate(diff = ifelse(is.na(pressure0), 0, obs - pressure0)) %>%
filter(sta_id > 0 & !isoutlier) %>%
group_by(sta_id) %>%
mutate(sta_id = paste0(sta_id, " (SD=", round(sd(diff), 2), " ; N=", n(), ")")) %>%
ggplot(aes(x = diff)) +
geom_histogram(aes(y = (..count..) / tapply(..count.., ..PANEL.., sum)[..PANEL..]), binwidth = .2) +
facet_wrap(~sta_id) +
scale_x_continuous(name = "Pressure Geolocator - best match ERA5 (hPa)") +
scale_y_continuous(name = "Normalized histogram")