analysis/compile_wind_data.R

In khufkens/swiftr: Process bird (swift) GPS position tracker data and meta-data.

library(raster)
library(swiftr)
library(maps)
library(mapview)
library(sf)

# read position data, remove 0 values
df <- sr_read_pos("data-raw/Obs010617_050807_Tag13819 - Copy - Combined.pos")
df <- df[which(df$latitude != 0),]

# convert altitude to pressure in hPa
df$pressure <- meter_to_hPa(df$altitude)
df <- df[1:150,]

# calculate heading + distance
out <- do.call("rbind", lapply(seq(1,nrow(df)), function(i){
  values <- try(great_circle_distance(df$latitude[i],
                            df$longitude[i],
                            df$latitude[i+1],
                            df$longitude[i+1]))

  # format the date of starting position
  date1 <- strptime(sprintf("20%02d-%s-%s %s:%s",
                           df$year[i],
                           df$month[i],
                           df$day[i],
                           df$hour[i],
                           df$min[i]), "%Y-%m-%d %H:%M")

  # set end position
  date2 <- strptime(sprintf("20%02d-%s-%s %s:%s",
                           df$year[i+1],
                           df$month[i+1],
                           df$day[i+1],
                           df$hour[i+1],
                           df$min[i+1]), "%Y-%m-%d %H:%M")

  # calculate difference in seconds
  seconds <- as.numeric(difftime(date2,date1, units = "secs"))

  # calculate speed (distance in meters / seconds)
  speed <- values$distance/seconds

  # grab local wind speed from ERA-5 data
  filename <- sprintf("data/wind/20%02d_%02d_%02d_%02d.nc",
                      df$year[i],
                      df$month[i],
                      df$day[i],
                      df$hour[i])

  # extract the value of wind speed and direction from
  # the raster for the start location
  m <- wind_direction(filename)
  location <- SpatialPoints(cbind(
    df$longitude[i], df$latitude[i]),
    CRS("+init=epsg:4326"))
  wind <- extract(m, location)

  # return everything
  return(data.frame(values, speed, seconds, wind_speed = wind[1],
                    wind_direction = wind[2]))
}))

# sort wind directions
out$direction <- ifelse(out$heading > 90 & out$heading < 270,
                         "forward",
                         "backward")

# bind it in a dataframe
df <- cbind(df,out)

# convert to an sf object TODO rename variable
test <- st_as_sf(df, coords = c("longitude","latitude"),
                 crs = 4326)

points <- test$geometry

# Number of total linestrings to be created
n <- length(points) - 1

# Build linestrings for the complete flight path
linestrings <- lapply(X = 1:n, FUN = function(x) {
  pair <- st_combine(c(points[x], points[x + 1]))
  line <- st_cast(pair, "LINESTRING")
  return(line)
})
line <- st_multilinestring(do.call("rbind", linestrings))
line <- st_sfc(line, crs = "+proj=longlat +datum=WGS84")

# combine everything in an mapview visualization
m <- mapview(line) + mapview(test, zcol = "direction",
                        popup = popupTable(test))
print(m)