harpSpatial Test: harpSpatialTeam version

This notebook is a test of new developments in the harpSpatial package, currently in the harpSpatialTeam repository on Github. Here spatial verification is run on the same data that used for the current version of harpSpatial that is available in the harphub repository on Github.

Begin by

library(harp)
library(here)
obs_dir      <- here("..", "data", "met_analysis")
obs_template <- "met_analysis_1_0km_nordic_{YYYY}{MM}{DD}T{HH}Z.nc"
obs_opts     <- netcdf_opts(proj4_var = "projection_lcc", lon_var = NULL, lat_var = NULL)

fcst_dir      <- here("..", "data", "meps")
fcst_template <- "meps_det_2_5km_{YYYY}{MM}{DD}T{HH}Z.nc"
fcst_opts     <- netcdf_opts(ref_time_var = NA, z_var = "height0")

start_date <- 2022081500
end_date   <- start_date

Set the verification domain by reading in a forecast field and subgridding. Note that here the harpSpatialTeam version of harpVis is used that doesn’t include the get_domain() function, so the domain is extracted a slightly more cumbersome way.

dom <- read_forecast(
  date_times          = 2022081500,
  fcst_model          = "meps", 
  parameter           = "Pcp",
  lead_time           = 0,
  file_path           = fcst_dir,
  file_template       = fcst_template,
  file_format_opts    = fcst_opts,
  transformation      = "subgrid",
  transformation_opts = subgrid_opts(172, 472, 180, 480),
  return_data         = TRUE 
)[[1]][["meps_det"]][[1]] %>% 
  meteogrid::as.geodomain()

Plot the verif domain, just to make sure we got the correct one

plot(dom)

Now it’s time to do the verification using verifiy_spatial(). Note that param_find needs to be used to map to the variable names in the netcdf files.

verif <- verify_spatial(
  start_date       = start_date,
  end_date         = end_date, 
  parameter        = "AccPcp1h",
  model            = "meps",
  lead_time        = seq(0, 48),
  fc_file_path     = fcst_dir, 
  fc_file_template = fcst_template,
  fc_options = modify_opts(
    fcst_opts, 
    param_find = list(AccPcp1h = "precipitation_amount_acc")
  ),
  ob_file_path     = obs_dir,
  ob_file_template = obs_template,
  ob_options       = modify_opts(
    obs_opts, 
    param_find = list(AccPcp1h = "precipitation_amount")
  ),
  ob_accumulation  = "1h",
  verif_domain     = dom
)

Now that the verification is done, the output can be inspected.

verif

$bias
# A tibble: 48 × 5
   model prm          fcdate leadtime  bias
   <chr> <chr>         <dbl>    <dbl> <dbl>
 1 meps  AccPcp1h 1660521600     3600    NA
 2 meps  AccPcp1h 1660521600     7200    NA
 3 meps  AccPcp1h 1660521600    10800    NA
 4 meps  AccPcp1h 1660521600    14400    NA
 5 meps  AccPcp1h 1660521600    18000    NA
 6 meps  AccPcp1h 1660521600    21600    NA
 7 meps  AccPcp1h 1660521600    25200    NA
 8 meps  AccPcp1h 1660521600    28800    NA
 9 meps  AccPcp1h 1660521600    32400    NA
10 meps  AccPcp1h 1660521600    36000    NA
# … with 38 more rows

$mse
# A tibble: 48 × 5
   model prm          fcdate leadtime   mse
   <chr> <chr>         <dbl>    <dbl> <dbl>
 1 meps  AccPcp1h 1660521600     3600    NA
 2 meps  AccPcp1h 1660521600     7200    NA
 3 meps  AccPcp1h 1660521600    10800    NA
 4 meps  AccPcp1h 1660521600    14400    NA
 5 meps  AccPcp1h 1660521600    18000    NA
 6 meps  AccPcp1h 1660521600    21600    NA
 7 meps  AccPcp1h 1660521600    25200    NA
 8 meps  AccPcp1h 1660521600    28800    NA
 9 meps  AccPcp1h 1660521600    32400    NA
10 meps  AccPcp1h 1660521600    36000    NA
# … with 38 more rows

$mae
# A tibble: 48 × 5
   model prm          fcdate leadtime   mae
   <chr> <chr>         <dbl>    <dbl> <dbl>
 1 meps  AccPcp1h 1660521600     3600    NA
 2 meps  AccPcp1h 1660521600     7200    NA
 3 meps  AccPcp1h 1660521600    10800    NA
 4 meps  AccPcp1h 1660521600    14400    NA
 5 meps  AccPcp1h 1660521600    18000    NA
 6 meps  AccPcp1h 1660521600    21600    NA
 7 meps  AccPcp1h 1660521600    25200    NA
 8 meps  AccPcp1h 1660521600    28800    NA
 9 meps  AccPcp1h 1660521600    32400    NA
10 meps  AccPcp1h 1660521600    36000    NA
# … with 38 more rows

$SAL
# A tibble: 48 × 7
   model prm          fcdate leadtime       S       A     L
   <chr> <chr>         <dbl>    <dbl>   <dbl>   <dbl> <dbl>
 1 meps  AccPcp1h 1660521600     3600  0.287  -0.138  0.134
 2 meps  AccPcp1h 1660521600     7200 -0.342  -0.129  0.157
 3 meps  AccPcp1h 1660521600    10800  0.161   0.224  0.353
 4 meps  AccPcp1h 1660521600    14400  0.140   0.0358 0.143
 5 meps  AccPcp1h 1660521600    18000 -0.353  -0.0966 0.197
 6 meps  AccPcp1h 1660521600    21600  0.220   0.827  0.236
 7 meps  AccPcp1h 1660521600    25200  0.0660  0.790  0.526
 8 meps  AccPcp1h 1660521600    28800 -0.0956  0.340  0.101
 9 meps  AccPcp1h 1660521600    32400  0.0148  0.598  0.216
10 meps  AccPcp1h 1660521600    36000  0.144   0.530  0.527
# … with 38 more rows

$FSS
# A tibble: 1,344 × 7
   model prm          fcdate leadtime threshold scale   fss
   <chr> <chr>         <dbl>    <dbl>     <dbl> <dbl> <dbl>
 1 meps  AccPcp1h 1660521600     3600       0.1     0 0.664
 2 meps  AccPcp1h 1660521600     3600       0.1     1 0.744
 3 meps  AccPcp1h 1660521600     3600       0.1     2 0.782
 4 meps  AccPcp1h 1660521600     3600       0.1     4 0.830
 5 meps  AccPcp1h 1660521600     3600       0.1     8 0.886
 6 meps  AccPcp1h 1660521600     3600       0.1    12 0.920
 7 meps  AccPcp1h 1660521600     3600       0.1    20 0.954
 8 meps  AccPcp1h 1660521600     3600       1       0 0.526
 9 meps  AccPcp1h 1660521600     3600       1       1 0.670
10 meps  AccPcp1h 1660521600     3600       1       2 0.736
# … with 1,334 more rows

$NACT
# A tibble: 1,344 × 10
   model prm          fcdate leadt…¹ thres…² scale       a       b       c     d
   <chr> <chr>         <dbl>   <dbl>   <dbl> <dbl>   <dbl>   <dbl>   <dbl> <dbl>
 1 meps  AccPcp1h 1660521600    3600     0.1     0 0.0185  0.00806 0.0107  0.963
 2 meps  AccPcp1h 1660521600    3600     0.1     1 0.0187  0.00782 0.0105  0.963
 3 meps  AccPcp1h 1660521600    3600     0.1     2 0.0190  0.00751 0.0102  0.963
 4 meps  AccPcp1h 1660521600    3600     0.1     4 0.0196  0.00685 0.00951 0.964
 5 meps  AccPcp1h 1660521600    3600     0.1     8 0.0208  0.00560 0.00827 0.965
 6 meps  AccPcp1h 1660521600    3600     0.1    12 0.0219  0.00450 0.00718 0.966
 7 meps  AccPcp1h 1660521600    3600     0.1    20 0.0231  0.00320 0.00592 0.968
 8 meps  AccPcp1h 1660521600    3600     1       0 0.00217 0.00189 0.00203 0.994
 9 meps  AccPcp1h 1660521600    3600     1       1 0.00228 0.00178 0.00192 0.994
10 meps  AccPcp1h 1660521600    3600     1       2 0.00240 0.00165 0.00180 0.994
# … with 1,334 more rows, and abbreviated variable names ¹​leadtime, ²​threshold

There may be some issues with the likes of bias, mse, rmse etc. as they are all NA.

For plotting, the harpSpatialTeam version of harpVis should be updated to be compatible with these outputs.

#plot_spatial_verif(verif, FSS)
#plot_spatial_verif(verif, SAL)

The problem is in extracting the correct data frame (it is done with verif_data$score_name instead of verif_data[[score_name]]). There is a missing function when it comes to plotting the FSS, but SAL works if you just send it the data frame! However, as some of the functions used in the code are not properly namespaced, we need to make sure that the dplyr library is attached.

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

plot_spatial_verif(verif$SAL, SAL) +
  coord_equal()

Multiple leadtimes found: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Plotting score: SAL

Plots complete!

Coordinate system already present. Adding new coordinate system, which will replace the existing
one.

For FSS we have to go back to ggplot

ggplot(
  verif$FSS, 
  aes(
    factor(threshold), 
    factor(scale), 
    fill = fss, 
    label = sprintf("%1.2f", fss)
  )
) +
  geom_raster() +
  geom_text(size = 3) +
  facet_wrap(vars(leadtime), ncol = 8) +
  scico::scale_fill_scico(
    "FSS", palette = "bam", limits = c(0, 1)
  ) +
  coord_equal(expand = FALSE) +
  labs(
    x = "Threshold [mm]",
    y = "Neighbourhood Radius [grid boxes]"
  )

Note that here, the spatial scale is expressed as the neighbourhood radius rather than the neighbourhood length as in the harphub version. Also the lead time is always in seconds, which isn’t necessarily ideal but could easily be converted in the aesthetic setting (e.g. factor(leadtime / 3600))

The alternative method is to use a configuration file… There is a configuration file under harpspatialteam/conf

harpSpatial_read_config(here("conf", "my_conf.R"))

$by
[1] "12h"

$fc_accumulation
NULL

$fc_domain
NULL

$fc_file_format
NULL

$fc_file_path
[1] "/home/andrewts/harpSpatial_ww/harp-ww-feb2023/harpspatialteam/../data/meps"

$fc_file_template
[1] "meps_det_2_5km_{YYYY}{MM}{DD}T{HH}Z.nc"

$fc_interp_method
[1] "closest"

$fc_options
$fc_options$options_set
[1] "none"

$fc_options$param_find
$fc_options$param_find$AccPcp1h
[1] "precipitation_amount_acc"

$fc_options$param_find$AccPcp3h
[1] "precipitation_amount_acc"

$fc_options$param_find$AccPcp6h
[1] "precipitation_amount_acc"

$fc_options$param_find$AccPcp12h
[1] "precipitation_amount_acc"

$fc_options$param_find$AccPcp24h
[1] "precipitation_amount_acc"


$fc_options$proj4_var
[1] "projection_lambert"

$fc_options$proj4_att
[1] "proj4"

$fc_options$proj4
NULL

$fc_options$x_dim
[1] "x"

$fc_options$y_dim
[1] "y"

$fc_options$lon_var
[1] "longitude"

$fc_options$lat_var
[1] "latitude"

$fc_options$x_rev
[1] FALSE

$fc_options$y_rev
[1] FALSE

$fc_options$dx
NULL

$fc_options$dy
NULL

$fc_options$z_var
[1] "height0"

$fc_options$member_var
[1] NA

$fc_options$time_var
[1] "time"

$fc_options$ref_time_var
[1] NA

$fc_options$force_param_name
[1] FALSE


$lead_time
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
[26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

$lt_unit
[1] "h"

$members
NULL

$ob_accumulation
[1] "1h"

$ob_domain
NULL

$ob_file_format
NULL

$ob_file_path
[1] "/home/andrewts/harpSpatial_ww/harp-ww-feb2023/harpspatialteam/../data/met_analysis"

$ob_file_template
[1] "met_analysis_1_0km_nordic_{YYYY}{MM}{DD}T{HH}Z.nc"

$ob_interp_method
[1] "closest"

$ob_options
$ob_options$options_set
[1] "none"

$ob_options$param_find
$ob_options$param_find$AccPcp1h
[1] "precipitation_amount"

$ob_options$param_find$AccPcp3h
[1] "precipitation_amount"

$ob_options$param_find$AccPcp6h
[1] "precipitation_amount"

$ob_options$param_find$AccPcp12h
[1] "precipitation_amount"

$ob_options$param_find$AccPcp24h
[1] "precipitation_amount"


$ob_options$proj4_var
[1] "projection_lcc"

$ob_options$proj4_att
[1] "proj4"

$ob_options$proj4
NULL

$ob_options$x_dim
[1] "x"

$ob_options$y_dim
[1] "y"

$ob_options$lon_var
NULL

$ob_options$lat_var
NULL

$ob_options$x_rev
[1] FALSE

$ob_options$y_rev
[1] FALSE

$ob_options$dx
NULL

$ob_options$dy
NULL

$ob_options$z_var
[1] NA

$ob_options$member_var
[1] NA

$ob_options$time_var
[1] "time"

$ob_options$ref_time_var
[1] NA

$ob_options$force_param_name
[1] FALSE


$sal_options
$sal_options$thresh_scale
[1] 15

$sal_options$min_rain
[1] 0.1

$sal_options$same_threshold
[1] FALSE

$sal_options$maxobj
[1] 1000


$sqlite_file
NULL

$sqlite_path
NULL

$thresholds
[1]  0.1  1.0  5.0 10.0

$use_mask
[1] FALSE

$verif_domain
301 x 301 domain
Projection summary:
proj= lcc 
centre = ( 10.55944 , 58.63803 )

$window_sizes
[1]  0  1  2  4  8 12 20

12h

/home/andrewts/harpSpatial_ww/harp-ww-feb2023/harpspatialteam/../data/meps

meps_det_2_5km_{YYYY}{MM}{DD}T{HH}Z.nc

closest

none

precipitation_amount_acc

precipitation_amount_acc

precipitation_amount_acc

precipitation_amount_acc

precipitation_amount_acc

projection_lambert

proj4

x

y

longitude

latitude

height0

time

h

1h

/home/andrewts/harpSpatial_ww/harp-ww-feb2023/harpspatialteam/../data/met_analysis

met_analysis_1_0km_nordic_{YYYY}{MM}{DD}T{HH}Z.nc

closest

none

precipitation_amount

precipitation_amount

precipitation_amount

precipitation_amount

precipitation_amount

projection_lcc

proj4

x

y

time

verif_from_conf <- verify_spatial(
  start_date = start_date,
  end_date   = end_date,
  parameter  = "AccPcp1h",
  model      = "meps"
)

We can check that we got the same result using identical()

identical(verif, verif_from_conf)

[1] TRUE