Read and transform examples - part 2
Andrew Singleton
16/02/2022
library(harp)
library(here)
library(dplyr)
library(scico)
library(tidyr)
library(meteogrid)
library(forcats)Show the built in harp parameters
- Note that for FA files, parameters are hardcoded.
show_param_defs()## # A tibble: 60 × 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 sst Sea surface temperature
## 36 t Air temperature
## 37 t0m Skin temperature
## 38 t2m Air temperature at 2m above the ground
## 39 td Dew point temperature
## 40 td2m Dew point temperature at 2m above the ground
## 41 tmax Maximum air temperature at 2m above the ground
## 42 tmin Minimum air temperature at 2m above the ground
## 43 topo Height of topography above sea level
## 44 u Wind speed in u direction
## 45 u10m Wind speed in u direction at 10m above the ground
## 46 ugust Wind gust in U direction
## 47 ugust10m Wind gust at 10m above the ground in U direction
## 48 v Wind speed in v direction
## 49 v10m Wind speed in v direction at 10m above the ground
## 50 vgust Wind gust in V direction
## 51 vgust10m Wind gust at 10m above the ground in V direction
## 52 vis Visibility in air
## 53 w Vertical (upward) wind speed
## 54 wd Wind direction calculated from U and V winds
## 55 wd10m Wind direction at 10m above the ground calculated from U and V…
## 56 wg Wind gust calculated from U and V gusts
## 57 wg10m Wind gust at 10m above the ground calculated from U gust and V…
## 58 ws Wind speed calculated from U and V winds
## 59 ws10m Wind speed at 10m above the ground calculated from U and V win…
## 60 z Geopotentialshow_param_defs("grib")## # A tibble: 50 × 2
## name grib_name
## <chr> <chr>
## 1 caf tcc
## 2 cape cape
## 3 cchigh hcc
## 4 cclow lcc
## 5 ccmed mcc
## 6 cctot tcc
## 7 d wdir
## 8 d10m wdir
## 9 g fg
## 10 g10m fg
## 11 gh gh
## 12 lsm lsm
## 13 pcp tp
## 14 pmsl msl
## 15 pressure pres
## 16 psfc pres
## 17 q q
## 18 q2m c(2q, q)
## 19 rh r
## 20 rh2m c(2r, r)
## 21 s ws
## 22 s10m c(10si, ws)
## 23 sea_ice icec
## 24 sfc_geo z
## 25 sst sst
## 26 t t
## 27 t0m c(0t, t)
## 28 t2m c(2t, t)
## 29 td td
## 30 td2m td
## 31 tmax c(mx2t, mxt)
## 32 tmin c(mn2t, mnt)
## 33 topo orog
## 34 u u
## 35 u10m c(10u, u)
## 36 ugust ugst
## 37 ugust10m ugst
## 38 v v
## 39 v10m c(10v, v)
## 40 vgust vgst
## 41 vgust10m vgst
## 42 vis vis
## 43 w w
## 44 wd list(u = u, v = v)
## 45 wd10m list(u = c(10u, u), v = c(10v, v))
## 46 wg list(u = ugst, v = vgst)
## 47 wg10m list(u = ugst, v = vgst)
## 48 ws list(u = u, v = v)
## 49 ws10m list(u = c(10u, u), v = c(10v, v))
## 50 z zshow_param_defs("netcdf")## # A tibble: 51 × 2
## name netcdf_name
## <chr> <chr>
## 1 caf cloud_area_fraction
## 2 cape specific_convective_available_potential_energy
## 3 cbase cloud_base_altitude
## 4 cchigh high_type_cloud_area_fraction
## 5 cclow low_type_cloud_area_fraction
## 6 ccmed medium_type_cloud_area_fraction
## 7 cctot cloud_area_fraction
## 8 cin atmosphere_convective_inhibition
## 9 d wind_direction
## 10 d10m wind_direction
## 11 g wind_speed_of_gust
## 12 g10m wind_speed_of_gust
## 13 gh geopotential_height
## 14 lsm land_binary_mask
## 15 pcp precipitation_amount_acc
## 16 pmsl air_pressure_at_sea_level
## 17 psfc surface_air_pressure
## 18 q specific_humidity
## 19 q2m specific_humidity_2m
## 20 rh relative_humidity
## 21 rh2m relative_humidity_2m
## 22 s wind_speed
## 23 s10m wind_speed
## 24 sea_ice ga_icec_102
## 25 sfc_geo surface_geopotential
## 26 sst sea_surface_temperature
## 27 t air_temperature
## 28 t0m air_temperature_0m
## 29 t2m air_temperature_2m
## 30 td dew_point_temperature
## 31 td2m dew_point_temperature_2m
## 32 tmax air_temperature_max
## 33 tmin air_temperature_min
## 34 topo altitude
## 35 u x_wind
## 36 u10m x_wind_10m
## 37 ugust x_wind_gust
## 38 ugust10m x_wind_gust_10m
## 39 v y_wind
## 40 v10m y_wind_10m
## 41 vgust y_wind_gust
## 42 vgust10m y_wind_gust_10m
## 43 vis visibility_in_air
## 44 w upward_air_velocity
## 45 wd list(u = x_wind, v = y_wind)
## 46 wd10m list(u = x_wind_10m, v = y_wind_10m)
## 47 wg list(u = x_wind_gust, v = y_wind_gust)
## 48 wg10m list(u = x_wind_gust_10m, v = y_wind_gust_10m)
## 49 ws list(u = x_wind, v = y_wind)
## 50 ws10m list(u = x_wind_10m, v = y_wind_10m)
## 51 z geopotentialshow_param_defs("v")## # A tibble: 35 × 2
## name v_name
## <chr> <chr>
## 1 accpcp12h PE12
## 2 accpcp1h PE1
## 3 accpcp24h PE24
## 4 accpcp3h PE3
## 5 accpcp6h PE6
## 6 cbase CH
## 7 cc_below_7500 N75
## 8 cchigh HC
## 9 cclow LC
## 10 ccmed MC
## 11 cctot NN
## 12 d DD
## 13 d10m DD
## 14 g10m GW
## 15 gh FI
## 16 gmax GX
## 17 pcp PE
## 18 pmsl PS
## 19 pressure PP
## 20 psfc PSS
## 21 q QQ
## 22 q2m QQ
## 23 rh RH
## 24 rh2m RH
## 25 s FF
## 26 s10m FF
## 27 smax WX
## 28 t TT
## 29 t2m TT
## 30 td TD
## 31 td2m TD
## 32 tmax TX
## 33 tmin c(TM, TN)
## 34 topo FI
## 35 vis VIshow_param_defs("wrf")## # A tibble: 22 × 2
## name wrf_name
## <chr> <chr>
## 1 caf CLDFRA
## 2 cctot CLDFRA
## 3 gh Z
## 4 lsm LANDMASK
## 5 pcp list(RAINNC, RAINC)
## 6 q QVAPOR
## 7 q2m Q2
## 8 sst SST
## 9 t T
## 10 t0m TSK
## 11 t2m T2
## 12 topo HGT
## 13 u U
## 14 u10m U10
## 15 v V
## 16 v10m V10
## 17 w W
## 18 wd list(u = U, v = V)
## 19 wd10m list(u = U10, v = V10)
## 20 ws list(u = U, v = V)
## 21 ws10m list(u = U10, v = V10)
## 22 z PHRead in U10m and V10m
and get the wind speed using “ws10m”
template <- "{fcst_model}/wind/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR3}/fc{YYYY}{MM}{DD}{HH}+{LDT3}grib"
ws10m <- read_forecast(
dttm = 2021012100,
fcst_model = "AAEPS",
parameter = "WS10m",
lead_time = 0,
members = 0,
file_path = here("data/grib"),
file_template = template,
return_data = TRUE
)
u10m <- read_forecast(
dttm = 2021012100,
fcst_model = "AAEPS",
parameter = "U10m",
lead_time = 0,
members = 0,
file_path = here("data/grib"),
file_template = template,
return_data = TRUE
)
v10m <- read_forecast(
dttm = 2021012100,
fcst_model = "AAEPS",
parameter = "V10m",
lead_time = 0,
members = 0,
file_path = here("data/grib"),
file_template = template,
return_data = TRUE
)
u10m## ::ensemble gridded forecast:: # A tibble: 1 × 11
## fcst_model fcst_dttm lead_time parameter valid_dttm
## <chr> <dttm> <dbl> <chr> <dttm>
## 1 AAEPS 2021-01-21 00:00:00 0 U10m 2021-01-21 00:00:00
## # ℹ 6 more variables: step_range <chr>, level_type <chr>, level <int>,
## # units <chr>, fcst_cycle <chr>, AAEPS_mbr000 <geolist>We can check the results by plotting
plot_field(u10m, plot_col = AAEPS_mbr000)
plot_field(v10m, plot_col = AAEPS_mbr000)
plot_field(ws10m, plot_col = AAEPS_mbr000)
plot_field(
mutate(
inner_join(
rename(u10m, u = AAEPS_mbr000),
select(rename(v10m, v = AAEPS_mbr000), -parameter)
),
AAEPS_mbr000 = geolist(sqrt(u ^ 2 + v ^ 2)),
parameter = "sqrt(u ^ 2 + v ^ 2)"
),
plot_col = AAEPS_mbr000
)## Joining with `by = join_by(fcst_model, fcst_dttm, lead_time, valid_dttm,
## step_range, level_type, level, units, fcst_cycle)`
Read a radar analysis (note this is over 1Gb of data - check you have enough memory)
template <- paste0(
"{fcst_model}/{YYYY}/{MM}/norway.mos.sri-acrr-1h.noclass-clfilter-",
"vpr-clcorr-block.utm33-1000.{YYYY}{MM}{DD}.nc"
)
radar <- read_analysis(
dttm = seq_dttm(2022021200, 2022021323),
analysis_model = "radar",
parameter = "lwe_precipitation_rate",
file_template = template,
file_path = here("data/netcdf"),
file_format_opts = netcdf_opts(
proj4_var = "projection_utm",
x_dim = "Xc",
y_dim = "Yc",
lon_var = "lon",
lat_var = "lat",
y_rev = TRUE
)
)plot_field(
radar$anl[[23]],
breaks = c(0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16),
palette = c(
"transparent", scico(255, palette = "oslo", direction = -1)
)
)
Read point observations
Use read_obs() in the same way as
read_forecast(). Currently only works for vobs, so the file
template is set to vobs by default - only need the dates and the
path.
obs <- read_obs(
dttm = seq_dates(2019021700, 2019021723),
file_path = here("data/vobs"),
return_data = TRUE
)
obs## $synop
## # A tibble: 80,598 × 24
## valid_dttm SID lat lon elev CCtot D10m S10m T2m Td2m
## <dttm> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2019-02-17 00:00:00 1001 70.9 -8.67 9 8 202 8.8 270. 268.
## 2 2019-02-17 00:00:00 1002 80.1 16.2 8 NA 330 13 248. 246
## 3 2019-02-17 00:00:00 1003 77 15.5 10 1.04 30 2 252. 245.
## 4 2019-02-17 00:00:00 1006 78.3 22.8 14 NA 20 3 247. 244.
## 5 2019-02-17 00:00:00 1007 78.9 11.9 8 NA 350 5 250. NA
## 6 2019-02-17 00:00:00 1008 78.2 15.5 27 7.04 298 7.6 251. 246
## 7 2019-02-17 00:00:00 1009 80.7 25.0 5 NA 330 12 244. 241.
## 8 2019-02-17 00:00:00 1010 69.3 16.1 13 NA 340 11.3 274. 268.
## 9 2019-02-17 00:00:00 1011 80.1 31.5 9 NA 352 11.4 240. 236.
## 10 2019-02-17 00:00:00 1013 78.1 13.6 -99 NA NA NA NA NA
## # ℹ 80,588 more rows
## # ℹ 14 more variables: RH2m <dbl>, Q2m <dbl>, Ps <dbl>, Pmsl <dbl>, vis <dbl>,
## # AccPcp3h <dbl>, AccPcp6h <dbl>, AccPcp24h <dbl>, N75 <int>, CClow <int>,
## # Cbase <int>, AccPcp1h <dbl>, Gmax <dbl>, AccPcp12h <dbl>
##
## $synop_params
## parameter accum_hours units
## 1 CCtot 0 oktas
## 2 D10m 0 degrees
## 3 S10m 0 m/s
## 4 T2m 0 K
## 5 Td2m 0 K
## 6 RH2m 0 percent
## 7 Q2m 0 kg/kg
## 8 Ps 0 hPa
## 9 Pmsl 0 hPa
## 10 vis 0 m
## 11 AccPcp3h 3 kg/m^2
## 12 AccPcp6h 6 kg/m^2
## 13 AccPcp24h 24 kg/m^2
## 14 N75 0 oktas
## 15 CClow 0 oktas
## 16 Cbase 0 m
## 17 AccPcp1h 1 kg/m^2
## 18 Gmax 0 m/s
## 19 AccPcp12h 12 kg/m^2
##
## $temp
## # A tibble: 37,872 × 13
## valid_dttm SID lat lon elev p Z T RH D
## <dttm> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2019-02-17 00:00:00 1001 70.9 -8.67 9 1000 64 269. 83.3 202
## 2 2019-02-17 00:00:00 1001 70.9 -8.67 9 925 672 264. 100 205
## 3 2019-02-17 00:00:00 1001 70.9 -8.67 9 850 1322 262. 86.4 161
## 4 2019-02-17 00:00:00 1001 70.9 -8.67 9 700 2792 254. 32 251
## 5 2019-02-17 00:00:00 1001 70.9 -8.67 9 500 5224 240. 8.7 270
## 6 2019-02-17 00:00:00 1001 70.9 -8.67 9 400 6762 230. 28.1 275
## 7 2019-02-17 00:00:00 1001 70.9 -8.67 9 300 8656 218 38.6 269
## 8 2019-02-17 00:00:00 1001 70.9 -8.67 9 250 NA NA NA NA
## 9 2019-02-17 00:00:00 1001 70.9 -8.67 9 200 11220 218. 2.9 262
## 10 2019-02-17 00:00:00 1001 70.9 -8.67 9 150 13048 217. 2.5 257
## # ℹ 37,862 more rows
## # ℹ 3 more variables: S <dbl>, Q <dbl>, Td <dbl>
##
## $temp_params
## parameter accum_hours units
## 1 p 0 hPa
## 2 Z 0 m
## 3 T 0 K
## 4 RH 0 percent
## 5 D 0 degrees
## 6 S 0 m/s
## 7 Q 0 kg/kg
## 8 Td 0 KWe can easily plot some synop data on a map
ggplot(drop_na(obs$synop, T2m)) +
geom_polygon(
aes(x = long, y = lat, group = group),
map_data("world"),
fill = NA,
colour = "grey30"
) +
geom_point(aes(x = lon, y = lat, colour = T2m)) +
coord_equal(
xlim = c(0, 40),
ylim = c(40, 80)
) +
scale_colour_viridis_c(limits = c(243, 303))
Transformations
harp’s read functions allow you to transform data as it is read in.
That is to say, it can interpolate gridded data to geographic points,
regrid and reproject gridded data and make a vertical cross section from
3d gridded data. These are activated through the
transformation argument to read_grid(),
read_forecast() and read_analysis(). Special
options for each transformation are passed via the
transformation_opts argument, with each transforamtion
having it’s own function to set those options.
Interpolation to points.
By default this is done using harp’s built in station list
station_list## # 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 rowsAlso, if there is sufficient data available (surface geopotential or topography height), 2m temperature can be corrected to take account of the difference in model elevation and the actual elevation of the station.
Use transformation = "interpolate" and
interpolate_opts(). All options have default values, so it
isn’t always necessary to set anything
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] FALSEWe can simply run read_forecast() as before and
interpolate to stations to the default stations using `transformation =
“interpolate”
template <- "{fcst_model}/{YYYY}/{MM}/{DD}/{HH}/mbr{MBR3}/fc{YYYY}{MM}{DD}{HH}+{LDT3}grib"
t2m <- read_forecast(
dttm = seq_dttm(2021012100, 2021012200, "1d"),
fcst_model = "AAEPS",
parameter = "T2m",
lead_time = seq(0, 24, 3),
members = seq(0, 3),
transformation = "interpolate",
file_path = here("data/grib"),
file_template = template,
return_data = TRUE
)## Warning: 'transformation_opts' not set for transformation = 'interpolate'.
## Using default interpolate_opts()## No stations specified. Using default stations: 'station_list'## Initializing interpolate weights using 'sfc_geo'.## Reading /home/andrewts/harp-documentation/md/harp_training_2022/data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000grib## 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/harp-documentation/md/harp_training_2022/data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000grib## Error: 2m temperature height correction is selected, but cannot find 'sfc_geo' in /home/andrewts/harp-documentation/md/harp_training_2022/data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000grib.
## You probably need to set transformation_opts = interpolate_opts(clim_file = '<clim_file>', clim_param = '<clim_param>'), or switch off 2m temperature correction with transformation_opts = interpolate_opts(correct_t2m = FALSE).We don’t actually have the correct data available to do the height
correction of 2m temperature so we can switch it off via
interpolate_opts()
t2m <- read_forecast(
dttm = seq_dttm(2021012100, 2021012200, "1d"),
fcst_model = "AAEPS",
parameter = "T2m",
lead_time = seq(0, 3, 3),
members = seq(0, 1),
transformation = "interpolate",
transformation_opts = interpolate_opts(
correct_t2m = FALSE
),
file_path = here("data/grib"),
file_template = template,
return_data = TRUE
)## Error : None of the requested data could be read from grib file: /home/andrewts/harp-documentation/md/harp_training_2022/data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000gribt2m## ::ensemble point forecast:: # A tibble: 48 × 15
## fcst_model fcst_dttm lead_time parameter valid_dttm
## <chr> <dttm> <dbl> <chr> <dttm>
## 1 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 2 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 3 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 4 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 5 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 6 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 7 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 8 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 9 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## 10 AAEPS 2021-01-21 00:00:00 0 T2m 2021-01-21 00:00:00
## # ℹ 38 more rows
## # ℹ 10 more variables: step_range <chr>, level_type <chr>, level <int>,
## # units <chr>, SID <int>, lat <dbl>, lon <dbl>, fcst_cycle <chr>,
## # AAEPS_mbr000 <dbl>, AAEPS_mbr001 <dbl>Regridding to a different grid
We can use any geofield (or geodomain) as a reference to regrid the
data to. Now we set transformation = "regrid") and we use
regrid_opts() to set the options. Note that
regrid_opts() must have a new_domain
argument.
Let’s read the radar data in again, but this time regrid it to the domain from the netcdf file we had with precipitation.
domain <- get_domain(
read_grid(
here("data/netcdf/meps/2022/02/12/meps_lagged_6_h_subset_2_5km_20220212T00Z.nc"),
"pcp",
lead_time = 0,
members = 0,
file_format_opts = netcdf_opts(
ref_time_var = "forecast_reference_time",
member_var = "ensemble_member",
z_var = "height0"
)
)
)## Reading /home/andrewts/harp-documentation/md/harp_training_2022/data/netcdf/meps/2022/02/12/meps_lagged_6_h_subset_2_5km_20220212T00Z.ncregrid_opts(new_domain = domain)## $new_domain
## 121 x 206 domain
## Projection summary:
## proj= lcc
## NE = ( 9.38741 , 62.7829 )
## SW = ( 5.104273 , 57.8871 )
##
## $method
## [1] "nearest"
##
## $clim_file
## NULL
##
## $clim_file_format
## NULL
##
## $clim_file_opts
## NULL
##
## $clim_param
## [1] "sfc_geo"
##
## $weights
## NULL
##
## $keep_raw_data
## [1] FALSENow we can simply regrid the radar data as we read it in
template <- paste0(
"{fcst_model}/{YYYY}/{MM}/norway.mos.sri-acrr-1h.noclass-clfilter-",
"vpr-clcorr-block.utm33-1000.{YYYY}{MM}{DD}.nc"
)
radar <- read_analysis(
dttm = seq_dttm(2022021300, 2022021323),
analysis_model = "radar",
parameter = "lwe_precipitation_rate",
transformation = "regrid",
transformation_opts = regrid_opts(new_domain = domain),
file_template = template,
file_path = here("data/netcdf"),
file_format_opts = netcdf_opts(
proj4_var = "projection_utm",
x_dim = "Xc",
y_dim = "Yc",
lon_var = "lon",
lat_var = "lat",
y_rev = TRUE
)
)## Error : None of the requested parameters found in file: /home/andrewts/harp-documentation/md/harp_training_2022/data/netcdf/radar/2022/02/norway.mos.sri-acrr-1h.noclass-clfilter-vpr-clcorr-block.utm33-1000.20220213.ncAnd now we can plot
ggplot(radar) +
geom_georaster(aes(geofield = anl)) +
geom_polygon(
aes(x, y, group = group),
get_map(domain),
fill = NA,
colour = "grey70"
) +
scale_fill_scico(
"mm",
palette = "oslo",
direction = -1,
limits = c(0.125, NA),
na.value = NA,
trans = "log",
breaks = c(0.125, 0.25, 0.5, 1, 2, 4, 8, 16),
labels = c(0.125, 0.25, 0.5, 1, 2, 4, 8, 16)
) +
facet_wrap(vars(fct_inorder(format(valid_dttm, "%a %H:%M")))) +
coord_equal(expand = FALSE) +
theme_void() +
theme(
plot.margin = unit(c(1, 1, 1, 1), "cm"),
panel.background = element_rect(fill = NA, colour = "grey30")
)## Warning: Transformation introduced infinite values in discrete y-axis## Warning: Removed 384478 rows containing missing values (`geom_raster()`).
Vertical Cross Sections
Vertical cross sections are extracted using
transformation = "xsection" and
transformation_opts = xsection_opts(a = ..., b = ...) where
a and b are the lon-lat locations of the ends
of the cross section.
xs <- read_grid(
here("data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000grib"),
"T",
vertical_coordinate = "model",
transformation = "xsection",
transformation_opts = xsection_opts(
a = c(7.76, 78.49), b = c(28.25, 79.29)
)
)## Reading /home/andrewts/harp-documentation/md/harp_training_2022/data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000grib## Computing interpolation weights for xsectioning.A data frame is returned (or in the case of
read_forecast a data frame with a column having a list of
data frames), so can easily plot with geom_raster()
ggplot(xs, aes(distance, level_number, fill = value)) +
geom_raster() +
scale_fill_scico("K", palette = "batlow") +
scale_y_reverse() +
coord_cartesian(expand = FALSE)
The cross sections are oriented from a to b, so if we reverse a and b everything will be the other way around.
xs <- read_grid(
here("data/grib/AAEPS/2021/01/21/00/mbr000/fc2021012100+000grib"),
"T",
vertical_coordinate = "model",
transformation = "xsection",
transformation_opts = xsection_opts(
b = c(7.76, 78.49), a = c(28.25, 79.29)
)
)ggplot(xs, aes(distance, level_number, fill = value)) +
geom_raster() +
scale_fill_scico("K", palette = "batlow") +
scale_y_reverse() +
coord_cartesian(expand = FALSE)