Doing GRIB in GrADS Mike Fiorino Program for Climate Model Diagnosis and Intercomparison Lawrence Livermore National Laboratory University of California P.O. Box 808 L-264 Livermore, CA 94551 510-423-8505 (voice) 510-422-7675(fax) 21 February,1996 1) INTRODUCTION One of the most powerful features of GrADS is its ability to work DIRECTLY with GRIB data. Unfortunately, "doing GRIB in GrADS" is not simple; it requires an understanding of GRIB and how GrADS works this type of data. This note attempts to provide the required understanding to use GRIB data in GrADS. Your comments, preferably, by email are requested. First, let's start with some basic resources: GrADS - http://grads/iges.org/grads/head.html - ftp://sprite.llnl.gov/grads GrADS doc - ftp://sprite.llnl.gov/grads/doc/gadoc151.* Note: Appendix B in gadoc151.* is another starting point for the GRIB-GrADS interface GRIB - ftp://ncardata.ucar.edu/docs/grb/guide.txt This is the "Gospel According to John" (Stackpole that is) GRIB-GrADS - http://grads.iges.org/grads/ncepdata.htm http://hera.gsfc.nasa.gov/grads_listserv/14Jun95.20920 2) WHAT IS GRIB? GRIB (GRIdded Binary) is a international, public, binary format for the efficient storage of meteorological/oceanographic variables. Typically, GRIB data consists of a sequence of 2-D (lon,lat) chunks of a, most generally, 4-D variable (e.g., u comp on the wind = f(lon,lat,level,time)). The sequence is often organized in files containing all variables at a particular time (i.e., 3-D (lon,lat,level) volume). A good example of this kind of organization is: ftp://nic.fb4.noaa.gov/avn/avn.00Z/gblav.T00Z.PGrbF00 which contains the analyses (F00) at 17 pressure levels from the 00Z AVN run of the NCEP MRF global model. 3) THE PROBLEM The problem for the user is how to "interface" or "display" GRIB file(s) to GrADS. The solution has two components. The first is to see what is in the GRIB file, i.e., variables, grids, times, etc and the second is to "map" the 2-D GRIB fields to higher dimensional structures available in GrADS. To make more concrete what is happening, I have set up a sample problem. See, ftp://sprite.llnl.gov/pub/fiorino/grads/grib where you will find two GRIB files (*.grb) from the NCEP reanalysis: grb2ctl.sh* ncep.reanl.mo.7901.grb ncep.reanl.mo.gmp ncep.reanl.mo.7902.grb ncep.reanl.mo.ctl the "7901" in the file name indicates the fields come from jan 1979 and "7902" feb 1979. The importance of the file name convention will become apparent... Also included is a helper unix shell script (grb2ctl.sh) and the "answer" ncep.reanl.mo.ctl and the associated "index" file ncep.reanl.mo.gmp which will be explained in a moment. 4) UP-FRONT LIMITATIONS There are some limitations on the kinds of GRIB data that can be interfaced/displayed in GrADS: a) lon,lat grids (NOT lat,lon) b) simple packing c) grid point data d) grids must be contiguous (no blocked or octet grids) Thus, "thinned" grids (non rectangular) and spectral coefficients are not supported. However, GRIB versions 1 AND 0 are supported, (GRIB 0 data must be filtered to GRIB 1 for wgrib (see ftp://wesley.ncep.noaa.gov/pub/grib0)). Further, it IS possible to display "preprojected" GRIB data (e.g., polar stereo fields, see http://grads.iges.org/grads/ncepdata.htm and http://grads.iges.org/grads/eta.ctl). 5) THE SOLUTION - PART 1 The first, and relatively straightforward, part of the solution is to find out what is in the GRIB file. I use two utilities: 1) gribscan (comes with GrADS); and 2) wgrib. While I contributed to Brian Doty's gribscan, my recommendation is that you go with wgrib written by Wesley Ebisuzaki of NCEP. The big virtue of wgrib is that the code is written in ANSI C and runs on all the major platforms from PC's to Cray. Although wgrib was designed to support the NCEP reanalysis project, it has been extended to handle other sources of GRIB data (e.g., ECMWF) and is my tool of choice. If I can't read the data in wgrib then I change the code and feed the changes to Wesley. See: ftp://wesley.ncep.noaa.gov/pub/wgrib and, ftp://nic.fb4.noa.gov/pub/info/reanl/wgrib Once you have wgrib running, wgrib ncep.reanl.mo.7901.grb yields, 1:0:d=79010100:UGRD:kpds5=33:kpds6=100:kpds7=850: TR=113:P1=0:P2=6:TimeU=1:850 mb:anl:ave@6hr:NAve=124 2:15852:d=79010100:UGRD:kpds5=33:kpds6=100:kpds7=500: TR=113:P1=0:P2=6:TimeU=1:500 mb:anl:ave@6hr:NAve=124 3:33018:d=79010100:UGRD:kpds5=33:kpds6=100:kpds7=200: TR=113:P1=0:P2=6:TimeU=1:200 mb:anl:ave@6hr:NAve=124 4:51498:d=79010100:VGRD:kpds5=34:kpds6=100:kpds7=850: TR=113:P1=0:P2=6:TimeU=1:850 mb:anl:ave@6hr:NAve=124 5:66036:d=79010100:VGRD:kpds5=34:kpds6=100:kpds7=500: TR=113:P1=0:P2=6:TimeU=1:500 mb:anl:ave@6hr:NAve=124 6:81888:d=79010100:VGRD:kpds5=34:kpds6=100:kpds7=200: TR=113:P1=0:P2=6:TimeU=1:200 mb:anl:ave@6hr:NAve=124 7:99054:d=79010100:PRES:kpds5=1:kpds6=102:kpds7=0: TR=113:P1=0:P2=6:TimeU=1:MSL:anl:ave@6hr:NAve=124 So we see 7 fields in the file valid at 00Z1jan1979 (d=79010100). for, wgrib ncep.reanl.mo.7902.grb we find, 1:0:d=79020100:UGRD:kpds5=33:kpds6=100:kpds7=850: TR=113:P1=0:P2=6:TimeU=1:850 mb:anl:ave@6hr:NAve=112 2:15852:d=79020100:UGRD:kpds5=33:kpds6=100:kpds7=500: TR=113:P1=0:P2=6:TimeU=1:500 mb:anl:ave@6hr:NAve=112 3:33018:d=79020100:UGRD:kpds5=33:kpds6=100:kpds7=200: TR=113:P1=0:P2=6:TimeU=1:200 mb:anl:ave@6hr:NAve=112 4:50184:d=79020100:VGRD:kpds5=34:kpds6=100:kpds7=850: TR=113:P1=0:P2=6:TimeU=1:850 mb:anl:ave@6hr:NAve=112 5:64722:d=79020100:VGRD:kpds5=34:kpds6=100:kpds7=500: TR=113:P1=0:P2=6:TimeU=1:500 mb:anl:ave@6hr:NAve=112 6:80574:d=79020100:VGRD:kpds5=34:kpds6=100:kpds7=200: TR=113:P1=0:P2=6:TimeU=1:200 mb:anl:ave@6hr:NAve=112 7:96426:d=79020100:PRES:kpds5=1:kpds6=102:kpds7=0: TR=113:P1=0:P2=6:TimeU=1:MSL:anl:ave@6hr:NAve=112 or the same fields as before, except they are valid at 00z1feb1979. To find out about the data grid use, wgrib -V ncep.reanl.mo.7901.grb and for the first record you will find: rec 1:pos 0:date 79020100 UGRD kpds5=33 kpds6=100 kpds7=850 levels=(3,82) grid=2 850 mb anl:ave@6hr: timerange 113 P1 0 P2 6 nx 144 ny 73 GDS grid 0 num_in_ave 112 missing 0 center 7 subcenter 0 process 80 latlon: lat 90.000000 to -90.000000 by 2.500000 long 0.000000 to -2.500000 by 2.500000, (144 x 73) scan 0 bdsgrid 1 min/max data -12.29 16.86 num bits 12 BDS_Ref -1229 DecScale 2 BinScale 0 The grid is the NCEP #2 grid (see the Gospel According to John) which is a global 2.5 deg grid. The grid has 144 points in x (lon) and 73 points in y (lat). The (1,1) point is located at 90N and 0E. 6) THE SOLUTION - PART 2 This is the hard part -- creating a relationship between the sequence of 2-D (lon,lat) GRIB fields in a file(s) and a GrADS, 4-D, external-to-the-data, spatio-temporal volume (lon,lat,level,time). In GrADS, the GRIB-to-4-D volume relationship is defined by the "data descriptor" or ".ctl" file. The actual relationship is created using the GrADS utility "gribmap" which generates an "index" or "map" between GrADS variables in the .ctl file and the GRIB data. Before describing the details of gribmap, first consider the unix shell script, grb2ctl.sh, originally written by Wesley which I have adapted to help me out. This script uses wgrib to first create a listing of fields in the GRIB data file and then parses the listing for times, variables and levels for the .ctl file. See, ftp://sprite.llnl.gov/grads/grib/grb2ctl.sh Wesley's version is ftp://wesley.ncep.noaa.gov/pub/wgrib.scripts/grib2ctl In our example, grb2ctl.sh ncep.reanl.mo.7901.grb yields to standard output, dset ^ncep.reanl.mo.7901.grb dtype grib options yrev index ^ncep.reanl.mo.7901.grb.gmp undef -9.99E+33 title ncep.reanl.mo.7901.grb xdef 144 linear 0 2.5 ydef 73 linear -90 2.5 zdef 3 levels 850 500 200 tdef 1 linear 00Z01jan79 1mo vars 3 PRES 0 1 ,102,0 Pressure [Pa] UGRD 3 33 ,100 u wind [m/s] VGRD 3 34 ,100 v wind [m/s] endvars How did grb2ctl.sh do this? First, only ONE grid was found in the GRIB data file (defined by dset) and the script had the grid geometry built in. Second, variables with multiply levels (3-D or lon,lat,level) where DEFINED to have a "level indicator" of 100 (more below). Third, there was only ONE time in the file and the script SET the time increment to 1mo. These conditions will often not occur. Thus, the output from grb2ctl.sh will likely have to be "tweeked" by the good 'ol trial and error method. In some cases the .ctl file may have to built "by hand" using the output from wgrib, so you'll have to know more about GRIB and how gribmap works in most situations. The key ingredients in the .ctl file are: a) grid geometry (xdef,ydef) b) starting time and time increment (tdef) c) variables and "units" parameter d) variable type - "level" or 3-D (zdef) or "surface" (2-D) The units parameter specifies the GRIB parameters of the variable in the .ctl to be used by gribmap for match GrADS variables to the fields in the GRIB files. This parameter consists of up to four, comma-delimited numbers: VV,LTYPE,(LEVEL),(TRI) where, VV - (Required) The GRIB parameter number (33 = u comp of the wind, table 2 in John:section 1 page 27, i.e., GRIB Edition 1 (FM 92)); [WNE: this is the kpds5 value from "wgrib (file)"] LTYPE - (Required) The level type indicator (100 = pressure level, in Table 3 in John:section 1 Page 33) [WNE: this is the kpds6 value from "wgrib (file)"] LEVEL - (optional) The value of the LTYPE (LTYPE 102 is mean sea level so LEVEL is 0 for where level is located (1,102,0, 0 is AT mean sea level) [WNE: this is the kpds7 value from "wgrib (file)"] TRI - (optional) The "time range indicator" for special applications (Table 5 in John:section 1 Page 37). [WNE: this is the TR value from "wgrib (file)"] Coming up with the units parameter, the grid geometry and the times is the trick. Fortunately, gribmap can tell us how well the .ctl mapped the GRIB data to the higher dimensional, GrADS data view. And, more importantly, the gribmap process does NOT depend on how the data are actually ordered in the GRIB file, in either level or variable. Let's redirect output from grb2ctl.sh to a .ctl file, e.g., grb2ctl.sh on ncep.reanl.mo.7901.grb > ncep.reanl.mo.ctl and then run gribmap. To reiterate, the gribmap utility compares each field in the GRIB file to each variable, at each level and for all times in the .ctl file and creates an index file telling GrADS WHERE the fields are (or are not) located in the GRIB data. With the verbose option on, gribmap -v -i ncep.reanl.mo.ctl we get, Scanning binary GRIB file(s): ncep.reanl.mo.7901.grb !!!!! MATCH: 1 15852 2 1 0 33 100 850 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 !!!!! MATCH: 2 33018 2 1 0 33 100 500 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 !!!!! MATCH: 3 51498 2 1 0 33 100 200 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 !!!!! MATCH: 4 66036 2 1 0 34 100 850 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 !!!!! MATCH: 5 81888 2 1 0 34 100 500 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 !!!!! MATCH: 6 99054 2 1 0 34 100 200 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 !!!!! MATCH: 7 116220 2 1 0 1 102 0 btim: 1979010100:00 tau: 0 dtim: 1979010100:00 Reached EOF We have succeeded!!! Each 2-D in the GRIB file has been mapped to a variable in ncep.reanl.mo.ctl. In the case of UGRD, a 3-D (lon,lat,level) variable which we can be sliced in the vertical with GrADS. However, failure to match will NOT stop GrADS from "working." If the data was NOT there, GrADS will return a grid with "undefined" values on display and this state can actually be tested... The "tweeking" is done by adjusting the .ctl file until we get a "!!!!! MATCH" for each GRIB field in the data file(s). I have added a number of options that finely control the mapping process in gribmap for NCEP. See the GrADS document for details (ftp://sprite.llnl.gov/grads/doc/gadoc151.*). Finally, let's adjust the ncep.reanl.mo.ctl file to take advantage of the file naming convention: dset ^ncep.reanl.mo.%y2%m2.grb dtype grib options yrev template index ^ncep.reanl.mo.gmp undef -9.99E+33 title ncep.reanl.mo.7901.grb xdef 144 linear 0 2.5 ydef 73 linear -90 2.5 zdef 3 levels 850 500 200 tdef 2 linear 00Z01jan79 1mo vars 3 PRES 0 1 ,102,0 Pressure [Pa] UGRD 3 33 ,100 u wind [m/s] VGRD 3 34 ,100 v wind [m/s] endvars I have changed the number of times to two and have used the template option. This tells GrADS to locate data in TIME by the file name (dset ^ncep.reanl.mo.%y2%m2.grb). Thus, I have two (or more) data files, but only ONE .ctl file and I can now work with the 2-D GRIB data as if they were 4-D (lon,lat,lev,time) in GrADS. I also changed the name of the index file to be reflective of the now 4-D data structure. To summarize the process: 1) use wgrib (or gribscan) to see if the data can be worked in GrADS; 2) use grb2ctl.sh or the output from wgrib to construct a .ctl file; 3) run gribmap in verbose mode (-v) to relate the GRIB data to the 4-D structure in the .ctl file, and to see how well the map worked; and 4) repeat steps 2) and 3) until you get "!!!! MATCH" from gribmap. 7) CONCLUSIONS There's no doubt about it, "Doing GRIB in GrADS" is not very straightforward, but the benefits are, in my opinion, immense for two reasons. First, we have avoided conversion of the GRIB to a form which supports higher dimensional data (e.g., netCDF). We've saved disk space and have minimized potential technical errors (every time you touch the data you have an opportunity to screw it up). Second, from a GrADS performance standpoint, GRIB is nearly as fast as other binary formats -- the cost in decompression on the fly is compensated by reduced I/O. In the end, GRIB-to-GrADS interface gives us the advantages of GRIB (efficient storage, self description and an open, international format) while overcoming the disadvantages of GRIB (2-D data and no means to organize to a higher dimension) via the GrADS 4-D data model. We get the best of both worlds, but only if we can make the .ctl file. Hopefully this document will help you do this.