analysis/gather_CDS_wind_data.R

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

# download the latest development version
# of github!!!!!
library(ecmwfr)

# format a base CDS request for ERA-5 data
# this request will be altered to suite different
# sampling locations, keeping most variables
# static, only altering time / space parameters
base_request <- list(
  'dataset' = 'reanalysis-era5-pressure-levels',
  'product_type' = 'reanalysis',
  'format' = 'netcdf',
  'variable' = c('u_component_of_wind','v_component_of_wind'), # + temperature c('u_component_of_wind','v_component_of_wind', 'temperature')
  'pressure_level' = '850',
  'year' = '2016',
  'month' = '01',
  'day' = '01',
  'time' = '00:00',
  "area" = "51/-0/50/10",
  'target' = 'download.nc'
)

# 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)

# set the different ERA-5 pressure levels to pick from
levels <- c(1,2,3,5,7,10,20,30,50,70,100,125,150,
  175,200,225,250,300,350,400,450,
  500,550,600,650,700,750,775,800,
  825,850,875,900,925,950,975,1000)

# select the closest match in ERA-5 pressure level
# for a given altitude (pressure)
df$level <- do.call("rbind", lapply(df$pressure, function(pressure){
  abs_diff <- abs(levels - pressure)
  levels[which(abs_diff == min(abs_diff))]
}))

# now loop over all the locations altering the base_request string
# for set variables, in particular:
# - location (padded by a degree in lat lon)
# - the time
# - pressure level (altitude)
# - date
# of the query

apply(df, 1, function(x){

  # split out variables for readability
  # formatting purposes
  year <- sprintf("20%02d",x['year'])
  month <- sprintf("%02d",x['month'])
  day <- sprintf("%02d",x['day'])
  time <- sprintf("%02d:00",x['hour'])
  hour <- sprintf("%02d",x['hour'])
  pressure_level <- as.character(x['level'])
  lat <- x['latitude']
  lon <- x['longitude']
  target <- paste0(year,"_",month,"_",day,"_",hour,".nc")

  # format the region of interest
  area <- paste0(lat + 0.5,"/",lon - 0.5,"/",lat - 0.5,"/", lon + 0.5)

  # modify the request
  new_request <- wf_modify_request(base_request,
                    year = year,
                    month = month,
                    day = day,
                    time = time,
                    pressure_level = pressure_level,
                    area = area,
                    target = target)

  wf_request(new_request,
             user = "2088",
             transfer = TRUE,
             path = "data/wind/")
  Sys.sleep(1)
})