library(harp)
library(here)Reading forecast data
In this section you will learn how to use read_grid(), read_forecast() and read_analysis() to read gridded data from files of various formats.
Before we begin, make sure that you are in your project directory for the training course and attach the packages that we are going to need.
You should have copied the data for the course into your data directory. If not, follow the instructions here.
Important: Helper packages
In order to read some file formats, specific packages are required. These can be installed as follows:
Rgrib2 remotes::install_github("harphub/Rgrib2")
ncdf4 install.packages("ncdf4")
Rfa remotes::install_github("harphub/Rfa")
read_grid()
Basic Usage
read_grid() is used to read data from a single file. You give it the file name and the parameter you want to read as well as other options to describe the data you want to read and the output.
Let’s begin by reading 2m temperature from a grib file.
read_grid(here("data/grib/exp1/mbr001/fc2022071012+000grib2_fp"), "t2m")eidb : Temperature
Time:
2022/07/10 12:00 +0
Domain summary:
161 x 211 domain
Projection summary:
proj= lcc
NE = ( 13.73206 , 57.99135 )
SW = ( 7.971006 , 53.08321 )
Data summary:
286.6064 289.4569 290.5917 290.6924 292.0204 300.3515 The output of read_grid() is an object with a class of geofield. A geofield is a 2d array with information about the grid dimensions and its co-ordinate reference system. More about this later. When printing a geofield to the screen, some information about the domain is shown as well as a summary of the data. The values are the minimum, the first quartile, the median, the mean, the third quartile and the maximum of the data in the geofield.
Parameter names
The second argument to read_grid() is the name of the parameter to read. Here we use "t2m", which is the parameter name that harp uses for 2m temperature. You can see all of the parameter names used by harp with show_param_defs()
show_param_defs()# A tibble: 61 × 2
name description
<chr> <chr>
1 accpcp12h 12-hour accumulated precipitation
2 accpcp1h 1-hour accumulated precipitation
3 accpcp24h 24-hour accumulated precipitation
4 accpcp3h 3-hour accumulated precipitation
5 accpcp6h 6-hour accumulated precipitation
6 caf Cloud area fraction at vertical levels
7 cape Convective available potential energy
8 cbase Height of cloud base
9 cc_below_7500 Cloud cover below 7500m
10 cchigh High level cloud cover
11 cclow Low level cloud cover
12 ccmed Medium level cloud cover
13 cctot Total integrated cloud cover
14 cin Convective inhibition
15 d Wind direction
16 d10m Wind direction at 10m above the ground
17 g Wind gust
18 g10m Wind gust at 10m above the ground
19 gh Geopotential height
20 gmax Maximum wind gust at 10m above the ground
21 lsm land-sea mask
22 pcp Accumulated precipitaion
23 pmsl Air pressure at mean sea level
24 pressure Atmospheric air pressure
25 psfc Surface air pressure
26 q Specific humidity of air
27 q2m Specific humidity of air at 2m above the ground
28 rh Relative humidity of air
29 rh2m Relative humidity of air at 2m above the ground
30 s Wind speed
31 s10m Wind speed at 10m above the ground
32 sea_ice Sea ice concentration
33 sfc_geo Surface geopotential
34 smax Maximum wind speed at 10m above the ground
35 snow Snow depth
36 sst Sea surface temperature
37 t Air temperature
38 t0m Skin temperature
39 t2m Air temperature at 2m above the ground
40 td Dew point temperature
41 td2m Dew point temperature at 2m above the ground
42 tmax Maximum air temperature at 2m above the ground
43 tmin Minimum air temperature at 2m above the ground
44 topo Height of topography above sea level
45 u Wind speed in u direction
46 u10m Wind speed in u direction at 10m above the ground
47 ugust Wind gust in U direction
48 ugust10m Wind gust at 10m above the ground in U direction
49 v Wind speed in v direction
50 v10m Wind speed in v direction at 10m above the ground
51 vgust Wind gust in V direction
52 vgust10m Wind gust at 10m above the ground in V direction
53 vis Visibility in air
54 w Vertical (upward) wind speed
55 wd Wind direction calculated from U and V winds
56 wd10m Wind direction at 10m above the ground calculated from U and V…
57 wg Wind gust calculated from U and V gusts
58 wg10m Wind gust at 10m above the ground calculated from U gust and V…
59 ws Wind speed calculated from U and V winds
60 ws10m Wind speed at 10m above the ground calculated from U and V win…
61 z Geopotential You can see how harp translates the parameter name for a particular file format with get_param_def()
get_param_def("t2m", "grib")$name
[1] "2t" "t"
$level_type
[1] "heightAboveGround" "surface"
$level
[1] 2You can see that with parameter = "t2m" harp will get the grib message with a shortName of t or 2t with a levelType of heightAboveGround or surface.
Other file formats
harp has built in functionality to read Grib, NetCDF, FA and vfld, vobs and obsoul files. Note that although the latter three do not contain gridded data, they can still be read by read_grid().
read_grid(
here("data/netcdf/meps_lagged/2024/02/15/07/mbr03/meps_sfc_24_20240215T07Z.nc"),
"pcp"
) : pcp kg/m^2
Time:
2024/02/16 07:00
Domain summary:
201 x 201 domain
Projection summary:
proj= lcc
NE = ( 13.82773 , 62.55298 )
SW = ( 5.582523 , 57.69936 )
Data summary:
1.470703 9.265625 13.08203 15.30311 19.65625 63.47949
Tip: vfld files need a
lead_time argument
vfld files do not include enough metadata to get the lead time from the contents of the files. Therefore you need to give read_grid() the lead time in order to fully populate the output. You wouldn’t normally read vfld files directly with read_grid(), but would use read_forecast() instead.
read_grid(
here("data/vfld/MEPS_preop/vfldMEPS_preopmbr000202402190003"),
"T2m",
lead_time = 3
) # A tibble: 2,083 × 9
SID lat lon model_elevation parameter station_data members lead_time
<int> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
1 1001 70.9 -8.67 9 T2m 273. NA 3
2 1010 69.3 16.1 4.7 T2m 275. NA 3
3 1014 69.2 17.9 46.4 T2m 272. NA 3
4 1015 69.6 17.8 0.6 T2m 274. NA 3
5 1018 69.2 16.0 153. T2m 273. NA 3
6 1021 70.0 18.7 24.4 T2m 275. NA 3
7 1022 70.0 18.7 24.4 T2m 275. NA 3
8 1023 69.1 18.5 66.2 T2m 271. NA 3
9 1025 69.7 18.9 3.9 T2m 272. NA 3
10 1026 69.7 18.9 1.3 T2m 272. NA 3
# ℹ 2,073 more rows
# ℹ 1 more variable: units <chr>Since vfld files contain point data, the output is a data frame with columns for other metadata. The data are in the station_data column.
Multiple geofields for one parameter
Files can contain multiple entries for a single parameter. For example a file could contain multiple lead times, multiple ensemble members, multiple vertical levels, or any combination of these. When many entries for a single parameter exist, a geolist is returned.
read_grid(
here("data/netcdf/arome_arctic/2024/02/19/arome_arctic_2_5km_20240219T12Z.nc"),
"t2m"
)<harp_geolist[13]>
[[1]] <numeric geofield [200 x 200] Min = 257.195 Max = 277.968 Mean = 270.880>
[[2]] <numeric geofield [200 x 200] Min = 257.907 Max = 278.076 Mean = 270.867>
[[3]] <numeric geofield [200 x 200] Min = 257.632 Max = 278.233 Mean = 270.622>
[[4]] <numeric geofield [200 x 200] Min = 256.557 Max = 278.322 Mean = 270.353>
[[5]] <numeric geofield [200 x 200] Min = 255.851 Max = 278.176 Mean = 270.146>
[[6]] <numeric geofield [200 x 200] Min = 254.634 Max = 278.187 Mean = 269.967>
[[7]] <numeric geofield [200 x 200] Min = 253.723 Max = 278.262 Mean = 269.810>
[[8]] <numeric geofield [200 x 200] Min = 253.297 Max = 278.318 Mean = 269.667>
[[9]] <numeric geofield [200 x 200] Min = 252.770 Max = 278.224 Mean = 269.533>
[[10]] <numeric geofield [200 x 200] Min = 252.329 Max = 277.918 Mean = 269.387>
# 3 more geofields
# Use `print(n = ...)` to see moreA geolist is simply a list of geofields that are on the same domain. To get more metadata about each geofield you can set data_frame = TRUE to get the output in the form of a data frame.
read_grid(
here("data/netcdf/arome_arctic/2024/02/19/arome_arctic_2_5km_20240219T12Z.nc"),
"t2m",
data_frame = TRUE
) # A tibble: 13 × 8
fcst_dttm valid_dttm lead_time gridded_data units
<dttm> <dttm> <dbl> <geolist> <chr>
1 2024-02-19 12:00:00 2024-02-19 12:00:00 0 [200 × 200] K
2 2024-02-19 12:00:00 2024-02-19 13:00:00 1 [200 × 200] K
3 2024-02-19 12:00:00 2024-02-19 14:00:00 2 [200 × 200] K
4 2024-02-19 12:00:00 2024-02-19 15:00:00 3 [200 × 200] K
5 2024-02-19 12:00:00 2024-02-19 16:00:00 4 [200 × 200] K
6 2024-02-19 12:00:00 2024-02-19 17:00:00 5 [200 × 200] K
7 2024-02-19 12:00:00 2024-02-19 18:00:00 6 [200 × 200] K
8 2024-02-19 12:00:00 2024-02-19 19:00:00 7 [200 × 200] K
9 2024-02-19 12:00:00 2024-02-19 20:00:00 8 [200 × 200] K
10 2024-02-19 12:00:00 2024-02-19 21:00:00 9 [200 × 200] K
11 2024-02-19 12:00:00 2024-02-19 22:00:00 10 [200 × 200] K
12 2024-02-19 12:00:00 2024-02-19 23:00:00 11 [200 × 200] K
13 2024-02-19 12:00:00 2024-02-20 00:00:00 12 [200 × 200] K
# ℹ 3 more variables: parameter <chr>, level_type <chr>, level <dbl>Other arguments to read_grid() allow you to specify which geofields for a parameter to get from the file. For example, to get lead times of 6 and 12 hours you would do this:
read_grid(
here("data/netcdf/arome_arctic/2024/02/19/arome_arctic_2_5km_20240219T12Z.nc"),
"t2m",
lead_time = c(6, 12),
data_frame = TRUE
) # A tibble: 2 × 8
fcst_dttm valid_dttm lead_time gridded_data units parameter
<dttm> <dttm> <dbl> <geolist> <chr> <chr>
1 2024-02-19 12:00:00 2024-02-19 18:00:00 6 [200 × 200] K t2m
2 2024-02-19 12:00:00 2024-02-20 00:00:00 12 [200 × 200] K t2m
# ℹ 2 more variables: level_type <chr>, level <dbl>File format options
Grib
Sometimes there isn’t sufficient information in the file, or the defaults are incorrect. Take, for example, a grib2 file that uses non standard parameter numbers for total precipitation.
read_grid(
here("data/grib/exp1/mbr001/fc2022071012+006grib2_fp"),
"pcp"
)Warning: Parameter "pcp" (shortName: "tp", typeOfLevel: "heightAboveGround" /
"surface" for level(s) 0) not found in grib file.Error: None of the requested data could be read from grib file: /home/andrewts/R-projects/harp_training_2024/data/grib/exp1/mbr001/fc2022071012+006grib2_fpWe see that it is looking for a parameter with shortName "tp", but cannot find it. Yet we know the file contains total precipitation. If you have access to the grib tables used to encode this grib file you can look up the correct information to get the data. In this case total precipitation uses the grib2 parameterNumber 8.
Options for the file format are passed through the file_format_opts argument, and those options can be generated with helper functions for the file format. For grib files those options can be generated by grib_opts().
grib_opts()$meta
[1] TRUE
$multi
[1] FALSE
$param_find
NULL
$level_find
NULL
$step_find
NULLHere, one of the arguments is param_find and we use that in conjunction with the use_grib_*() function use_grib_parameterNumber(). The param_find argument takes a named list, with the name the same as the parameter so that it knows which parameter to which to apply those options.
read_grid(
here("data/grib/exp1/mbr001/fc2022071012+006grib2_fp"),
"pcp",
file_format_opts = grib_opts(
param_find = list(pcp = use_grib_parameterNumber(8))
)
)eidb : Total precipitation
Time:
2022/07/10 12:00 +0
Domain summary:
161 x 211 domain
Projection summary:
proj= lcc
NE = ( 13.73206 , 57.99135 )
SW = ( 7.971006 , 53.08321 )
Data summary:
0 0 0 0.0005843104 0 0.07226562 NetCDF
Reading NetCDF files with harp used to require the user to always specify the names of the dimensions and the projection variable in the file. This has improved in recent versions, and as shown above you can often read from NetCDF files without providing extra options. However, in some cases you still have to pass some information about the content of the file. This happens when the x and y dimensions are not named “x” and “y”, and/or the projection is not a lambert projection. Additionally if the data are stored in reverse order this must be specified. An example of this is would be forecasts from a global model like ECMWF’s IFS.
In the below example, we tell the function that the x and y dimensions are “longitude” and “latitude”, the projection variable is “projection_regular_ll” and that the data are stored in reverse order in the y dimension (north - south as opposed to south - north).
read_grid(
here("data/netcdf/ifsens/ifsens_20240219T000000Z.nc"),
"t2m",
lead_time = 6,
members = 0,
file_format_opts = netcdf_opts(
x_dim = "longitude",
y_dim = "latitude",
y_rev = TRUE,
proj4_var = "projection_regular_ll"
)
) : t2m K
Time:
2024/02/19 00:00
Domain summary:
81 x 80 domain
Projection summary:
proj= longlat
NE = ( -4 , 57.9 )
SW = ( -12 , 50 )
Data summary:
276.9579 281.7704 283.1766 282.7288 283.8329 285.3329 read_forecast()
read_grid() is used to read data from a single file. In many cases, you will need to read data from multiple files. This is where read_forecast() comes in.
read_forecast() takes at a minimum the date-times you want to read, the name of the forecast model to give to the data, the parameter you want to read, the path of the files and a file template. Taking the first file we read in this tutorial, here("data/grib/exp1/mbr001/fc2022071012+000grib2_fp"), we can read the same data using read_forecast() with the following expression:
read_forecast(
dttm = 2022071012,
fcst_model = "exp1",
parameter = "t2m",
lead_time = 0,
members = 1,
file_path = here("data/grib"),
file_template = "{fcst_model}/mbr{MBR3}/fc{YYYY}{MM}{DD}{HH}+{LDT3}grib2_fp",
return_data = TRUE
)::ensemble gridded forecast:: # A tibble: 1 × 11
fcst_dttm valid_dttm lead_time parameter exp1_mbr001
<dttm> <dttm> <dbl> <chr> <geolist>
1 2022-07-10 12:00:00 2022-07-10 12:00:00 0 t2m [161 × 211]
# ℹ 6 more variables: fcst_model <chr>, step_range <chr>, level_type <chr>,
# level <int>, units <chr>, fcst_cycle <chr>
Tip: You need to tell
read_forecast() to return data
One of the primary functions of read_forecast() is to process large volumes of data to interpolate to points and write these interpolated data to new files. Such an operation would often lead to running out of memory if the data were returned to the global environment so the default behaviour is to not return any data. Therefore, if you want the data returned to your global environment, you must set return_data = TRUE. Note that this default behaviour is under review and may be changed in future versions.
File Name Templates
In the above example, the file names are generated by replacing everything that is inside braces with dynamic data. That is to say, values that can change depending on the date-time, the lead time, the ensemble member and the name of the forecast model. We refer to the embraced values as substitutions, and the available substitutions are listed below:
| Substitution | Description |
|---|---|
{YYYY} |
4 digit year |
{MM} |
2 digit month with leading zeros |
{M} |
Single digit month |
{DD} |
2 digit day with leading zeros |
{D} |
Single digit day |
{HH} |
2 digit hour with leading zeros |
{H} |
Single digit hour |
{mm} |
2 digit minute with leading zeros |
{m} |
Single digit minute |
{LDTx} |
Lead time |
{MBRx} |
Ensemble member |
In the above table, the substitutions {LDTx} and {MBRx} have an optional number of digits to use, where smaller values use leading zeros. The x should be replaced by the number of digits that the file name uses. Leaving off a value for x means that no leading zeros are used.
Getting the file template correct can often be quite a trying process, so harp includes some built in templates that can be seen with show_file_templates().
show_file_templates()# A tibble: 34 × 2
template_name template
<chr> <chr>
1 arome_arctic_extracted /lustre/storeB/immutable/archive/projects/metproducti…
2 arome_arctic_full /lustre/storeB/immutable/archive/projects/metproducti…
3 arome_arctic_sfx /lustre/storeB/immutable/archive/projects/metproducti…
4 fctable {file_path}/{fcst_model}/{YYYY}/{MM}/FCTABLE_{paramet…
5 fctable_det {file_path}/{det_model}/{YYYY}/{MM}/FCTABLE_{paramete…
6 fctable_eps {file_path}/{eps_model}/{YYYY}/{MM}/FCTABLE_{paramete…
7 fctable_eps_all_cycles {file_path}/{eps_model}/{YYYY}/{MM}/FCTABLE_{paramete…
8 fctable_eps_all_leads {file_path}/{eps_model}/{YYYY}/{MM}/FCTABLE_{paramete…
9 glameps_grib {file_path}/{eps_model}/{sub_model}/{YYYY}/{MM}/{DD}/…
10 harmoneps_grib {file_path}/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR3}/fc{YYYY}{…
11 harmoneps_grib_fp {file_path}/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR3}/fc{YYYY}{…
12 harmoneps_grib_sfx {file_path}/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR3}/fc{YYYY}{…
13 harmonie_grib {file_path}/{YYYY}/{MM}/{DD}/{HH}/fc{YYYY}{MM}{DD}{HH…
14 harmonie_grib_fp {file_path}/{YYYY}/{MM}/{DD}/{HH}/fc{YYYY}{MM}{DD}{HH…
15 harmonie_grib_sfx {file_path}/{YYYY}/{MM}/{DD}/{HH}/fc{YYYY}{MM}{DD}{HH…
16 meps_cntrl_extracted /lustre/storeB/immutable/archive/projects/metproducti…
17 meps_cntrl_sfx /lustre/storeB/immutable/archive/projects/metproducti…
18 meps_det /lustre/storeB/immutable/archive/projects/metproducti…
19 meps_extracted /lustre/storeB/immutable/archive/projects/metproducti…
20 meps_full /lustre/storeB/immutable/archive/projects/metproducti…
21 meps_lagged_6h_subset /lustre/storeB/immutable/archive/projects/metproducti…
22 meps_sfx /lustre/storeB/immutable/archive/projects/metproducti…
23 meps_subset /lustre/storeB/immutable/archive/projects/metproducti…
24 obsoul {file_path}/obsoul_1_xxxxxy_{country}_{YYYY}{MM}{DD}{…
25 obstable {file_path}/OBSTABLE_{YYYY}.sqlite
26 vfld {file_path}/{fcst_model}/vfld{fcst_model}{YYYY}{MM}{D…
27 vfld_det {file_path}/{det_model}/vfld{det_model}{YYYY}{MM}{DD}…
28 vfld_det_noexp {file_path}/{det_model}/vfld{YYYY}{MM}{DD}{HH}{LDT2}
29 vfld_eps {file_path}/{sub_model}/vfld{sub_model}mbr{MBR3}{YYYY…
30 vfld_eps_noexp {file_path}/{sub_model}/vfldmbr{MBR3}{YYYY}{MM}{DD}{H…
31 vfld_multimodel {file_path}/{sub_model}/vfld{sub_model}mbr{MBR3}{YYYY…
32 vfld_multimodel_noexp {file_path}/{sub_model}/vfldmbr{MBR3}{YYYY}{MM}{DD}{H…
33 vfld_noexp {file_path}/{fcst_model}/vfldmbr{MBR3}{YYYY}{MM}{DD}{…
34 vobs {file_path}/vobs{YYYY}{MM}{DD}{HH} Often the templates are a bit long to be seen on screen so a single template can be shown by selecting the number for that template.
show_file_templates(29)
template_name:
vfld_eps
template:
{file_path}/{sub_model}/vfld{sub_model}mbr{MBR3}{YYYY}{MM}{DD}{HH}{LDT2} This means we can, for example, read an ensemble vfld file with the "vfld_eps" template. Here we are setting parameter = NULLto read all parameters from the vfld file (this is only an option for vfld files).
read_forecast(
dttm = 2024021900,
fcst_model = "MEPS_preop",
parameter = NULL,
lead_time = 3,
members = 0,
file_path = here("data/vfld"),
file_template = "vfld_eps",
return_data = TRUE
)::ensemble point forecast:: # A tibble: 80,746 × 12
fcst_dttm lead_time parameter SID MEPS_preop_mbr000 units
<dttm> <dbl> <chr> <int> <dbl> <chr>
1 2024-02-19 00:00:00 3 CCtot 1001 8 oktas
2 2024-02-19 00:00:00 3 CCtot 1010 8 oktas
3 2024-02-19 00:00:00 3 CCtot 1014 8 oktas
4 2024-02-19 00:00:00 3 CCtot 1015 8 oktas
5 2024-02-19 00:00:00 3 CCtot 1018 8 oktas
6 2024-02-19 00:00:00 3 CCtot 1021 8 oktas
7 2024-02-19 00:00:00 3 CCtot 1022 8 oktas
8 2024-02-19 00:00:00 3 CCtot 1023 7 oktas
9 2024-02-19 00:00:00 3 CCtot 1025 8 oktas
10 2024-02-19 00:00:00 3 CCtot 1026 8 oktas
# ℹ 80,736 more rows
# ℹ 6 more variables: fcst_model <chr>, lat <dbl>, lon <dbl>, p <dbl>,
# valid_dttm <dttm>, fcst_cycle <chr>Data classes
Unlike read_grid(), read_forecast() always returns a data frame. These data frames have a harp_df class and then an attempt is made to assign a subclass depending on the data. In the above examples you will see one output shown as
::ensemble gridded forecast::
and the other as
::ensemble point forecast::
This means that the function has recognised that these are ensemble gridded forecasts and ensemble point forecasts and assigned the appropriate sub classes. Since in both of these cases only one ensemble member has been read in, they can be converted to deterministic forecasts using as_det()
read_forecast(
dttm = 2024021900,
fcst_model = "MEPS_preop",
parameter = NULL,
lead_time = 3,
members = 0,
file_path = here("data/vfld"),
file_template = "vfld_eps",
return_data = TRUE
) |>
as_det()::deterministic point forecast:: # A tibble: 80,746 × 12
fcst_dttm lead_time parameter SID fcst units fcst_model lat
<dttm> <dbl> <chr> <int> <dbl> <chr> <chr> <dbl>
1 2024-02-19 00:00:00 3 CCtot 1001 8 oktas MEPS_preop 70.9
2 2024-02-19 00:00:00 3 CCtot 1010 8 oktas MEPS_preop 69.3
3 2024-02-19 00:00:00 3 CCtot 1014 8 oktas MEPS_preop 69.2
4 2024-02-19 00:00:00 3 CCtot 1015 8 oktas MEPS_preop 69.6
5 2024-02-19 00:00:00 3 CCtot 1018 8 oktas MEPS_preop 69.2
6 2024-02-19 00:00:00 3 CCtot 1021 8 oktas MEPS_preop 70.0
7 2024-02-19 00:00:00 3 CCtot 1022 8 oktas MEPS_preop 70.0
8 2024-02-19 00:00:00 3 CCtot 1023 7 oktas MEPS_preop 69.1
9 2024-02-19 00:00:00 3 CCtot 1025 8 oktas MEPS_preop 69.7
10 2024-02-19 00:00:00 3 CCtot 1026 8 oktas MEPS_preop 69.7
# ℹ 80,736 more rows
# ℹ 4 more variables: lon <dbl>, p <dbl>, valid_dttm <dttm>, fcst_cycle <chr>The output is now labelled as ::deterministic point forecast:: and the data column is now simply fcst.
Tip:
|> The pipe operator
In the above we also use R’s pipe operator |>. The pipe operator takes the result of the function that comes before it and passes it to the function that comes after it as the first argument. It can be thought of as “and then” or “and send to”.
While the harp_df classes may not be something you need to know much about, many functions in harp rely on these classes in order to know what to do with the data.
Any data frame can be converted to a harp_df data frame as long as it has a valid_dttm column. That is to say a column that contains the valid date time for each row of data. Deterministic data frames are initially recognised by having a column name that ends with “_det” and ensemble columns are recognised by column names that end with “_mbrXXX”, where XXX is a 3 digit member number with leading zeros.
To demonstrate we can use some columns from harp’s built in test data to construct a data frame and give it a harp class using as_harp_df().
point_df <- data.frame(
fcst_dttm = det_point_df$fcst_dttm[1:5],
valid_dttm = det_point_df$valid_dttm[1:5],
lead_time = det_point_df$lead_time[1:5],
SID = det_point_df$SID[1:5],
point_det = det_point_df$fcst[1:5]
)
class(point_df)[1] "data.frame"point_df <- as_harp_df(point_df)
class(point_df)[1] "harp_det_point_df" "harp_point_df" "harp_df"
[4] "tbl_df" "tbl" "data.frame" The harp_df class and subclasses can be removed with deharp(). This can be important as a small number functions that take simple data frames will not work with harp_df data frames.
class(deharp(point_df))[1] "tbl_df" "tbl" "data.frame"Geographic transformations
read_forecast() and read_grid() have the ability to perform geographic transformations on the data at read time. These transformations include interpolating to point locations, regridding to another grid, taking a subset of a grid or pulling out a line section through a grid to create cross sections. Here we will concentrate on the interpolation to point locations, transformation = "interpolate"
Interpolation to geographic point locations
When selecting a transformation, there is an accompanying transformation_opts argument to pass the options for that transformation. Each transformation has a function that generates those options - in the case of the “interpolate” transformation, that is interpolate_opts(), which has a number of default options.
interpolate_opts()No stations specified. Using default stations: 'station_list'$stations
# A tibble: 13,417 × 5
SID lat lon elev name
<int> <dbl> <dbl> <dbl> <chr>
1 1001 70.9 -8.67 9.4 JAN MAYEN
2 1002 80.1 16.2 8 VERLEGENHUKEN
3 1003 77 15.5 11.1 HORNSUND
4 1004 78.9 11.9 8 NY-ALESUND II
5 1006 78.3 22.8 14 EDGEOYA
6 1007 78.9 11.9 7.7 NY-ALESUND
7 1008 78.2 15.5 26.8 SVALBARD AP
8 1009 80.7 25.0 5 KARL XII OYA
9 1010 69.3 16.1 13.1 ANDOYA
10 1011 80.1 31.5 10 KVITOYA
# ℹ 13,407 more rows
$method
[1] "nearest"
$correct_t2m
[1] TRUE
$keep_model_t2m
[1] FALSE
$lapse_rate
[1] 0.0065
$clim_file
NULL
$clim_file_format
NULL
$clim_file_opts
NULL
$clim_param
[1] "sfc_geo"
$use_mask
[1] FALSE
$weights
NULL
$keep_raw_data
[1] FALSEThe most important information for the “interpolate” transformation is the locations to which to interpolate the data. This is provided in the stations argument, and by default harp’s built in list of meteorological stations, station_list is used. If not using the default, stations needs to be a data frame with the following columns:
| Column name | Description |
|---|---|
| SID | A unique identifier for the location |
| lat | The latitude of the location in decimal degrees |
| lon | The longitude of the location in decimal degrees |
| elev* | The elevation of the location in meters |
| name* | A name for the location |
Column names marked with * are optional.
For a first example we will read in exp1 10m wind speed (“S10m”) for lead times 0 - 3 hours, and members 1 -3, and interpolate to the the default stations (in this case we do not need to set anything for transformation_opts as we will use the defaults).
read_forecast(
dttm = 2022071012,
fcst_model = "exp1",
parameter = "s10m",
lead_time = seq(0, 3),
members = seq(1, 3),
file_path = here("data/grib"),
file_template = "{fcst_model}/mbr{MBR3}/fc{YYYY}{MM}{DD}{HH}+{LDT3}grib2_fp",
transformation = "interpolate",
return_data = TRUE
)Warning: 'transformation_opts' not set for transformation = 'interpolate'.
Using default interpolate_opts()Warning: Parameter "sfc_geo" (shortName: "z", typeOfLevel: "heightAboveGround"
/ "surface" for level(s) 0) not found in grib file.Error : None of the requested data could be read from grib file: /home/andrewts/R-projects/harp_training_2024/data/grib/exp1/mbr001/fc2022071012+000grib2_fp::ensemble point forecast:: # A tibble: 628 × 16
fcst_dttm valid_dttm lead_time SID exp1_mbr001
<dttm> <dttm> <dbl> <int> <dbl>
1 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2512 5.29
2 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2513 4.98
3 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2516 4.96
4 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2517 4.28
5 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2518 5.92
6 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2526 5.84
7 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2527 5.85
8 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2531 3.90
9 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2536 5.45
10 2022-07-10 12:00:00 2022-07-10 12:00:00 0 2539 5.42
# ℹ 618 more rows
# ℹ 11 more variables: exp1_mbr002 <dbl>, exp1_mbr003 <dbl>, units <chr>,
# fcst_model <chr>, parameter <chr>, step_range <chr>, level_type <chr>,
# level <int>, lat <dbl>, lon <dbl>, fcst_cycle <chr>
Note: Warnings and errors
Interpolation weights need to be computed from one of the fields in the file. By default it tries to get these from the surface geopotential. If surface geopotential is not found in the file a warning and error are thrown, but it simply continues to compute the interpolation weights from the first parameter to be read. It’s nothing to worry about!
It is also possible to interpolate to any point location (as long as it exists inside the domain - in this case over Denmark) with a data frame to be sent to the stations argument. In this case, we will also use bilinear interpolation (the default is nearest neighbour).
my_stations <- data.frame(
SID = c("CPH", "AAR"),
lon = c(12.64, 10.62),
lat = c(55.61, 56.31)
)
Tip: Station IDs (SID) do not have to be numbers
read_forecast(
dttm = 2022071012,
fcst_model = "exp1",
parameter = "s10m",
lead_time = seq(0, 3),
members = seq(1, 3),
file_path = here("data/grib"),
file_template = "{fcst_model}/mbr{MBR3}/fc{YYYY}{MM}{DD}{HH}+{LDT3}grib2_fp",
transformation = "interpolate",
transformation_opts = interpolate_opts(
stations = my_stations,
method = "bilinear"
),
return_data = TRUE
)Warning: No 'elev' column found in stations, and correct_t2m = TRUE. Setting
correct_t2m = FALSEWarning: Parameter "sfc_geo" (shortName: "z", typeOfLevel: "heightAboveGround"
/ "surface" for level(s) 0) not found in grib file.Error : None of the requested data could be read from grib file: /home/andrewts/R-projects/harp_training_2024/data/grib/exp1/mbr001/fc2022071012+000grib2_fp::ensemble point forecast:: # A tibble: 8 × 16
fcst_dttm valid_dttm lead_time SID exp1_mbr001
<dttm> <dttm> <dbl> <chr> <dbl>
1 2022-07-10 12:00:00 2022-07-10 12:00:00 0 CPH 6.23
2 2022-07-10 12:00:00 2022-07-10 12:00:00 0 AAR 8.32
3 2022-07-10 12:00:00 2022-07-10 13:00:00 1 CPH 6.65
4 2022-07-10 12:00:00 2022-07-10 13:00:00 1 AAR 8.36
5 2022-07-10 12:00:00 2022-07-10 14:00:00 2 CPH 6.73
6 2022-07-10 12:00:00 2022-07-10 14:00:00 2 AAR 8.14
7 2022-07-10 12:00:00 2022-07-10 15:00:00 3 CPH 6.93
8 2022-07-10 12:00:00 2022-07-10 15:00:00 3 AAR 7.94
# ℹ 11 more variables: exp1_mbr002 <dbl>, exp1_mbr003 <dbl>, units <chr>,
# fcst_model <chr>, parameter <chr>, step_range <chr>, level_type <chr>,
# level <int>, lat <dbl>, lon <dbl>, fcst_cycle <chr>Once the data are interpolated they can be output to files in SQLite format. SQLite files allow the data to be filtered and read much more quickly for further use, for example in verification. This is done with the argument output_file_opts and the options can be set with fctable_opts(), with the most important of those options being the path to which to write the files.
read_forecast(
dttm = 2022071012,
fcst_model = "exp1",
parameter = "s10m",
lead_time = seq(0, 3),
members = seq(1, 3),
file_path = here("data/grib"),
file_template = "{fcst_model}/mbr{MBR3}/fc{YYYY}{MM}{DD}{HH}+{LDT3}grib2_fp",
transformation = "interpolate",
transformation_opts = interpolate_opts(
stations = my_stations,
method = "bilinear"
),
output_file_opts = fctable_opts(path = here("data/FCTABLE"))
)Warning: No 'elev' column found in stations, and correct_t2m = TRUE. Setting
correct_t2m = FALSEWarning: Parameter "sfc_geo" (shortName: "z", typeOfLevel: "heightAboveGround"
/ "surface" for level(s) 0) not found in grib file.Error : None of the requested data could be read from grib file: /home/andrewts/R-projects/harp_training_2024/data/grib/exp1/mbr001/fc2022071012+000grib2_fpWhen interpolating 2m temperature to geographic point locations an attempt is made to correct the temperature at the topographic height in the model to the actual topographic heights of the geographic point locations. This is done using a simple lapse rate conversion of 0.0065 K.m-1, although this can be set in interpolate_opts(). In order for this correction to happen, the stations data frame needs to have an “elev” column that has the contains the elevation of the station in meters. Furthermore, there needs to be information about the model elevation. By default the model elevation is inferred from the surface geopotential (“sfc_geo”), but it can be set to topographic elevation by setting the clim_param argument to interpolate_opts() to “topo”. This information about model elevation can be in the same file as the forecast, or in a separate file, the path to which can be set in the clim_file argument to interpolate_opts().
In the following example, the surface geopotential is in the same file as the forecast, but it requires different options to read the NetCDF variable (there is no time dimension), which we can use modify_opts() to modify the existing options. Since we need to set clim_file_opts, as things stand this means we also have to populate the clim_file argument. Since we are going to use the file format options repeatedly, it makes sense to write them to a variable.
ifsens_fmt_opts <- netcdf_opts(
x_dim = "longitude",
y_dim = "latitude",
y_rev = TRUE,
proj4_var = "projection_regular_ll"
)First we will read the height corrected 2m temperature
t2m <- read_forecast(
dttm = 2024021900,
fcst_model = "ifsens",
parameter = "t2m",
lead_time = 0,
members = 0,
file_path = here("data/netcdf"),
file_template = "{fcst_model}/{fcst_model}_{YYYY}{MM}{DD}T{HH}{mm}00Z.nc",
file_format_opts = ifsens_fmt_opts,
transformation = "interpolate",
transformation_opts = interpolate_opts(
clim_file = here("data/netcdf/ifsens/ifsens_20240219T000000Z.nc"),
clim_file_opts = modify_opts(ifsens_fmt_opts, time_var = NA)
),
return_data = TRUE
)And now the uncorrected, by setting correct_t2m = FALSE in interpolate_opts()
t2m_uncorrected <- read_forecast(
dttm = 2024021900,
fcst_model = "ifsens",
parameter = "t2m",
lead_time = 0,
members = 0,
file_path = here("data/netcdf"),
file_template = "{fcst_model}/{fcst_model}_{YYYY}{MM}{DD}T{HH}{mm}00Z.nc",
file_format_opts = ifsens_fmt_opts,
transformation = "interpolate",
transformation_opts = interpolate_opts(correct_t2m = FALSE),
return_data = TRUE
)Error in nc_id[["var"]][[nc_param]][["dim"]][[idx]] :
attempt to select less than one element in get1indexAnd we can then compare the impact of the temperature correction.
t2m$ifsens_mbr000 - t2m_uncorrected$ifsens_mbr000 [1] 0.1458409286 0.0414019990 -0.3619363903 1.3419084521 0.1547063988
[6] 0.0507633395 0.4325378796 -3.9913955694 -0.1321044335 -4.7310049222
[11] 1.4207357056 1.5684432188 1.9425989706 1.4803179003 0.3218133717
[16] -0.0740768766 0.5275350526 0.4473404526 1.9519767322 0.5423946145
[21] 0.5653087660 0.3477330038 -0.1038059902 0.9476696306 0.5350001595
[26] -0.5841755696 -0.5483523687 0.2738574474 -0.2281737633 0.1697904205
[31] -1.4340309741 -0.4900098497 0.6684409930 -0.2920713961 -0.1685245416
[36] 0.0447920530 -0.0175268444 -0.3850984770 0.6235608309 -0.0081554177
[41] -0.0016554177 0.0040337500 0.0945613069 -0.1910567042 -0.3376315915
[46] -0.1266898831 -0.0464718076 0.2549739661 0.7939446467 0.1764446467
[51] 0.2554832418 0.5228946145 -0.0522917743 0.2616012528 0.3904855446
[56] -0.3578620560 0.0427422151 0.1355783016 0.0090695662 -0.3385631367
[61] 0.1972646904 0.4642036200 0.0509845927 0.1260929935 0.6516448525
[66] -0.0276085632 -0.0617536552 0.4499873715 -0.6436678067 -0.1589722066
[71] 0.2974116092 -0.0125245416 -0.1436168226 -0.0905580551 -0.0002582609
[76] 0.2583546686 -0.4707912983 -0.0034010500 0.1112770795 -0.8161154717
[81] -0.0619781631 0.1267503825 0.1419074795 -0.1062954280Reading SQLite files
Once SQLite files have been written by read_forecast(), they can then be read using read_point_forecast(), which works in much the same way as read_forecast(). The main difference is that you have to specify whether you are reading ensemble (“eps”) or deterministic (“det”) forecasts via the fcst_type argument.
read_point_forecast(
dttm = 2022071012,
fcst_model = "exp1",
fcst_type = "eps",
parameter = "s10m",
file_path = here("data/FCTABLE")
)::ensemble point forecast:: # A tibble: 4 × 12
fcst_dttm valid_dttm lead_time SID exp1_mbr001
<dttm> <dttm> <int> <chr> <dbl>
1 2022-07-10 12:00:00 2022-07-10 12:00:00 0 AAR 8.32
2 2022-07-10 12:00:00 2022-07-10 12:00:00 0 CPH 6.23
3 2022-07-10 12:00:00 2022-07-10 15:00:00 3 AAR 7.94
4 2022-07-10 12:00:00 2022-07-10 15:00:00 3 CPH 6.93
# ℹ 7 more variables: exp1_mbr002 <dbl>, exp1_mbr003 <dbl>, units <chr>,
# fcst_model <chr>, parameter <chr>, z <int>, fcst_cycle <chr>read_point-forecast() will be explained further in the Point Verification Workflow tutorial
Lagged ensembles
Many ensembles are generated with time lags with the goal of maximising the number of ensemble members while making most efficient use of computing resources. This means that output files for different members are produced for different forecast times. This creates a challenge in reading in the data using read_forecast()’s templating system. However, the lags argument allows you to specify which members have output files for which forecast times. This is done by having a vector that is the same length as the vector for members with the corresponding lag for each member.
MEPS, the ensemble run by MetCoOp (a collaboration between MET Norway, SMHI and FMI) is one such example. It produces 5 ensemble members every hour, resulting in 15 independent members every 3 hours. If we take the 12:00 UTC 15 member ensemble as an example, the members are allocated as follows:
| Time (UTC) | Members |
|---|---|
| 12:00 | 0, 1, 2, 9, 12 |
| 11:00 | 5, 6, 8, 11, 14 |
| 10:00 | 3, 4, 7, 10, 13 |
This means that the members at 11:00 UTC have a lag of 1 hour and those at 10:00 UTC have a lag of 2 hours.
Rather than have to write the members and lags out as a vector, we can make a nested list of the lagging information and a function to generate those vectors.
meps_mbrs <- list(
list(lag = 0, members = c(0, 1, 2, 9, 12)),
list(lag = 1, members = c(5, 6, 8, 11, 14)),
list(lag = 2, members = c(3, 4, 7, 10, 13))
)
get_mbrs <- function(x, mbr = "members") {
unlist(lapply(x, \(d) d[[mbr]]))
}
get_lags <- function(x, lag = "lag", mbr = "members") {
unlist(lapply(x, \(d) rep(d[[lag]], length(d[[mbr]]))))
}
get_mbrs(meps_mbrs) [1] 0 1 2 9 12 5 6 8 11 14 3 4 7 10 13get_lags(meps_mbrs) [1] 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2And now we can run read_forecast() for those lags and members
read_forecast(
dttm = 2024021512,
fcst_model = "meps",
parameter = "pcp",
lead_time = 0,
members = get_mbrs(meps_mbrs),
lags = get_lags(meps_mbrs),
file_path = here("data/netcdf"),
file_template = "{fcst_model}_lagged/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR2}/{fcst_model}_sfc_{LDT2}_{YYYY}{MM}{DD}T{HH}Z.nc",
file_format_opts = netcdf_opts(ref_time_var = "forecast_reference_time"),
return_data = TRUE
)::ensemble gridded forecast:: # A tibble: 1 × 24
fcst_dttm valid_dttm lead_time parameter meps_mbr000
<dttm> <dttm> <dbl> <chr> <geolist>
1 2024-02-15 12:00:00 2024-02-15 12:00:00 0 pcp [201 × 201]
# ℹ 19 more variables: meps_mbr001 <geolist>, meps_mbr002 <geolist>,
# meps_mbr009 <geolist>, meps_mbr012 <geolist>, meps_mbr005_lag1h <geolist>,
# meps_mbr006_lag1h <geolist>, meps_mbr008_lag1h <geolist>,
# meps_mbr011_lag1h <geolist>, meps_mbr014_lag1h <geolist>,
# meps_mbr003_lag2h <geolist>, meps_mbr004_lag2h <geolist>,
# meps_mbr007_lag2h <geolist>, meps_mbr010_lag2h <geolist>,
# meps_mbr013_lag2h <geolist>, units <chr>, fcst_model <chr>, …
Errors and Warnings
You may see warnings along the likes of Ensemble members were requested for 'pcp' but there is no member information. This occurs because the NetCDF files do not have an ensemble_member dimension since each file only contains one ensemble member.
Since these 15 members are produced by MEPS every three hours, then to read more than one ensemble the value for dttm should have a 3 hour time step, This can be achieved by using seq_dttm() to generate a sequence of date-times as will be seen in the next example.
Tip: Writing interpolated lagged forecasts to SQLite
When writing lagged forecasts to SQLite files the members are not all collected together to form the full ensemble as was returned in the above example. Rather the members for each lag for the full ensemble are collected together in SQLite files for each date-time that contributes to the ensemble. This allows for more flexibility in constructing lagged ensembles from SQLite files and has implications for how you specify lags for read_point_forecast() as will be seen below.
Here we will get the lagged ensemble for 2 forecasts 3 hours apart, interpolate the data to points and write the results to SQLite files.
read_forecast(
dttm = seq_dttm(2024021509, 2024021512, "3h"),
fcst_model = "meps",
parameter = "pcp",
lead_time = seq(0, 6),
members = get_mbrs(meps_mbrs),
lags = get_lags(meps_mbrs),
file_path = here("data/netcdf"),
file_template = "{fcst_model}_lagged/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR2}/{fcst_model}_sfc_{LDT2}_{YYYY}{MM}{DD}T{HH}Z.nc",
file_format_opts = netcdf_opts(ref_time_var = "forecast_reference_time"),
transformation = "interpolate",
output_file_opts = fctable_opts(path = here("data/FCTABLE"))
)Error : None of the requested parameters found in file: /home/andrewts/R-projects/harp_training_2024/data/netcdf/meps_lagged/2024/02/15/09/mbr00/meps_sfc_00_20240215T09Z.ncWhen reading the data back with read_point_forecast(), we only need to provide the lags (here we have to specify that they are in hours), since all of the members are collected together for each lag.
read_point_forecast(
dttm = seq_dttm(2024021509, 2024021512, "3h"),
fcst_model = "meps",
fcst_type = "eps",
parameter = "pcp",
lags = paste0(seq(0, 2), "h"),
file_path = here("data/FCTABLE")
)::ensemble point forecast:: # A tibble: 1,038 × 23
fcst_dttm valid_dttm lead_time SID meps_mbr000
<dttm> <dttm> <dbl> <int> <dbl>
1 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1203 0
2 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1207 0
3 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1209 0
4 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1218 0
5 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1230 0
6 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1233 0
7 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1238 0
8 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1239 0
9 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1250 0
10 2024-02-15 09:00:00 2024-02-15 09:00:00 0 1265 0
# ℹ 1,028 more rows
# ℹ 18 more variables: meps_mbr001 <dbl>, meps_mbr002 <dbl>, meps_mbr009 <dbl>,
# meps_mbr012 <dbl>, meps_mbr005 <dbl>, meps_mbr006 <dbl>, meps_mbr008 <dbl>,
# meps_mbr011 <dbl>, meps_mbr014 <dbl>, meps_mbr003 <dbl>, meps_mbr004 <dbl>,
# meps_mbr007 <dbl>, meps_mbr010 <dbl>, meps_mbr013 <dbl>, units <chr>,
# fcst_model <chr>, parameter <chr>, fcst_cycle <chr>Here the member column names do not include how much that member is lagged.
Since the data are precipitation, we can return the accumulated precipitation for given accumulation periods. This is done by prefixing the parameter name with “Acc” and following it with the accumulation period. So, 3 hour accumulated precipitation would be “Accpcp3h”. We can read any selection of stations by specifying their IDs (SID) in the stations argument.
read_point_forecast(
dttm = seq_dttm(2024021509, 2024021512, "3h"),
fcst_model = "meps",
fcst_type = "eps",
parameter = "Accpcp3h",
lags = paste0(seq(0, 2), "h"),
stations = c(1425, 1439),
file_path = here("data/FCTABLE")
)::ensemble point forecast:: # A tibble: 8 × 23
fcst_dttm valid_dttm lead_time SID meps_mbr000
<dttm> <dttm> <dbl> <int> <dbl>
1 2024-02-15 09:00:00 2024-02-15 12:00:00 3 1425 2.77
2 2024-02-15 09:00:00 2024-02-15 12:00:00 3 1439 3.48
3 2024-02-15 09:00:00 2024-02-15 15:00:00 6 1425 5.74
4 2024-02-15 09:00:00 2024-02-15 15:00:00 6 1439 5.28
5 2024-02-15 12:00:00 2024-02-15 15:00:00 3 1425 6.02
6 2024-02-15 12:00:00 2024-02-15 15:00:00 3 1439 5.47
7 2024-02-15 12:00:00 2024-02-15 18:00:00 6 1425 9.06
8 2024-02-15 12:00:00 2024-02-15 18:00:00 6 1439 5.98
# ℹ 18 more variables: meps_mbr001 <dbl>, meps_mbr002 <dbl>, meps_mbr009 <dbl>,
# meps_mbr012 <dbl>, meps_mbr005 <dbl>, meps_mbr006 <dbl>, meps_mbr008 <dbl>,
# meps_mbr011 <dbl>, meps_mbr014 <dbl>, meps_mbr003 <dbl>, meps_mbr004 <dbl>,
# meps_mbr007 <dbl>, meps_mbr010 <dbl>, meps_mbr013 <dbl>, units <chr>,
# fcst_model <chr>, parameter <chr>, fcst_cycle <chr>In the next tutorial we will going through the workflow for doing point verification.