Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

Code Block
languagepy
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 the 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 at a  particular location for each time period.We'll do it step by step.

Compute the 'period' precipitation from precip.grib

This is what we already did earlier, so it's done! Just make a copy your earlier macro, compute_precip, and call it precip_rate. Change the result variable name to precip_diff to make it more generic. Remove the return line, as we want to use this fieldset, not return it.

Construct a loop to go through the fields

Now, create an empty list (dates = nil). We will add each date variable to it when we loop through the fields.

...

Code Block
languagepy
dates = nil
for i = 1 to n do
   print(i) # we will put proper code here later!
end for

Extract the date and time from each field

You can get the date and time of a field as numbers like this, inside the loop, where i is the field index:

...

Info
The grib_get() functions are general-purpose functions to get pieces of meta-data from a GRIB field, specified by keys such as 'validityDate'. The Grib Examiner can help you find the available keys.

Convert these numbers into a date variable

You will need to divide the time variable (t) to convert it into a fraction of a day before adding it to the date variable (d). A time of 12:00 is returned as 1200, so we need to divide by 100 to get it into hours. The hours() function then converts it into a fraction of a day. The result can be converted with the date() function into a proper date; call the resulting variable dt and print it to check that it is a full date variable.

Code Block
languagepy
dt = date(d) + hour(t/100)

Add the date to the list

If the date variable is dt, then we add it to the list like this (inside the loop):

...

Code Block
languagepy
date_diffs_in_hours = date_diffs*24

Extract the point value for each field in precip_diff

Use the nearest_gridpoint() function on the precip_diff fieldset. It returns a list of values, one for each field. Choose a location with some high precipitation.

...

The result is a list of values, a 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 - do this to scale from metres to mm.

Compute precipitation rate in mm per hour

The final calculation requires converting the data values into mm per hour - divide the list of precipitation values by the list of time differences, which are in hours.

...