RPN standard files
The RPN standard file format is a file format developed by "Le Centre de recherche en prévision numérique" (RPN) to store weather and climate related data. Usually, there is one 2-D (latitude and longitude) record for each level, variable and timestep.
For the end-user, the appearance of the RPN is that of a flat-file database structure. A single file contains numerous records which hold fields on 1, 2 or 3D data grids and self-describing information about the fields (header).
A thorough set of FST-related routines are contained in the RPN Library (armnlib). These routines provide easy access to the file contents. They can be called from FORTRAN or C routines.
A complete description of the FST file format and the utilities which can be used to manipulate and view their contents is available online at www.cmc.ec.gc.ca/rpn/utilities.
Sometimes the files are archived in an cmc-archive file (*.ca). If this is the case they need to be unarchived first.
The maximum size of RPN file is a little more than 8 GB!!! If you are trying to write a file larger than the maximum size the file will not get closed. Meaning no application will be able to open it and it can be considered as "corrupted". The usual error message one get when trying to open such a file is "file (unit=...) currently used by another application in write mode".
Looking at RPN files
vl3, vl4, vl5 -> list variable with description
At UQAM, there is also a set of 'vl' (variable list) commands (one for each model version) which will give you a list of all the variables in an RPN file. The first description is always correct, but rarely contains the unit. For physics fields take the second description with a grain of salt. If it matches the first one, it's unit is correct as well. If it doesn't match, ignore it. However, to all units of accumulators a 'per second' needs to get added, for example:
precipitation is not in [m] but in [m/s]
radiation fields are not in [J/m2] but in [J/s/m2] so in [W/m2]
Usage:
voir -> list the content of an RPN file
'voir' is the most simple way to list the content of an RPN file (also xvoir)
It prints the headers for all of the records of the given RPN standard file:
Without the -style option the date displayed is for raw GEM model output the start date of the simulation, for averages/accumulators/minimum/maximum fields it is the start of the avg/accum/min/max period and levels are coded.
With the -style option the date displayed is for raw GEM model output the validity (actual) date of the fields, for averages/accumulators/minimum/maximum fields it is the end of the avg/accum/min/max period and levels are decoded.
-style key | levels | date | |
---|---|---|---|
instantaneous fields | avg/accum/min/max | ||
without -style | coded | start date of simulation | start date of avg/accum/min/max period |
with -style | decoded | actual date | end date of avg/accum/min/max period |
You can find more information about the output of 'voir' in this PDF: voir_output_explanation.pdf
You can also use the following aliases:
alias vo='voir -iment'
alias vos='voir -style -iment'
Since the output of 'voir' can be rather long I usually pipe it into less, for example:
vos filename | less
xrec -> graphical visualization tool
xrec is a X-window visualization program used to display 2D (horizontal or vertical) meteorological fields stored in FST files. Usage:
or
You can also use the following alias:
alias xr='xrec -imflds'
The above command will open 3 windows (so make sure you connected with 'ssh -YC ...'!!!):
- xrecControlPanel : Here you make general selections about what you want to display
- xrecRecordSelector : Here all the records you can see with 'voir' (see above) are listed. You can double click on them to be displayed.
- xrecDisplay : Results will be displayed in this window. You might want to change its size to match your domain
Do not trust the wind direction for any other vectors but 'UU' and 'VV'!!!
Usage
Among other things, xrec can get used to :
- display 2-D fields
- vertical cross sections and profiles
- time series for selected points.
- time animation
Check the full documentation (xrec-documentation) about how to do this or just click on the buttons and look at the menus.
Dictionary
xrec is using a dictionary to define the title, unit (including a multiplication factor) and contour intervals for each variable.
The name of the dictionary is:
~/.recrc
It needs to have exactly this name and it needs to be in your home. To have a dictionary to start from you are welcome to take a copy of mine ( cp ~winger/.recrc ~/. ) or ask a college if you can copy hers/his.
You can find more information about this dictionary by clicking on the following link: xrec-dictionary
Remove default contour lines
The first time you open any field with xrec you will see there are "Colors" and "Contours". If you do not want to have "Contours" by default (you can always add or remove them later) do the following:
- Open xrec with any file
- Click on "Display" and unselect "Contours"
- Click on "File" and on "Quit and save configuration"
- Open xrec with any file again
- Select any file. If the contours are gone, you did the above correctly. If the background is now gray and you would like to have a white one again...
- Click on "Options" and "Contours" - a new window will open. In that window...
- Click on "Background" and select "Color:" "white"
- Click on "File" and on "Quit and save configuration"
The next time you open xrec the default contour lines should be back and the background should be white again.
High resolution geography
If you want to use high resolution geography on our UQAM servers you need to export:
Set min/max
You can use the Min-Max tool to fix a minimum and maximum value for a specific field.
1) Open the Min-Max menu window by clicking on "Options" in the Control Panel and select "Min-Max / Missing values...".
2) In the "Code" box search for and select the field for which you want to fix min/max.
3) In the "Max" box replace "Auto" with your maximum value. Do the same for "Min".
4) Now click on "Fixed"
Full documentation
You can find the full xrec documentation by clicking on the following link: xrec-documentation
r.fstliste -> list the content of an RPN file
r.fstliste is available if you need a bit more information than what is provided by voir.
It can be used in batch mode to extract specific or all of the record parameters of a RPN standard file record(s) based on a criteria given through the selection keys:
r.fstliste -izfst RPN-file [-datev -vdatev -ip1 -typvar -nomvar -col ...]
Code/Decode levels
In RPN files levels are coded. The different types of levels are:
- 0 - height [m] (metres)
- 1 - sigma [sg] (0.0->1.0)
- 2 - pressure [mb] (millibars)
- 3 - arbitrary code
- 4 - height [M] (metres) with respect to ground level
- 5 - hybrid coordinates [hy] (0.0->1.0)
- 6 - theta [th]
r.ip1 (Shell)
r.ip1 is a shell command that can be used to code and decode RPN levels as coded in RPN standard files. For more information click on the following link (login/password are the usual 'science'):
Or simply execute the following command in a terminal:
r.ip1 -h
convertIp (Python)
'convertIp' is a function of the RPN Python module to code and decode RPN levels. To get access to it you first have to import the module, for example with:
import rpnpy.librmn.all as rmn
Once the module is loaded you can get a documentation of 'convertIp' with:
help(rmn.convertIp)
Example to convert pressure level to coded IP1:
import rpnpy.librmn.all as rmn # Module to read RPN files
p_lev = 500
ip1 = rmn.convertIp(rmn.CONVIP_ENCODE, p_lev, rmn.KIND_PRESSURE) # Code level 'p_lev' hPa into IP1 code
print ('Coded value for', p_lev, 'hPa is:', ip1)
Answer:
Coded value for 500 hPa is: 41394464
Example to decode IP1:
import rpnpy.librmn.all as rmn
# Module to read RPN files
# Read one record filename = '...'
# Name of RPN file to read varname = '...'
# Name of variable to read fid = rmn.fstopenall(filename,rmn.FST_RO)
# Open the file rec = rmn.fstlir(fid,nomvar=varname)
# Read the full record of variable 'varname' ip1 = rec['ip1']
# Read the IP1 code of the level
(val, kind) = rmn.convertIp(rmn.CONVIP_DECODE, ip1)
# Decode IP1 into level value and unitprint ('Decoded value of', ip1, 'is:', val, rmn.kindToString(kind))
rmn.fstcloseall(fid)
# Close the RPN file
Answer:
Decoded value of 41394464 is: 500.0 mb
convip (Fortran)
'convip' is a Fortran function to code and decode RPN levels. 'convip' is part of the RPN Fortran library. To get access to it you need to load the library 'librmn.a' when creating the executable. You can do that by added the key '-lrmn' to the compile command. For example:
s.f90 source_code.f90 -o executable.Abs -lrmn
For more information click on the following link (login/password are the usual 'science'):
Code/Decode time stamps
Dates and times are coded in RPN files. This also means that it is not possible to have any times in RPN files:
- Between 1900 01 01 and 1980 01 01 the time resolution is 1 hour.
- Between 1980 01 01 and 2236 01 01 the time resolution is 5 seconds.
- Outside these time periods the time resolution is 3 hours.
r.date (Shell)
'r.date' is a shell command that can be used to code and decode RPN dates/times as coded in RPN standard files. For more information click on the following link (login/password are the usual 'science'):
Or simply execute the following command in a terminal:
r.date -h
newdate (Python)
'newdate' is a Python function to code and decode RPN dates/times. To get access to it you first have to import the module, for example with:
import rpnpy.librmn.all as rmn
Once the module is loaded you can get a documentation of 'newdate' with:
help(rmn.newdate)
Example to convert date & time to coded timestamp:
date = 20241106 # YYYYMMDD
time = 13300000 # hhmmsshh => 13h30
timestamp = rmn.newdate(rmn.NEWDATE_PRINT2STAMP, date, time) # Code date and time into timestamp
print ('Coded value for date', date, 'at', time, 'is:', timestamp)
Answer:
Coded value of date 20241106 at 13300000 is: 477041750
Example to convert coded timestamp to date & time:
# Read one record
filename = '...' # Name of RPN file to read
varname = '...' # Name of variable to read
fid = rmn.fstopenall(filename,rmn.FST_RO) # Open the file
rec = rmn.fstlir(fid,nomvar=varname) # Read the full record of variable 'varname'
timestamp = rec['datev'] # Read the code of the validity date
(date, time) = rmn.newdate(rmn.NEWDATE_STAMP2PRINT, timestamp) # Decodes timestamp to date (YYYYMMDD) and time (hhmmsshh)
print ('Decoded value of timestamp', timestamp, 'is:', date, 'at', time)
rmn.fstcloseall(fid) # Close the RPN file
Answer:
Decoded value of timestamp 477041750 is: 20241106 at 13300000
newdate (Fortran)
'newdate' is a Fortran function to code and decode RPN dates/times. 'newdate' is part of the RPN Fortran library. To get access to it you need to load the library 'librmn.a' when creating the executable. You can do that by added the key '-lrmn' to the compile command. For example:
s.f90 source_code.f90 -o executable.Abs -lrmn
For more information click on the following link (login/password are the usual 'science'):
Manipulating RPN files
The two main tools to modify FST files are called editfst and r.diag. And of course, you can also write your own C or Fortran programs.
editfst
editfst (also called editfst2000) is a utility used for editing and copying records from RPN standard files into a new or existing RPN standard file. It can do a straight (complete) copy of the input file or it can copy records selectively as indicated from the standard input or from a file of directives named in the "-i" key. One can select (via desire), exclude (via exclure) or zap records by their:
- TYPVAR: type of variable; 'A'(analysis),'P'(forecast),'X'
- NOMVAR: 4 caracter variable name
- ETIKET: 12 caracter label recorded from the original model run
- DATE: CMC date timestamp
- IP1: the level of the field (in pressure,eta,meters,other)
- IP2: the hour of the forecast (rounded off if not exact)
- IP3: value of 0 unless otherwise modified
Example:
You could use it to select i.e. temperature (TT) and winds (UU,VV) at 500 hPa:
echo "desire (-1,['UU','VV'], -1 ,-1 ,[500.,MBAR])" | editfst -s input_file -d output_file
Note:
editfst appends to its output file. This means that if editfst is used more than once writing into the same output file, the output of the second set of commands will be appended to the output of the first set.
r.diag
r.diag (also called r.diag2000) is a toolkit used to manipulate FST or CCCma files. As editfst, it can be used to select or exclude specific records but in addition to this it also contains a large set of functions which can be used to manipulate the files.
'r.diag' is always used in the way:
r.diag function input_file(s) [ output_file(s) key(s) ]
To get information about a specific function simply type:
r.diag function
Notes:
- r.diag generally overwrites its output files. There are a few functions, eg. select, which have an option to append to the output file instead of overwriting it.
- By default, output files will be written under ${TMPDIR}!!!
If there is no path written in front of the output file name and if the output file does not yet exist in the local directory, it will be written in one of the ${BIG_TMPDIR} or ${TMPDIR} directories. Also, if the input file is not in the local directory, r.diag will look for it in the ${BIG_TMPDIR} and ${TMPDIR} directories.
So if you want to have the output file in a certain directory you have to specify said directory name before the file name. Fir the current directory you can just preceed the file name by './'.
Examples:
'select' - Select fields according to name and level
For example, temperature (TT) and winds (UU,VV) at 500 hPa can get selected with 'select' with the following command:
r.diag select input_file output_file_tt output_file_uu output_file_vv -name TT UU VV -lv1 500
'gsapl' - Interpolate model to pressure levels
Fields on model levels can get interpolated to a set of pressure levels with 'gsapl'.
Input files needed
Files can have any names, but to facilitate the explanation the following names are chosen for the example below:
- mod_lev : input file on model levels
- pres : input file containing surface pressure for the same time steps which are in 'mod_lev'
- niveau.dat : file containing list if pressure levels to interpolate to. See more below!!!
Also just to facilitate the example below, the following parameters should get set:
- plv : parameter set to number of pressure levels in 'niveau.dat'
- kind : parameter set to type of model levels in the input file. Depending on the model version with which the output got produced set:
kind=GEM2 : sigma model level ('hyb' in gem_settings.nml) starting NOT with 0.000
kind=GEM3 : hybrid model level ('hyb' in gem_settings.nml) starting with 0.000
kind=GEM4 : GEM4 staggered fields with '!!' record
kind=GEM5 : GEM5 staggered fields with '!!' record
File 'niveau.dat'
You need to create this file which needs to contain the list of all the pressure levels you want to interpolate to. This file needs to have the following format:
- max 16 levels per row
- 5 digits per level
- right bound
- all integer numbers!!!
Here are some example to put pressure levels smaller than 1:
'-1100' stands for '0.1'
'-1200' stands for '0.2'
:
' -100' is the same as '1'
For example:
# Create file 'niveau.dat' (instead of with 'cat' you can also create it with an editor)
cat > niveau.dat <<+ -1100-1200-1400-1600 -100 2 4 6 10 20 30 50 70 100 150 200
250 300 400 500 600 700 850 925 1000
+
kind=GEM5 # When output is from GEM5
plv=25 # Number of pressure levels in niveau.dat
# As input, 'gsapl' needs the log of the surface pressure. So you need to calculate this first, for example:
# Create the log of the surface pressure
r.diag loge pres log_pres
# Change variable name to 'LP' (required name by 'gsapl')
r.diag newnam log_pres lp -name LP
# Now the file 'lp' contains the log of the surface pressure and the variable has the name 'LP'.
# IF(!) there is more than 1 variable on model levels in your input file, you need to create a set of surface pressure for each variable!!!
# To do this, do the following. Otherwise skip the next four lines of code!
# Create a file 'mlp' containing as many records of the same ln(P0) as there are variables in the input file 'mod_lev'.
# And of course there must be one "set" of ln(P0) for each time step.
r.diag select cubet univeau -lv1 1.0
r.diag xlin univeau un -a 0 -b 1 -name LP
r.diag gmlt un lp mlp
mv mlp lp
# Interpolating all fields on model levels of the input file 'mod_lev' to pressure levels
# Make sure the input file 'mode_lev' contains the level record: 'HY' for GEM3 and '!!' for GEM4 and up.
r.diag gsapl mod_lev lp pres_lev -plv ${plv} -kind ${kind} < niveau.dat
The file 'pres_lev' should now contain the same fields and timesteps as the input file 'mod_lev' but on the pressure level set in 'niveau.dat'.
Writing your own C and Fortran programs
There is a number of functions which can be called in a C or Fortran program to read and write FST files, and they are documented in the RPN Standard File Functions. These functions are part of the library 'librmn'.
There are functions to:
open and close an FST file
search for a record
read a record
read the field descriptor of a record
- write a recordcdf2rpn
Have a look here on how to compile.
Conversion: RPN ↔ NetCDF format
The tool 'cdf2rpn' which can get used to convert RPN to NetCDF format and vice versa. You can get the documentation of cdf2rpn by simply executing the command in a terminal:
The general usage is the following:
Convert from NetCDF to RPN format:
cdf2rpn -rpn RPN_standard_file -cdf NetCDF_file -attr attribut_list
Convert from RPN to NetCDF format - one needs to add the key '-dir':
cdf2rpn -rpn RPN_standard_file -cdf NetCDF_file -attr attribut_list -dir
Notes:
- If you do not add a path to the output file name the file will end up in $TMPDIR!!!
- All fields in a file must have the same "type" of level.
You cannot have i.e. sigma, hybrid, arbitrary and pressure levels in the same file!!! If your RPN input file contains different types of fields you can select the fields your want to convert with 'r.diag' and/or 'editfst' - see above. - If converting RPN to NetCDF make sure your RPN file contains all necessary grid descriptors ('^^', '>>', '^>', 'HY','!!')
- All variables you want to convert need to be in the attribute list.
You can find an example file on the servers under: ~winger/Scripts/NetCDF_converter/attribut_netcdf.dat
Make sure that the factors are correct!!! If one or more of your variables are not in this file, just take a copy and add the missing variables.
2 commentaires
Anonyme dit :
Once the module is loaded you can get a documentation of 'convertIp' with:
help(rmn.newdate)
In the newdate (Python) section, it is 'newdate' and not 'convertIp'
Winger, Katja dit :
Thank you! Changed!
Ajouter un commentaire