diff --git a/Resources/GAMs.pptx b/Resources/GAMs.pptx new file mode 100644 index 00000000..a5244363 Binary files /dev/null and b/Resources/GAMs.pptx differ diff --git a/Resources/YAPS.pptx b/Resources/YAPS.pptx new file mode 100644 index 00000000..e3f2570d Binary files /dev/null and b/Resources/YAPS.pptx differ diff --git a/_episodes/12-basic-animation.md b/_episodes/12-basic-animation.md index 6a48ea79..43385bed 100644 --- a/_episodes/12-basic-animation.md +++ b/_episodes/12-basic-animation.md @@ -7,7 +7,7 @@ questions: - "How do I animate my animal movements?" --- -Static plots are excellent tools and are appropriate a lot of the time, but there are instances where something extra is needed to demonstrate your interesting fish movements: This is where plotting animated tracks can be a useful tool. In this lesson we will explore how to take data from your [OTN-style detection extract documents](https://members.oceantrack.org/data/otn-detection-extract-documentation-matched-to-animals) and animate the journey of one fish between stations. +Static plots are excellent tools and are appropriate a lot of the time, but there are instances where something extra is needed to demonstrate your interesting fish movements. This is where plotting animated tracks can be a useful tool. In this lesson we will explore how to take data from your [OTN-style detection extract documents](https://members.oceantrack.org/data/otn-detection-extract-documentation-matched-to-animals) and animate the journey of one fish between stations. ### Getting our Packages @@ -33,10 +33,12 @@ For the purposes of this lesson we will assume that any detection that did not p Finally, we use the `detection_events` function with station as the `location_col` argument to get our detection events. ~~~ +unzip('nsbs_matched_detections_2022.zip', overwrite = TRUE) + detection_events <- #create detections event variable - read_otn_detections('proj58_matched_detections_2016.csv') %>% # reading detections + read_otn_detections('nsbs_matched_detections_2022/nsbs_matched_detections_2022.csv') %>% false_detections(tf = 3600) %>% #find false detections - filter(passed_filter != FALSE) %>% + dplyr::filter(passed_filter != FALSE) %>% detection_events(location_col = 'station', time_sep=3600) ~~~ {: .language-r} @@ -49,10 +51,10 @@ plot_data <- detection_events %>% ~~~ {: .language-r} -Additionally, animating many animal tracks can be computationally intensive as well as create a potentially confusing plot, so for this lesson we will only be plotting one fish. We well subset our data by filtering where the `animal_id` is equal to `PROJ58-1218508-2015-10-13`. +Additionally, animating many animal tracks can be computationally intensive as well as create a potentially confusing plot, so for this lesson we will only be plotting one fish. We well subset our data by filtering where the `animal_id` is equal to `NSBS-1393342-2021-08-10`. ~~~ -one_fish <- plot_data[plot_data$animal_id == "PROJ58-1218508-2015-10-13",] +one_fish <- plot_data[plot_data$animal_id == "NSBS-1393342-2021-08-10",] ~~~ {: .language-r} @@ -61,7 +63,7 @@ one_fish <- plot_data[plot_data$animal_id == "PROJ58-1218508-2015-10-13",] Now that we have our data we can begin to create our plot. We will start with creating a static plot and then once happy with that, we will animate it. -The first thing we will do for our plot is download the basemap, this will be the background to our plot. To do this we will use the `get_stamenmap` function from `ggmap`. This function gets a Stamen Map based off a bounding box that we provide. To create the bounding box we will pass a vector of four values to the argument `bbox` ; those four values represent the left, bottom, right, and top boundaries of the map. +The first thing we will do for our plot is download the basemap. This will provide the background for our plot. To do this we will use the `get_stadiamap` function from `ggmap`. This function gets a Stamen Map based off a bounding box that we provide. "Stamen" is the name of the service that provides the map tiles, but it was recently bought by Stadia, so the name of the function has changed. To create the bounding box we will pass a vector of four values to the argument `bbox` ; those four values represent the left, bottom, right, and top boundaries of the map. To determine which values are needed we will use the `min` and `max` function on the `mean_longitude` and `mean_latitude` columns of our `one_fish` variable. `min(one_fish$mean_longitude)` will be our left-most bound, `min(one_fish$mean_latitude)` will be our bottom bound, `max(one_fish$mean_longitude)` will be our right-most bound, and `max(one_fish$mean_latitude)` will be our top bound. This gives most of what we need for our basemap but we can further customize our plot with `maptype` which will change what type of map we use, `crop` which will crop raw map tiles to the specified bounding box, and `zoom` which will adjust the zoom level of the map. @@ -76,55 +78,55 @@ To determine which values are needed we will use the `min` and `max` function on ~~~ basemap <- - get_stamenmap( + get_stadiamap( bbox = c(left = min(one_fish$mean_longitude), bottom = min(one_fish$mean_latitude), right = max(one_fish$mean_longitude), top = max(one_fish$mean_latitude)), - maptype = "toner-lite", + maptype = "stamen_toner_lite", crop = FALSE, - zoom = 8) + zoom = 7) ggmap(basemap) ~~~ {: .language-r} -Now that we have our basemap ready we can create our static plot. We will store our plot in a variable called `act.plot` so we can access it later on. +Now that we have our basemap ready we can create our static plot. We will store our plot in a variable called `otn.plot` so we can access it later on. -To start our plot we will call the `ggmap` function and pass it our basemap as an argument. To make our detection locations we will then call `geom_point`, supplying `one_fish` as the data argument and for the aesthetic will make the `x` argument equal to `mean_longitude` and the `y` argument will be `mean_latitude`. +To start our plot we will call the `ggmap` function and pass it our basemap as an argument. To make our detection locations we will then call `geom_point`, supplying `one_fish` as the data argument. For the aesthetic we will make the `x` argument equal to `mean_longitude` and the `y` argument will be `mean_latitude`. This will orient our map and data properly. -We will then call `geom_path` to connect those detections supplying `one_fish` as the data argument and for the aesthetic `x` will again be `mean_longitude` and `y` will be `mean_latitude`. +We will then call `geom_path` to connect those detections supplying `one_fish` as the data argument. The aesthetic `x` will again be `mean_longitude` and `y` will be `mean_latitude`. -Lastly, we will use the `labs` function to add context to our plot such as adding a `title`, a label for the `x` axis, and a label for the `y` axis. We are then ready to view our graph by calling `ggplotly` with `act.plot` as the argument! +Lastly, we will use the `labs` function to add context to our plot including a `title`, a label for the `x` axis, and a label for the `y` axis. We are then ready to view our graph by calling `ggplotly` with `otn.plot` as the argument! ~~~ -act.plot <- - ggmap(base) + - geom_point(data = one_fish2, aes(x = mean_longitude, y = mean_latitude, group = animal_id, color = animal_id), size = 2) + - geom_path(data = one_fish2, aes(x = mean_longitude, y = mean_latitude, group = animal_id, color = animal_id)) + - labs(title = "ACT animation", +otn.plot <- + ggmap(basemap) + + geom_point(data = one_fish, aes(x = mean_longitude, y = mean_latitude), size = 2) + + geom_path(data = one_fish, aes(x = mean_longitude, y = mean_latitude)) + + labs(title = "NSBS Animation", x = "Longitude", y = "Latitude", color = "Tag ID") -ggplotly(act.plot) +ggplotly(otn.plot) ~~~ {: .language-r} ### Animating our Static Plot -Once we have a static plot we are happy with we are ready for the final step of animating it! We will use the `gganimate` package for this, since it integrates nicely with `ggmap`. +Once we have a static plot we are happy with, we are ready for the final step of animating it! We will use the `gganimate` package for this, since it integrates nicely with `ggmap`. -To animate our plot we update our `act.plot` variable by using it as our base, then add a label for the dates to go along with the animation. We then call `transition_reveal`, which is a function from `gganimate` that determines how to create the transitions for the animations. There are many transitions you can use for animations with `gganimate` but `transition_reveal` will calculate intermediary values between time observations. For our plot we will pass `transition_reveal` the `first_detection` information. We will finally use the functions `shadow_mark` with the arguments of `past` equal to `TRUE` and `future` equal to `FALSE`. This makes the animation continually show the previous data (a track) but not the future data yet to be seen (allowing it to be revealed as the animation progresses). +To animate our plot we update our `otn.plot` variable by using it as our base, then add a label for the dates to go along with the animation. We then call `transition_reveal`, which is a function from `gganimate` that determines how to create the transitions for the animations. There are many transitions you can use for animations with `gganimate` but `transition_reveal` will calculate intermediary values between time observations. For our plot we will pass `transition_reveal` the `first_detection` information. We will finally use the functions `shadow_mark` with the arguments of `past` equal to `TRUE` and `future` equal to `FALSE`. This makes the animation continually show the previous data (a track) but not the future data yet to be seen (allowing it to be revealed as the animation progresses). -Finally, to see our new animation we call the `animate` function with `act.plot` as the argument. +Finally, to see our new animation we call the `animate` function with `otn.plot` as the argument. ~~~ -act.plot <- - act.plot + +otn.plot <- + otn.plot + labs(subtitle = 'Date: {format(frame_along, "%d %b %Y")}') + transition_reveal(first_detection) + shadow_mark(past = TRUE, future = FALSE) -animate(act.plot) +animate(otn.plot) ~~~ {: .language-r} diff --git a/_episodes/13-animation-with-pathroutr.md b/_episodes/13-animation-with-pathroutr.md index a9051cd2..78bd4406 100644 --- a/_episodes/13-animation-with-pathroutr.md +++ b/_episodes/13-animation-with-pathroutr.md @@ -25,44 +25,51 @@ library(pathroutr) library(ggspatial) library(sp) library(raster) +library(geodata) detection_events <- #create detections event variable - read_otn_detections('proj58_matched_detections_2016.csv') %>% # reading detections + read_otn_detections('nsbs_matched_detections_2022.csv') %>% # reading detections false_detections(tf = 3600) %>% #find false detections - filter(passed_filter != FALSE) %>% + dplyr::filter(passed_filter != FALSE) %>% detection_events(location_col = 'station', time_sep=3600) plot_data <- detection_events %>% dplyr::select(animal_id, mean_longitude,mean_latitude, first_detection) -one_fish <- plot_data[plot_data$animal_id == "PROJ58-1218518-2015-09-16",] +one_fish <- plot_data[plot_data$animal_id == "NSBS-1393342-2021-08-10",] -one_fish <- one_fish %>% filter(mean_latitude < 38.90 & mean_latitude > 38.87) %>% - slice(155:160) +~~~ +{: .language-r} +There is one small tweak we are going to make that is not immediately intuitive, and which we're only doing for the sake of this lesson. The blue sharks in our dataset have not given us many opportunities to demonstrate `pathroutr`'s awareness of coastlines. In order to give you a fuller demonstration of the package, we are going to cheat and shift the data 0.5 degrees to the west, which brings it more into contact with the Nova Scotian coastline and lets us show off `pathroutr` more completely. You do not need to do this with your real data. + +~~~ +one_fish_shifted <- one_fish %>% mutate(mean_longitude_shifted = mean_longitude-0.5) ~~~ {: .language-r} ### Getting our Shapefile -The first big difference between our basic animation lesson and this lesson is that we will need a shapefile of the study area, so `pathroutr` can determine where the landmasses are located. To do this we will use the `getData` function from the `raster` library which gets geographic data for anywhere in the world. The first argument we will pass to `getData` is the name of the dataset we want to use, for our visualization we will use `GADM`. When you use `GADM` with `getData` you need to provide a `country` argument, which is specified by a country's 3 letter ISO code. In our case we will pass `USA`. Last, we pass an argument for `level`, this is how much we would like the data to be subdivided. Since we are only interested in a portion of the country we will pass `level` the value of `1` which means that the data will be divided by states. When we run this code we will get back a `SpatialPolygonsDataFrame` with the shapefiles for the US states. +The first big difference between our basic animation lesson and this lesson is that we will need a shapefile of the study area, so `pathroutr` can determine where the landmasses are located. To do this we will use the `gadm` function from the `geodata` library which gets administrative boundaries (i.e, borders) for anywhere in the world. The first argument we will pass to `gadm` is the name of the country we wish to get, in this case, Canada. We will specify `level` as 1, meaning we want our data to be subdivided at the first level after 'country' (in this case, provinces). 0 would get us a single shapefile of the entire country; 1 will get us individual shapefiles of each province. We must also provide a path for the downloaded shapes to be stored (`./geodata` here), and optionally a resolution. `gadm` only has two possible values for resolution: 1 for 'high' and 2 for 'low'. We'll use low resolution here because as we will see, for this plot it is good enough and will reduce the size of the data objects we download. + +This is only one way to get a shapefile for our coastlines- you may find you prefer a different method. Regardless, this is the one we'll use for now. ~~~ -USA<-getData('GADM', country='USA', level=1) +CAN<-gadm('CANADA', level=1, path="./geodata", resolution=2) ~~~ {: .language-r} -Since we only need one state we will have to filter out the states we don't need. We can do this by filtering the data frame using the same filtering methods we have explored in previous lessons. +We only need one province, which we can select using the filtering methods common to R. ~~~ -shape_file <- USA[USA$NAME_1 == 'Maryland',] +shape_file <- CAN[CAN$NAME_1 == 'Nova Scotia',] ~~~ {: .language-r} This shapefile is a great start, but we need the format to be an `sf` `multipolygon`. To do that we will run the `st_as_sf` function on our shapefile. We also want to change the coordinate reference system (CRS) of the file to a projected coordinate system since we will be mapping this plot flat. To do that we will run `st_transform` and provide it the value `5070`. ~~~ -md_polygon <- st_as_sf(single_poly) %>% st_transform(5070) +ns_polygon <- st_as_sf(single_poly) %>% st_transform(5070) ~~~ {: .language-r} @@ -73,7 +80,7 @@ We will also need to make some changes to our detection data as well, in order t Using the `SpatialPoints` function we will pass our new `path` variable and `CRS("+proj=longlat +datum=WGS84 +no_defs")` for the `proj4string` argument. Just like for our shapefile we will need to turn our path into an `sf` object by using the `st_as_sf` function and change the CRS to a projected coordinate system because we will be mapping it flat. ~~~ -path <- one_fish %>% dplyr::select(mean_longitude,mean_latitude) +path <- one_fish_shifted %>% dplyr::select(mean_longitude_shifted,mean_latitude) path <- SpatialPoints(path, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs")) @@ -86,7 +93,7 @@ We can do a quick plot to just check how things look at this stage and see if th ~~~ ggplot() + - ggspatial::annotation_spatial(md_polygon, fill = "cornsilk3", size = 0) + + ggspatial::annotation_spatial(ns_polygon, fill = "cornsilk3", size = 0) + geom_point(data = path, aes(x=unlist(map(geometry,1)), y=unlist(map(geometry,2)))) + geom_path(data = path, aes(x=unlist(map(geometry,1)), y=unlist(map(geometry,2)))) + theme_void() @@ -107,14 +114,14 @@ track_pts <- st_sample(plot_path, size = 10000, type = "regular") The first `pathroutr` function we will use is `prt_visgraph`. This creates a visibility graph that connects all of the vertices for our shapefile with a Delaunay triangle mesh and removes any edges that cross land. You could think of this part as creating the viable routes an animal could swim through (marking the "water" as viable). ~~~ -vis_graph <- prt_visgraph(md_polygon, buffer = 150) +vis_graph <- prt_visgraph(ns_polygon, buffer = 100) ~~~ {: .language-r} To reroute our paths around the landmasses we will call the `prt_reroute` function. Passing `track_pts`, `md_polygon`, and `vis_graph` as arguments. To have a fully updated path we can run the `prt_update_points` function, passing our new path `track_pts_fix` with our old path `track_pts`. ~~~ -track_pts_fix <- prt_reroute(track_pts, land_barrier, vis_graph, blend = TRUE) +track_pts_fix <- prt_reroute(track_pts, ns_polygon, vis_graph, blend = TRUE) track_pts_fix <- prt_update_points(track_pts_fix, track_pts) ~~~ @@ -126,7 +133,7 @@ For `geom_point` and `geom_path` we will pass in `track_pts_fix` for the `data` ~~~ pathroutrplot <- ggplot() + - ggspatial::annotation_spatial(md_polygon, fill = "cornsilk3", size = 0) + + ggspatial::annotation_spatial(ns_polygon, fill = "cornsilk3", size = 0) + geom_point(data = track_pts_fix, aes(x=unlist(map(geometry,1)), y=unlist(map(geometry,2)))) + geom_path(data = track_pts_fix, aes(x=unlist(map(geometry,1)), y=unlist(map(geometry,2)))) + theme_void() @@ -145,6 +152,18 @@ pathroutrplot.animation <- transition_reveal(fid) + shadow_mark(past = TRUE, future = FALSE) -gganimate::animate(pathroutrplot.animation) +gganimate::animate(pathroutrplot.animation, nframes=100, detail=2) ~~~ -{: .language-r} \ No newline at end of file +{: .language-r} + +> ## A Note on Detail +> +> You'll note that the animation we've generated still crosses the landmass at certain points. This is a combination of several factors: our coastline polygon is not very high-res, our animation does not have many frames, and what frames it does have are not rendered in great detail. We can increase all of these and get a more accurate plot. For example: +> - We can specify `resolution=1` when downloading our shapefile from GADM. +> - We can increase the `nframes` variable in our call to `gganimate::animate`. +> - We can pass `detail = 2` or higher to the call to `gganimate::animate`. +> +> All of these will give us an animation that more scrupulously respects the landmass, however, they will all bloat the runtime of the code significantly. This may not be a consideration when you create your own animations, but they do make it impractical for this workshop. +> Embedded below is an animation created with high-resolution polygons and animation parameters to show an example of the kind of animation we could create with more time and processing power. +> ![High-resolution Pathroutr animation](../files/highres_pathroutr.gif) +{: .callout} \ No newline at end of file diff --git a/_episodes/15-introducing-gams.md b/_episodes/15-introducing-gams.md new file mode 100644 index 00000000..f30daccc --- /dev/null +++ b/_episodes/15-introducing-gams.md @@ -0,0 +1,402 @@ +--- +title: Spatial and Temporal Modelling with GAMs +exercises: 0 +questions: + - What are GAMs? + - How can I use GAMs to visualise my data? +--- + +GAM is short for 'Generalized Additive Model', a type of statistical model. In this lesson, we will be using GAMs to visualise our data. + +While we are still developing this lesson as templated text like our other lessons, we can provide the Powerpoint slides for the GAMs talk given by Dr. Lennox at the CANSSI Early Career Researcher workshop. You can access the slides [here](../Resources/GAMs.pptx). + +The following code is meant to be run alongside this presentation. + +~~~ +require(BTN) # remotes::install_github("robertlennox/BTN") +require(tidyverse) +require(mgcv) +require(lubridate) +require(Thermimage) +require(lunar) +library(sf) +library(gratia) + +theme_set<-theme_classic()+ + theme(text=element_text(size=20), axis.text=element_text(colour="black"))+ + theme(legend.position="top") # set a theme for ourselves + +aur<-BTN::aurland %>% + st_transform(32633) %>% + slice(2) %>% + ggplot()+ + geom_sf() + +#Load RData +load("YOUR/PATH/TO/data/otn/troutdata.RDS") + +# first plot of the data + +troutdata %>% + ggplot(aes(dt, Data, colour=Data))+ + geom_point()+ + theme_classic()+ + scale_colour_gradientn(colours=Thermimage::flirpal) + +# mapping out the data + +aur+ + geom_point(data=troutdata %>% + group_by(lon, lat) %>% + dplyr::summarise(m=mean(Data)), + aes(lon, lat, colour=m))+ + scale_colour_gradientn(colours=Thermimage::flirpal)+ + theme(legend.position="top", legend.key.width=unit(3, "cm"))+ + theme_bw() + +# 4. going for smooth +a<-troutdata %>% + mutate(h=hour(dt)) %>% + bam(Data ~ s(h, bs="cc", k=5), data=., method="fREML", discrete=T) + +b<-troutdata %>% + mutate(h=hour(dt)) %>% + bam(Data ~ s(h, bs="tp", k=5), data=., , method="fREML", discrete=T) + +c<-troutdata %>% + mutate(h=hour(dt)) %>% + bam(Data ~ h, data=., , method="fREML", discrete=T) + +tibble(h=c(0:23)) %>% + mutate(circular=predict.gam(a, newdata=.)) %>% + mutate(thin_plate=predict.gam(b, newdata=.)) %>% + # mutate(v3=predict.gam(c, newdata=.)) %>% + gather(key, value, -h) %>% + ggplot(aes(h, value, colour=key))+ + geom_line(size=2)+ + theme_set+ + labs(x="Hour", y="Predicted Depth", colour="Model")+ + scale_y_reverse(limits=c(20, 0))+ + geom_hline(yintercept=0)+ + coord_polar() + +#6. model fitting vehicle + +m1<-troutdata %>% + mutate(h=hour(dt)) %>% + mutate(foid=factor(oid)) %>% + gam(Data ~ s(h, k=5)+s(foid, bs="re"), data=., method="REML") + +m2<-troutdata %>% + mutate(h=hour(dt)) %>% + mutate(foid=factor(oid)) %>% + bam(Data ~ s(h, k=5)+s(foid, bs="re"), data=., method="fREML", discrete=T) + +tibble(h=c(0:23)) %>% + mutate(foid=1) %>% + mutate(gam=predict.gam(m1, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(bam=predict.gam(m2, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i=gam-bam) %>% + gather(key, value, -h, -foid, -i) %>% + ggplot(aes(h, value, colour=i))+ + geom_line(size=2)+ + theme_set+ + facet_wrap(~key)+ + labs(x="Hour", y="Predicted temperature", colour="Difference between predictions")+ + theme(legend.key.width=unit(3, "cm"))+ + scale_colour_gradientn(colours=Thermimage::flirpal) + + +#8. check your knots + +k1<-troutdata %>% + mutate(h=hour(dt)) %>% + bam(Data ~ s(h, bs="cc", k=5), data=., method="fREML", discrete=T) + +k2<-troutdata %>% + mutate(h=hour(dt)) %>% + bam(Data ~ s(h, bs="cc", k=15), data=., , method="fREML", discrete=T) + +tibble(h=c(0:23)) %>% + mutate("k=5"=predict.gam(k1, newdata=., type="response")) %>% + mutate("k=15"=predict.gam(k2, newdata=., type="response")) %>% + gather(key, value, -h) %>% + ggplot(aes(h, value/10, colour=key))+ + geom_line(size=2)+ + theme_set+ + labs(y="Temperature", x="Hour", colour="model") + +#9. temporal dependency + + +t1<-troutdata %>% + mutate(h=hour(dt), yd=yday(dt), foid=factor(oid)) %>% + group_by(foid, dti=round_date(dt, "1 hour")) %>% + dplyr::filter(dt==min(dt)) %>% + bam(Data ~ s(h, k=5, bs="cc")+s(yd, k=10)+s(foid, bs="re"), data=., method="fREML", discrete=T) + +t2<-troutdata %>% + mutate(h=hour(dt), yd=yday(dt), foid=factor(oid)) %>% + group_by(foid, dti=round_date(dt, "1 hour")) %>% + bam(Data ~ s(h, k=5, bs="cc")+s(yd, k=10)+s(foid, bs="re"), data=., method="fREML", discrete=T) + +expand_grid(h=c(12), + yd=c(32:60), + foid=1) %>% + mutate(partial_series=predict.gam(t1, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(full_series=predict.gam(t2, newdata=., type="response", exclude=c("s(foid)"))) %>% + gather(key, value, -h, -foid, -yd) %>% + ggplot(aes(yd, value, colour=key))+ + geom_point(data=troutdata %>% + mutate(h=hour(dt), yd=yday(dt), foid=factor(oid)), + aes(yday(dt), Data), inherit.aes=F)+ + geom_path(size=2)+ + theme_set+ + labs(x="Date", y="Temperature", colour="Model") + +# 10. spatial dependency + +aur+ + geom_point(data=troutdata %>% + group_by(lon, lat) %>% + dplyr::summarise(m=mean(Data)), + aes(lon, lat, colour=m))+ + scale_colour_gradientn(colours=Thermimage::flirpal)+ + theme(legend.position="top", legend.key.width=unit(3, "cm"))+ + theme_bw()+ + theme_set+ + theme(legend.position="top", legend.key.width=unit(3, "cm"))+ + labs(colour="mean temperature") + +#11. interactions + +mi<-troutdata %>% + mutate(h=hour(dt), yd=yday(dt), foid=factor(oid)) %>% + bam(Data ~ te(h, yd, bs=c("cc", "tp"), k=c(5, 10))+ + s(foid, bs="re"), data=., family=Gamma(link="log"), method="fREML", discrete=T) + +ms<-troutdata %>% + mutate(h=hour(dt), yd=yday(dt), foid=factor(oid)) %>% + bam(Data ~ s(h, bs="cc", k=5)+ + s(yd, bs="tp", k=10)+ + s(foid, bs="re"), data=., family=Gamma(link="log"), method="fREML", discrete=T) + + +p1<-expand_grid(h=c(0:23), yd=c(182:212)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(ms, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i="Simple model, AIC=426 801") %>% + ggplot(aes(yd, h, fill=value))+ + geom_raster()+ + scale_fill_viridis_c()+theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + labs(x="Date", y="Hour", fill="Predicted temperature")+ + facet_wrap(~i) + +p2<-expand_grid(h=c(0:23), yd=c(182:212)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(mi, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i="Interaction model, AIC=425 805") %>% + ggplot(aes(yd, h, fill=value))+ + geom_raster()+ + scale_fill_viridis_c()+theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + labs(x="Date", y="Hour", fill="Predicted temperature")+ + facet_wrap(~i) + +AIC(mi, ms) +gridExtra::grid.arrange(p1, p2) + +# is it the moon? + +a<-expand_grid(h=c(0:23), yd=c(182:212)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(ms, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i="Simple model, AIC= -262 975") %>% + mutate(dt=ymd("2022-12-31")+days(yd)) %>% + mutate(l=lunar::lunar.illumination(dt)) %>% + ggplot(aes(yd, h, fill=value))+ + geom_raster()+ + scale_fill_viridis_c()+theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + labs(x="Date", y="Hour", fill="Predicted temperature")+ + facet_wrap(~i)+ + geom_point(data=expand_grid(h=c(0:23), yd=c(182:212)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(ms, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i="Simple model, AIC= -262 975") %>% + mutate(dt=ymd("2022-12-31")+days(yd)) %>% + mutate(l=lunar::lunar.illumination(dt)) %>% + distinct(dt, yd, l), + aes(yd, 10, size=l), inherit.aes=F, colour="white") + +b<-expand_grid(h=c(0:23), yd=c(182:212)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(mi, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i="Interaction model, AIC=425 805") %>% + ggplot(aes(yd, h, fill=value))+ + geom_raster()+ + scale_fill_viridis_c()+theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + labs(x="Date", y="Hour", fill="Predicted temperature")+ + facet_wrap(~i)+ + geom_point(data=expand_grid(h=c(0:23), yd=c(182:212)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(ms, newdata=., type="response", exclude=c("s(foid)"))) %>% + mutate(i="Simple model, AIC=426 801") %>% + mutate(dt=ymd("2022-12-31")+days(yd)) %>% + mutate(l=lunar::lunar.illumination(dt)) %>% + distinct(dt, yd, l), + aes(yd, 10, size=l), inherit.aes=F, colour="white") + + +gridExtra::grid.arrange(a, b) + + +########### the worked example + +troutdata %>% + mutate(lun=lunar::lunar.illumination(dt)) %>% + ggplot(aes(lun, Data))+ + geom_point()+ + theme_set+ + labs(x="Lunar illumination", y="Temperature") + +troa<-troutdata %>% + mutate(foid=factor(oid)) %>% + mutate(lun=lunar::lunar.illumination(dt)) + +m0<-troa %>% + bam(Data ~ lun, data=., family=Gamma(link="log")) + +tibble(lun=seq(0,1, by=0.1)) %>% + mutate(p=predict.gam(m0, newdata=., type="response")) %>% + ggplot(aes(lun, p))+ + geom_point(data=troa, aes(lun, Data))+ + geom_line(colour="red")+ + theme_set+ + labs(x="Moonlight", y="Temperature") + +m0<-troa %>% + bam(Data ~ lun+s(foid, bs="re"), data=., family=Gamma(link="log")) + +tibble(lun=seq(0,1, by=0.1)) %>% + mutate(foid=1) %>% + mutate(p=predict.gam(m0, newdata=., type="response", exclude=c("s(foid)"))) %>% + ggplot(aes(lun, p))+ + geom_point(data=troa, aes(lun, Data, colour=factor(foid)))+ + geom_line(colour="red")+ + theme_set+ + guides(colour=F)+ + labs(x="Moonlight", y="Temperature") + +m0<-troa %>% + bam(Data ~ s(lun, k=7)+s(foid, bs="re"), data=., family=Gamma(link="log")) + +tibble(lun=seq(0,1, by=0.1)) %>% + mutate(foid=1) %>% + mutate(p=predict.gam(m0, newdata=., type="response", exclude=c("s(foid)"))) %>% + ggplot(aes(lun, p))+ + geom_point(data=troa, aes(lun, Data, colour=factor(foid)))+ + geom_line(colour="red")+ + theme_set+ + guides(colour=F)+ + labs(x="Moonlight", y="Temperature") + +m01<-troa %>% + mutate(yd=yday(dt), h=hour(dt)) %>% + bam(Data ~ s(lun, k=7)+s(foid, bs="re")+ + s(h, bs="cc", k=5)+ + s(lon, lat, k=15), data=., family=Gamma(link="log")) + +BTN::aurland %>% + st_transform(32633) %>% + slice(2) + +tibble(x=c(7.26, 7.32), + y=c(60.85, 60.87)) %>% + st_as_sf(., coords=c("x", "y")) %>% + st_set_crs(4326) %>% + st_transform(32633) + +sp<-expand_grid(lon=seq(80077, 83583, by=10), + lat=seq(6770907, 6772740, by=10)) %>% + st_as_sf(., coords=c("lon", "lat")) %>% + st_set_crs(32633) %>% + st_intersection(BTN::aurland %>% + slice(1) %>% + st_transform(32633)) %>% + as(., "Spatial") %>% + as_tibble %>% + dplyr::rename(lon=coords.x1, lat=coords.x2) %>% + dplyr::select(lon, lat) + +sp %>% + slice(1) %>% + expand_grid(., lun=seq(0, 1, by=0.1), h=c(0:23)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(m01, newdata=., exclude=c("s(foid)"))) %>% + ggplot(aes(h, lun, fill=value))+ + geom_raster()+ + scale_fill_gradientn(colours=Thermimage::flirpal)+ + theme_classic()+ + theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + labs(x="Hour", y="Moonlight", fill="Predicted Temperature") + +expand_grid(sp, lun=seq(0, 1, by=0.1), h=c(0, 12), yd=200) %>% + mutate(foid=1) %>% + mutate(p=predict.gam(m01, newdata=., type="response", exclude=c("s(foid)", "s(lon,lat)", "s(h)"))) %>% + ggplot(aes(lun, p))+ + geom_point(data=troa, aes(lun, Data, colour=factor(foid)))+ + geom_line()+ + theme_set+ + guides(colour=F)+ + labs(x="Moonlight", y="Temperature") + + +sp %>% + expand_grid(., lun=seq(0, 1, by=0.3)) %>% + mutate(foid=1, h=1) %>% + mutate(value=predict.gam(m01, newdata=., type="response", exclude=c("s(foid)"))) %>% + ggplot(aes(lon, lat, fill=value))+ + geom_raster()+ + scale_fill_gradientn(colours=Thermimage::flirpal)+ + theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + labs(x="UTM (x)", y="UTM (y)", fill="Predicted Temperature")+ + coord_fixed(ratio=1)+ + facet_wrap(~lun) + +gratia::draw(m01) + +sp %>% + expand_grid(., lun=seq(0, 1, by=0.05), h=c(0:23)) %>% + mutate(foid=1) %>% + mutate(value=predict.gam(m01, newdata=., type="response", exclude=c("s(foid)", "s(lon,lat)"))) %>% + ggplot(aes(lun, value, colour=h))+ + geom_point()+ + theme_set+ + theme(legend.key.width=unit(3, "cm"))+ + scale_fill_gradientn(colours=Thermimage::flirpal) + +coef(m01) %>% + data.frame %>% + rownames_to_column() %>% + as_tibble %>% + dplyr::filter(grepl("foid", rowname)) %>% + bind_cols(troa %>% distinct(foid)) %>% + left_join(troutdata %>% + group_by(oid) %>% + dplyr::summarise(m=mean(Data)) %>% + dplyr::rename(foid=oid) %>% + mutate(foid=factor(foid))) %>% + dplyr::rename(value=2) %>% + ggplot(aes(reorder(foid, value), value, colour=m))+ + geom_point()+ + coord_flip()+ + theme_set+ + labs(x="ID", y="Random intercept of temperature", size="Length (mm)", colour="True mean temperature")+ + scale_colour_viridis_c() +~~~ +{: .language-r} \ No newline at end of file diff --git a/_episodes/16-introducing-miscYaps.md b/_episodes/16-introducing-miscYaps.md new file mode 100644 index 00000000..cfd87f17 --- /dev/null +++ b/_episodes/16-introducing-miscYaps.md @@ -0,0 +1,70 @@ +--- +title: Introduction to the miscYAPS package +teaching: 0 +exercises: 0 +questions: + - What is YAPS? + - For what applications is YAPS well-suited? + - How can I use the miscYAPS package to make working with YAPS easier? +--- + +[YAPS](https://github.com/baktoft/yaps) (short for 'Yet Another Position Solver') is a package originally presented in a 2017 [paper](https://nature.com/articles/s41598-017-14278-z.pdf) by Baktoft, Gjelland, Økland & Thygesen. YAPS represents a way of positioning aquatic animals using statistical analysis (specifically, maximum likelihood analysis) in conjunction with time of arrival data, movement models, and state-space models. Likewise, [`miscYaps()`](https://github.com/robertlennox/miscYAPS) is a package developed by Dr. Robert Lennox to provide more intuitive wrappers for YAPS functions, to ease the analysis process. + +While we are still developing this lesson as templated text, as with the other lessons, we can provide the Powerpoint slides for the miscYaps talk given by Dr. Lennox at the CANSSI Early Career Researcher workshop. You can access the slides [here](../Resources/YAPS.pptx). + +The following code is meant to be run alongside this powerpoint. + +~~~ +remotes::install_github("robertlennox/miscYAPS") +require(yaps) +require(miscYAPS) + +remotes::install_github("robertlennox/BTN") +require(BTN) + +require(tidyverse) +require(lubridate) +require(data.table) + +data(boats) +dets<-boats %>% pluck(1) +hydros<-dets %>% + dplyr::distinct(serial, x, y, sync_tag=sync) %>% + mutate(idx=c(1:nrow(.)), z=1) +detections<-dets %>% + dplyr::filter(Spp=="Sync") %>% + dplyr::select(ts, epo, frac, serial, tag) +ss_data<-boats %>% pluck(2) %>% + dplyr::rename(ts=dt) %>% + setDT + +############ +require(miscYAPS) +sync_model<-sync(hydros, + detections, + ss_data, + keep_rate=0.5, + HOW_THIN=100, + ss_data_what="data", + exclude_self_detections=T, + fixed=NULL) +plotSyncModelResids(sync_model) +sync_model<-sync(hydros, + detections, + ss_data, + keep_rate=0.5, + HOW_THIN=100, + ss_data_what="data", + exclude_self_detections=T, + fixed=c(1:9, 11:20)) +fish_detections<-dets %>% + dplyr::filter(Spp!="Sync") %>% + mutate(tag=factor(tag)) %>% + dplyr::select(ts, epo, frac, tag, serial) +tr<-swim_yaps(fish_detections, runs=3, rbi_min=60, rbi_max=120) +data(aur) +btnstorel <- BTN::storel +raster::plot(btnstorel) +points(tr$x, tr$y, pch=1) +~~~ +{: .language-r} \ No newline at end of file diff --git a/_episodes/15-using-gitlab.md b/_episodes/17-using-gitlab.md similarity index 100% rename from _episodes/15-using-gitlab.md rename to _episodes/17-using-gitlab.md diff --git a/_episodes/16-other-curriculums.md b/_episodes/18-other-curriculums.md similarity index 100% rename from _episodes/16-other-curriculums.md rename to _episodes/18-other-curriculums.md diff --git a/files/highres_pathroutr.gif b/files/highres_pathroutr.gif new file mode 100644 index 00000000..b7456ee7 Binary files /dev/null and b/files/highres_pathroutr.gif differ