...
As an exercise to put all of this together, we will write a new macro to do the followingcompute the precipitation rate in mm per hour over a sub-area of the globe for each time step. The steps will be:
- 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_diff
- 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)
- compute the mean value over a sub-area for each field in
precip_diff
(use theintegrate()
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
...
Code Block | ||
---|---|---|
| ||
XXXX 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_diff = 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 = grib_get_long(precip_diff[i], 'validityDate')
t = grib_get_long(precip_diff[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]
means = means * 1000 # remember to scale up from m to mm!
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) |
...