...
Info | ||
---|---|---|
| ||
Multiple data provided in several data icons are overlaid according to the overlay setting in the current view. |
Hovmoeller?
Precipitation
Precipitation data provides an interesting challenge. Precipitation fields in MARS are stored as accumulated fields. Visualise the supplied precip.grib icon with the precip_shade visdef. The first field is empty (check using the Cursor Data). The first field has a step of 0, meaning that it contains the total precipitation accumulated between the run time and the run time plus step. Since these are the same, there is no accumulated precipitation! Subsequent steps show more and more precipitation (the amount accumulated over 3, 6, 9, etc hours).
...
Code Block | ||
---|---|---|
| ||
for d = 2015-01-01 to 2015-03-01 do print(d) # each step is 1 day end for for d = 2015-01-01 to 2015-03-01 by 2 do print(d) # each step is 2 days end for for d = 2015-01-01 to 2015-03-01 by hour(6) do print(d) # each step is 6 hours end for |
Computing mean precipitation rate
...
at a point
As an exercise to put all of this together, we will write a new macro to compute the precipitation rate in mm per hour over a sub-area of the globe at a particular location 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 (which was initialised to
nil
before the loop)
- 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 extract the mean point value over a sub-area for each field in
precip_diff
(use theintegratenearest_gridpoint()
function). Choose a location with some high precipitation - 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 | ||
---|---|---|
| ||
dates = nil for i = .... do dt = ..... # construct a date/time variable dates = dates & [dt] end for |
The integratenearest_gridpoint()
function can be called in a number of ways, but we will use it like this:
Code Block | ||
---|---|---|
| ||
means values = integratenearest_gridpoint(precip_diffs, [N,W,S,E])lat, lon) |
to compute the mean values over a sub-area of each field. The result is a list of values, a mean value 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 | ||
---|---|---|
| ||
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) |