Python/pflexpart
Introduction to Working with pflexible
See the instructions and documentation for working with "pflexible" at bitbucket.
Recipes
This is a place for keeping track of the current process for doing various things with pflexpart. If you are working through a plotting routine, or otherwise doing some processing with pflexpart, please feel free to contribute some notes here.
Getting output from a model run for a known time/species/level
H = pf.Header('path/to/output')
# read a grid for a known time / nspec
# NOTE:these are indexes to H.available_dates and H.species so 1 is not _001
H.read_grid(time_ret=0, nspec_ret=0 )
In the code above, time_ret and nspec_ret are indices to H.available_dates and H.species. It is important to note that nspec_ret is *not* equal to the number used in the filename of the flexpart output (e.g. _001, etc.).
H.read_grid will create a new attribute to the Header. This new object ('FD') is a Python dictionary. It is keyed by tuples that represent the species and time for the request. See H.FD.keys() to get an overview of what is available.
If you know approximately the date you want to plot, and quick way to get the index is to use the Header method:
H.closest_date('2009-02-13 10:00:00', fmt="%Y-%m-%d %H:%M:%S")
Note you can pass any format string representation of the time you like, so long as you pass an accompanying string format [fmt] argument. Also, you can simply pass a datetime object.
The next step is to get the data from the FD (which stands for Flexpart Data, by the way). Each 'grid' for each species/date pair will have shape (X , Y, Z, nspec) but these days nspec is almost always 1, and *not* related to above nspec_ret (sorry, just the way it is -- legacy).
So, we want to get the 'grid' for our request:
flexdata = H.FD[ (nspec, date) ]
The next step is to calculate the total column and break the grid into 'slabs' that we want for our analysis, and add it to the flexdata dictionary:
flexdata['slabs'] = pf.get_grid(H, flexdata['grid'])
Now 'slabs' is another dictionary. This time keyed by {0 ... N} where N is the number of outheights+1. The flexdata['slabs'][0] grid will be the Total Column, flexdata['slabs'][1] is the footprint:
tc = flexdata['slabs'][0] tc.shape
