library(tidyverse)
library(sf)
library(mapdeck)
library(lubridate)
library(data.table)

mb_key <-'your_key_goes_here'

# my key is in this file that has not been uploaded to Github
source("scripts/mapbox_key.R")

# reading in the geojson file
heli_big2 <- read_sf("data/heli_points.geojson")

glimpse(heli_big2)
## Observations: 30,446
## Variables: 26
## $ timestamp                 <dttm> 2020-05-31 20:00:01, 2020-05-31 20:00:04, …
## $ altitude                  <int> 1025, 1025, 1025, 1000, 1000, 975, 975, 975…
## $ speed                     <dbl> 44.1, 43.9, 44.8, 46.5, 48.9, 48.1, 45.1, 4…
## $ track.x                   <dbl> 123.0, 133.2, 141.3, 151.8, 157.1, 163.1, 1…
## $ vert_rate                 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ baro_rate                 <int> 128, 64, 0, -192, -192, -128, -128, 0, 0, 0…
## $ code.x                    <chr> "A1DFFB", "A1DFFB", "A1DFFB", "A1DFFB", "A1…
## $ join_id                   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, …
## $ type                      <chr> NA, "adsb_icao", NA, NA, NA, "adsb_icao", N…
## $ flight                    <chr> NA, "N22PP", NA, NA, NA, "N22PP", NA, NA, N…
## $ alt_geom                  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ track.y                   <dbl> NA, 133.16, NA, NA, NA, 163.08, NA, NA, NA,…
## $ code.y                    <chr> NA, "A1DFFB", NA, NA, NA, "A1DFFB", NA, NA,…
## $ date                      <date> NA, 2020-06-01, NA, NA, NA, 2020-06-01, NA…
## $ start_time                <dbl> 1590969602, 1590969602, 1590969602, 1590969…
## $ flightnum                 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ timestr                   <dttm> 2020-05-31 20:00:01, 2020-05-31 20:00:04, …
## $ altitude_modeled          <dbl> 1286.176, 1286.176, 1286.176, 1261.084, 126…
## $ altitude_adjusted         <dbl> 1284.18, 1284.18, 1284.18, 1259.18, 1259.18…
## $ altitude_modeled_meters   <dbl> 392.0265, 392.0265, 392.0265, 384.3785, 384…
## $ altitude_adjusted_meters  <dbl> 391.4182, 391.4182, 391.4182, 383.7982, 383…
## $ bare_earth_elevation      <dbl> 26.628, 26.040, 25.739, 26.902, 27.806, 24.…
## $ dsm_elevation             <dbl> 50.06, 38.16, 54.23, 26.94, 34.44, 24.66, 2…
## $ bare_earth_elevation_feet <dbl> 87.362208, 85.433077, 84.445542, 88.261159,…
## $ dsm_elevation_feet        <dbl> 164.238855, 125.196854, 177.919952, 88.3858…
## $ geometry                  <POINT [°]> POINT (-77.02939 38.90662), POINT (-7…
# Can choose which helicopters to focus on and what time
# This code below will filter only the two National Guard Lakotas
# And show where they traveled between 9 and 10 p.m. on June 1
# AE0BED is the Black Hawk, if you'd like to include that

heli_one9 <- heli_big2 %>% 
  #filter(code.x=="AE1F45" | code.x=="AE1FE3" | code.x=="AE0BED") %>% 
  filter(code.x=="AE1F45" | code.x=="AE1FE3") %>% 
  mutate(day=day(timestamp)) %>% 
  mutate(hour=hour(timestamp)) %>% 
  mutate(minute=minute(timestamp)) %>% 
  filter(day==1) %>% 
  filter(hour==21) %>% 
  filter(minute>00 & minute <=59) 

heli_one10 <- heli_big2 %>% 
  #filter(code.x=="AE1F45" | code.x=="AE1FE3" | code.x=="AE0BED") %>% 
  filter(code.x=="AE1F45" | code.x=="AE1FE3") %>% 
  mutate(day=day(timestamp)) %>% 
  mutate(hour=hour(timestamp)) %>% 
  mutate(minute=minute(timestamp)) %>% 
  filter(day==1) %>% 
  filter(hour==22) %>% 
  filter(minute>00 & minute <=59) 

heli_one <- rbind(heli_one9, heli_one10) %>% 
  arrange(code.x,timestamp)

## Quick chart for altitude

heli_one <- heli_one %>% 
  mutate(heli=case_when(
    code.x=="AE1FE3" ~ "NG Lakota",
    code.x=="AE1F45" ~ "NG Lakota Red Cross",
    code.x=="AE0BED" ~ "Black Hawk"
  ))

ggplot(heli_one, aes(x=timestamp, y=altitude_adjusted, colour=heli)) +
  geom_line() +
  theme_minimal() +
  labs(title="How low the two Lakotas traveled in DC on June 1")

## Mapping the path with Mapbox/Mapdeck
## We need to take the flight data data frame and
## transform the coordinate sinto a sfc_LINESTRING geometry 
## so it will be compatible with Mapbox/Mapdeck code

heli_path <- heli_one 


## grab the helicopter path coordinates only
dt <- sf::st_coordinates( heli_path )

## transform the matrix into a data table
dt <- as.data.table( dt )


# Extracting the altitude into an array
hel_alt <- heli_one %>% 
  mutate(day=day(timestamp)) %>% 
  pull(altitude)

# Extracting the timestamps into an array
hel_time <- heli_one %>% 
  mutate(day=day(timestamp)) %>% 
  pull(timestamp)

# Extracting a numeric value for the helicopter ID (1 and 2)
hel_id <- heli_one %>% 
  mutate(day=day(timestamp)) %>% 
  mutate(id=case_when(
    code.x=="AE1F45" ~ 1,
    code.x=="AE1FE3" ~ 2,
   # code.x=="AE0BED" ~ 3, # in case there was interest in mapping the Black Hawk
    TRUE ~ 0
  )) %>% 
  pull(id)


## This is a lot of work just to structure our 
## lat, lon, altitude, timestamp, and helicopter id
## juuust right

## Need Z and M attributes
dt[, Z := hel_alt ]   ## setting Z as our altitude
dt[, M := hel_time]   ## setting M as our timestamp
dt[, L1 := hel_id]    ## setting L1 as our helicopter id number

## now we can convert back to sf LINESTRING
heli_trips <- sfheaders::sf_linestring(
  obj = dt
  , x = "X"
  , y = "Y"
  , z = "Z"
  , m = "M"
  , linestring_id = "L1"
)

Mapbox visualization

mapdeck(
  token=mb_key # if you don't have an mb_key you won't see any street tiles
  , location = c(-77.024058, 38.904742)
  , zoom = 13
  , pitch = 200
  , style = mapdeck_style("dark")
) %>%
  add_trips(
    data = heli_trips
    , stroke_colour = "L1"
    , stroke_width = 30
    , animation_speed = 50
    , trail_length = 1000
  )

Click and drag the map around or zoom in. You can also hold control or command to tilt the camera and get a better glimpse at the elevation of the paths.