harp news… and some features you might not know about

Andrew Singleton

MET Norway

New read function

read_forecast()

read_forecast(
  start_date,
  end_date,
  fcst_model,
  parameter
)

Replaces read_eps_interpolate() and read_det_interpolate()

…and does a whole lot more!

Default behaviour

Read data “as is”

  • For vfld, this means to just read the data
  • For gridded data, the gridded fields are returned

vfld

read_forecast(
  start_date  = 2019021700,
  end_date    = 2019021718,
  fcst_model  = "AROME_Arctic_prod",
  file_path   = system.file("vfld", package = "harpData"),
  return_data = TRUE
)
## â—Ź AROME_Arctic_prod
## # A tibble: 65,004 x 5
##    fcdate              lead_time   SID parameter AA_det
##    <dttm>                  <dbl> <int> <chr>      <dbl>
##  1 2019-02-17 00:00:00         0  1001 D10m      215.  
##  2 2019-02-17 00:00:00         0  1002 D10m      335.  
##  3 2019-02-17 00:00:00         0  1003 D10m       14.9 
##  4 2019-02-17 00:00:00         0  1004 D10m       11.1 
##  5 2019-02-17 00:00:00         0  1006 D10m      238.  
##  6 2019-02-17 00:00:00         0  1007 D10m        9.70
##  7 2019-02-17 00:00:00         0  1008 D10m      348.  
##  8 2019-02-17 00:00:00         0  1009 D10m      344.  
##  9 2019-02-17 00:00:00         0  1010 D10m      346.  
## 10 2019-02-17 00:00:00         0  1011 D10m        6.67
## # … with 64,994 more rows

grib

read_forecast(
  2018071000,
  2018071000,
  "AROME_Arctic",
  "T2m",
  file_path     = system.file("grib/AROME_Arctic", package = "harpData"),
  file_template = "harmonie_grib_fp",
  return_data   = TRUE
)
## â—Ź AROME_Arctic
## # A tibble: 2 x 4
##   fcdate              lead_time parameter AA_det                
##   <dttm>                  <dbl> <chr>     <geolist>             
## 1 2018-07-10 00:00:00         0 T2m       <geofield [100 Ă— 100]>
## 2 2018-07-10 00:00:00         3 T2m       <geofield [100 Ă— 100]>

Transformations

  • interpolate : interpolate gridded data to points
  • regrid : interpolate / reproject gridded data to a new grid
  • xsection : extract a vertical cross section of the data

transformation = "interpolate"

read_forecast(
    2018071000,
    2018071000,
    "AROME_Arctic",
    "T2m",
    transformation = "interpolate",
    file_path      = system.file("grib/AROME_Arctic", package = "harpData"),
    file_template  = "harmonie_grib_fp",
    return_data    = TRUE
)
## â—Ź AROME_Arctic
## # A tibble: 14 x 5
##   fcdate              lead_time   SID parameter AA_det
##   <dttm>                  <dbl> <int> <chr>      <dbl>
## 1 2018-07-10 00:00:00         0  1003 T2m         277.
## 2 2018-07-10 00:00:00         0  1004 T2m         278.
## 3 2018-07-10 00:00:00         0  1006 T2m         278.
## 4 2018-07-10 00:00:00         0  1007 T2m         278.
## 5 2018-07-10 00:00:00         0  1008 T2m         280.
## 6 2018-07-10 00:00:00         0  1017 T2m         277.
## 7 2018-07-10 00:00:00         0 20107 T2m         278.
## # … with 7 more rows

transformation_opts = interpolate_opts()

my_stations <- data.frame(
  SID = c(1003, 1004, 1006),
  lat = c(77.0000, 78.9167, 78.2506),
  lon = c(15.5000, 11.9331, 22.8225)
)
read_forecast(
  2018071000,
  2018071000,
  "AROME_Arctic",
  "T2m",
  transformation      = "interpolate",
  transformation_opts = interpolate_opts(
    stations    = my_stations,
    correct_t2m = FALSE,
    method      = "bilinear"
  ),
  file_path           = system.file("grib/AROME_Arctic", package = "harpData"),
  file_template       = "harmonie_grib_fp",
  return_data         = TRUE
)

transformation_opts = interpolate_opts()

## â—Ź AROME_Arctic
## # A tibble: 6 x 5
##   fcdate              lead_time   SID parameter AA_det
##   <dttm>                  <dbl> <dbl> <chr>      <dbl>
## 1 2018-07-10 00:00:00         0  1003 T2m         277.
## 2 2018-07-10 00:00:00         0  1004 T2m         277.
## 3 2018-07-10 00:00:00         0  1006 T2m         278.
## 4 2018-07-10 00:00:00         3  1003 T2m         277.
## 5 2018-07-10 00:00:00         3  1004 T2m         276.
## 6 2018-07-10 00:00:00         3  1006 T2m         278.

transformation = "regrid"

library(meteogrid)
my_domain <- Make.domain(
  "lambert",
  clonlat = c(15.57905, 78.21638),
  nxny    = c(15, 15), 
  dxdy    = c(10000, 10000),
  reflat  = 77.5,   
  reflon  = -25
)
read_forecast(
  2018071000,
  2018071000,
  "AROME_Arctic",
  "S10m",
  transformation      = "regrid",
  transformation_opts = regrid_opts(new_domain = my_domain),
  file_path           = system.file("grib/AROME_Arctic", package = "harpData"),
  file_template       = "harmonie_grib_fp",
  return_data         = TRUE
)

transformation = "regrid"

## â—Ź AROME_Arctic
## # A tibble: 2 x 4
##   fcdate              lead_time parameter AA_det              
##   <dttm>                  <dbl> <chr>     <geolist>           
## 1 2018-07-10 00:00:00         0 S10m      <geofield [15 Ă— 15]>
## 2 2018-07-10 00:00:00         3 S10m      <geofield [15 Ă— 15]>

transformation = "xsection"

read_forecast(
  20180710, 
  20180710, 
  "AAEPS", 
  "T", 
  members             = c(0, 1), 
  vertical_coordinate = "model", 
  transformation      = "xsection", 
  transformation_opts = xsection_opts(a = c(11, 78), b = c(14, 78.25))
  file_path           = system.file("netcdf", package = "harpData"), 
  file_template       = "{fcst_model}/fc{YYYY}{MM}{DD}{HH}.nc", 
  file_format_opts    = netcdf_opts("met_norway_eps"), 
  return_data         = TRUE
)

transformation = "xsection"

## â—Ź AAEPS
## # A tibble: 130 x 7
##    fcdate              lead_time parameter level_type   level AAEPS_mbr000
##    <dttm>                  <dbl> <chr>     <chr>        <dbl> <list>      
##  1 2018-07-10 00:00:00         0 T         model      0.00987 <dbl [30]>  
##  2 2018-07-10 00:00:00         0 T         model      0.0296  <dbl [30]>  
##  3 2018-07-10 00:00:00         0 T         model      0.0494  <dbl [30]>  
##  4 2018-07-10 00:00:00         0 T         model      0.0691  <dbl [30]>  
##  5 2018-07-10 00:00:00         0 T         model      0.0890  <dbl [30]>  
##  6 2018-07-10 00:00:00         0 T         model      0.109   <dbl [30]>  
##  7 2018-07-10 00:00:00         0 T         model      0.129   <dbl [30]>  
##  8 2018-07-10 00:00:00         0 T         model      0.149   <dbl [30]>  
##  9 2018-07-10 00:00:00         0 T         model      0.169   <dbl [30]>  
## 10 2018-07-10 00:00:00         0 T         model      0.190   <dbl [30]>  
## # … with 120 more rows, and 1 more variable: AAEPS_mbr001 <list>

Save to SQLite

(Interpolated / point data only)

read_forecast(
  start_date       = 2019021700,
  end_date         = 2019021718,
  fcst_model       = c("MEPS_prod", "CMEPS_prod"),
  parameter        = c("T", "Td"),
  by               = 3,
  lead_time        = seq(0, 12, 3),
  members          = list(
    MEPS_prod  = seq(0, 10),
    CMEPS_prod = c(0, 1, 3, 4, 5, 6)
  ),
  lags             = list(CMEPS_prod = c(0, 0, 2, 2, 1, 1)),
  file_path        = system.file("vfld", package = "harpData"),
  file_template    = "vfld_eps",
  output_file_opts = sqlite_opts(
    path = "/path/to/sqlite/data", template = "fctable"
  )
)

Vertical profile verification

fcst <- read_point_forecast(
  start_date          = 2019021700,
  end_date            = 2019021718,
  fcst_model          = c("CMEPS_prod", "MEPS_prod"),
  fcst_type           = "EPS",
  parameter           = "T", 
  vertical_coordinate = "pressure",
  lead_time           = seq(0, 12, 3),
  by                  = "6h",
  lags                = list(
    CMEPS_prod = paste0(seq(0, 5), "h"),
    MEPS_prod  = "0h"
  ),
  file_path           = system.file("FCTABLE/ensemble", package = "harpData")
) %>% common_cases()

Vertical profile verification

obs <- read_point_obs(
  first_validdate(fcst),
  last_validdate(fcst),
  "T",
  vertical_coordinate = "pressure",
  obs_path            = system.file("OBSTABLE", package = "harpData")
)

fcst <- join_to_fcst(fcst, obs)

verif <- ens_verify(fcst, T, groupings = c("leadtime", "p"))

plot_profile_verif()

plot_profile_verif(verif, spread_skill, filter_by = vars(leadtime == 3))

plot_profile_verif()

verif <- mutate_list(verif, leadtime = paste0("T + ", leadtime, "h"))
plot_profile_verif(verif, spread_skill, facet_by = vars(fct_inorder(leadtime)))

Grouped verification

fcst <- read_point_forecast(
  2019021700,
  2019021718,
  "MEPS_prod",
  "EPS",
  "T2m",
  file_path = system.file("FCTABLE/ensemble", package = "harpData")
)

obs <- read_point_obs(
  first_validdate(fcst),
  last_validdate(fcst),
  "T2m",
  stations = pull_stations(fcst),
  obs_path = system.file("OBSTABLE", package = "harpData")
)

Grouped verification

obs  <- mutate(obs, region = floor(SID / 1000))

fcst <- join_to_fcst(fcst, obs)

verif <- ens_verify(
  fcst, T2m, groupings = list("leadtime", c("leadtime", "region"))
)

verif <- mutate_list(
  verif, region = fct_inorder(paste0("Region: ", region))
)

Grouped verification

plot_point_verif(verif, crps, facet_by = vars(region))

Binding verifcation data

verif <- bind_point_verif(verif_T2m, verif_S10m)
plot_point_verif(verif, spread_skill, facet_by = vars(parameter), facet_scales = "free_y")

Gridded data

t2m <- read_forecast(
  20180710, 
  20180710, 
  "AAEPS", 
  "T2m", 
  lead_time           = seq(0, 12, 3),
  file_path           = system.file("netcdf", package = "harpData"), 
  file_template       = "{fcst_model}/fc{YYYY}{MM}{DD}{HH}.nc", 
  file_format_opts    = netcdf_opts("met_norway_eps"), 
  return_data         = TRUE
) %>% 
  scale_point_forecast(-273.15, new_units = "degC")

Gridded data

t2m
## â—Ź AAEPS
## # A tibble: 5 x 13
##   fcdate              lead_time parameter validdate           level_type level
##   <dttm>                  <dbl> <chr>     <dttm>              <chr>      <dbl>
## 1 2018-07-10 00:00:00         0 T2m       2018-07-10 00:00:00 height         2
## 2 2018-07-10 00:00:00         3 T2m       2018-07-10 03:00:00 height         2
## 3 2018-07-10 00:00:00         6 T2m       2018-07-10 06:00:00 height         2
## 4 2018-07-10 00:00:00         9 T2m       2018-07-10 09:00:00 height         2
## 5 2018-07-10 00:00:00        12 T2m       2018-07-10 12:00:00 height         2
## # … with 7 more variables: units <chr>, fcst_cycle <chr>,
## #   AAEPS_mbr000 <geolist>, AAEPS_mbr001 <geolist>, AAEPS_mbr002 <geolist>,
## #   AAEPS_mbr003 <geolist>, AAEPS_mbr004 <geolist>

Gridded data

ens_stats(t2m)
## â—Ź AAEPS
## # A tibble: 5 x 15
##   fcdate              lead_time parameter validdate           level_type level
##   <dttm>                  <dbl> <chr>     <dttm>              <chr>      <dbl>
## 1 2018-07-10 00:00:00         0 T2m       2018-07-10 00:00:00 height         2
## 2 2018-07-10 00:00:00         3 T2m       2018-07-10 03:00:00 height         2
## 3 2018-07-10 00:00:00         6 T2m       2018-07-10 06:00:00 height         2
## 4 2018-07-10 00:00:00         9 T2m       2018-07-10 09:00:00 height         2
## 5 2018-07-10 00:00:00        12 T2m       2018-07-10 12:00:00 height         2
## # … with 9 more variables: units <chr>, fcst_cycle <chr>,
## #   AAEPS_mbr000 <geolist>, AAEPS_mbr001 <geolist>, AAEPS_mbr002 <geolist>,
## #   AAEPS_mbr003 <geolist>, AAEPS_mbr004 <geolist>, ens_mean <geolist>,
## #   ens_spread <geolist>

Gridded data

ens_stats(t2m, min = TRUE, max = TRUE, prob_thresh = 0, nbh_radius = 7)
## â—Ź AAEPS
## # A tibble: 5 x 18
##   fcdate              lead_time parameter validdate           level_type level
##   <dttm>                  <dbl> <chr>     <dttm>              <chr>      <dbl>
## 1 2018-07-10 00:00:00         0 T2m       2018-07-10 00:00:00 height         2
## 2 2018-07-10 00:00:00         3 T2m       2018-07-10 03:00:00 height         2
## 3 2018-07-10 00:00:00         6 T2m       2018-07-10 06:00:00 height         2
## 4 2018-07-10 00:00:00         9 T2m       2018-07-10 09:00:00 height         2
## 5 2018-07-10 00:00:00        12 T2m       2018-07-10 12:00:00 height         2
## # … with 12 more variables: units <chr>, fcst_cycle <chr>,
## #   AAEPS_mbr000 <geolist>, AAEPS_mbr001 <geolist>, AAEPS_mbr002 <geolist>,
## #   AAEPS_mbr003 <geolist>, AAEPS_mbr004 <geolist>, ens_mean <geolist>,
## #   ens_spread <geolist>, ens_min <geolist>, ens_max <geolist>,
## #   prob_ge_0 <geolist>

Gridded data

t2m <- ens_stats(t2m, min = TRUE, max = TRUE, prob_thresh = 0)
plot_field(t2m, "AAEPS", ens_spread, 2018071000, 12)

Gridded data

plot_field(t2m, "AAEPS", prob_ge_0, 2018071000, 12)

Gridded data

t2m <- ens_stats(t2m, prob_thresh = 0, prob_inequality = `<=`)
plot_field(t2m, "AAEPS", prob_le_0, 2018071000, 12, breaks = seq(0.2, 1, 0.2))

Gridded data

t2m <- ens_stats(t2m, prob_thresh = 0, prob_inequality = `<=`, nbh_radius = 7)
plot_field(t2m, "AAEPS", prob_le_0, 2018071000, 12, breaks = seq(0.2, 1, 0.2), palette = scico::scico(10, palette = "tokyo", direction = -1))

Observation errors

For ensembles, observation errors can be accounted for by including an error component in the ensemble forecast distribution. This is done with the jitter_forecast function.

verif        <- ens_verify(fcst, T2m)
jitter_func  <- function(x) x + rnorm(1, 0, 1)
verif_jitter <- ens_verify(jitter_fcst(fcst, jitter_func), T2m)

verif <- mutate_list(verif, `Obs error` = "No Obs errors")
verif_jitter <- mutate_list(verif_jitter, `Obs error` = "Obs errors ncluded")

verif <- bind_point_verif(verif, verif_jitter)

Observation errors

plot_point_verif(verif, spread_skill, colour_by = `Obs error`)

Observation errors

plot_point_verif(
  verif, normalized_rank_histogram, colour_by = `Obs error`, 
  rank_is_relative = TRUE, rank_hist_type = "lollipop"
)

Score cards

See Article on harphub

What’s happening next year

  • Finalize inline documentation and websites with function reference
  • harp book - Need reviewers
  • Formal tests [this is a big job!]

What’s happening next year

  • Minor improvements
    • Profiling of the more complicated functions for speed and memory use
    • Harmonize function arguments
    • Verification against analysis [read_obs() / read_analysis() function to compliment read_forecast()]
    • Simpler plotting of fields and enhanced capability (e.g. faceting)
    • Wind speed and direction from U and V
    • Verification of wind direction
    • Combined and categorical probabilities

What’s happening next year

  • More options in shiny app
    • Plot vertical profile verification
    • Different options for top level menus
    • [Graphics better optimized for web display with popups etc.: d3.js]
  • Your requests