...
- compute the 'period' precipitation from precip.grib - this is what we already did earlier, so it's done! Just copy the code to your new macro and change the result variable name to
precip_perioddiff
- loop through the fields in the original fieldset and for each:
- get the date and time of the forecast step
- combine these into a Metview date variable
- add it to a list
- use syntax similar to the line of code used to compute the 'period' precipitation to find the differences between the times of adjacent fields (ok, we know it's 3 hours, but in theory it could be anything)extract a sub-area from the computed precipitation data
- compute the mean value over that a sub-area for each field in
precip_diff
(use theaverageintegrate()
function) - scale up from metres (as the data are stored) to mm by multiplying by 1000
- convert this into a rate, mm per hour, using the time differences computed earlier
...
You can get the date and time of a field as numbers like this:
Code Block | ||
---|---|---|
| ||
field_valid_dated = grib_get_stringlong(precip_3hdiff[i], 'validityDate') field_valid_timet = grib_get_stringlong(precip_3hdiff[i], 'validityTime') |
then combine those numbers into proper date variables.
...
Code Block | ||
---|---|---|
| ||
dates = nil for i = .... do ddt = ..... # construct a date/time variable dates = dates & [ddt] end for |
The integrate()
function can be called in a number of ways, but we will use it like thisTo extract a sub-area, use the following syntax (derived from using the GRIB Filter icon's Area parameter):
Code Block | ||
---|---|---|
| ||
f means = read(data: fintegrate(precip_diffs, area:[N,W,S,E]) |
The average()
function returns to compute the mean values over a sub-area of each field. The result is a list of values, a mean for each field. You can directly multiply a list variable by a number to obtain a new list where each element has been multiplied.
The final calculation requires converting the time intervals into hours (because if the time difference between two steps is 7 hours, then the rate of precip per hour is the mean precip value divided by 7).
Code Block | ||
---|---|---|
| ||
XXXXXXXX remove once tested! # read the original data file with accumulated precipitation precip = read("precip.grib") # compute the 'period' precipitation - the difference between steps n = count(precip) # the number of fields in the fieldset precip_perioddiff = precip[2, n] - precip[1, n-1] # extract the dates and times from the original fields and put into a list dates = nil for i = 1 to count(precip) do d_valid_date = grib_get_stringlong(precip_3hdiff[i], 'validityDate') t = grib_get_stringlong(precip_3hdiff[i], 'validityTime') # ASSUME it's going to be in hours dt = date(d) + hour(t/100) dates = dates & [dt] end for print(dates) # compute the differences between the valid times date_diffs = dates[2, n] - dates[1, n-1] print(date_diffs) means = integrate(precip_diff, [60,-15,50,5]) # [N,W,S,E] print(means) precip_rates_per_hour = means / (date_diffs*24) # date_diffs*24 gives the number of hours in each interval print(precip_rates_per_hour) |