GrADS Tutorial
April 2009
An Introduction to GrADS Luke Garde (room 330)
[email protected]
6th April 2009
Sample scripts and data files 1. 2. 3. 4. 5. 6.
Run Cygwin.exe (should be an icon on the Desktop) > cd /usr/X11R6/bi /usr/X11R6/bin n > startx > ssh –X cove (you should be in your home directory i.e. /home/user) > mkdir grads_tute > cd grads_tute 7. > cp /home/lgarde/grads_tutorial/grads_pack.tar.gz . 8. > gunzip grads_ pack.tar.gz 9. > tar xvf grads_ pack.tar This will copy a collection of files to your home directory and extract them. Between tasks you can simply exit GrADS (quit) or reinitialize GrADS (reinit).
Note: commands written with “>” are ar e UNIX commands, whilst “ga ->” are GrADS commands All information will be available for download via my Earth Science’s personal page:
http://www.earthsci.unimelb.edu.au/~lgarde/grads/ Task 1
Command-Line
1. Execute GrADS > grads 2. Open mslp.nc ga-> sdfopen mslp.nc 3. Display mean sea level pressure (notice contours are in Pa) ga-> d mslp 4. Clear display ga-> c 5. Display mean sea level pressure again, but this time in hPa ga-> d mslp/100 6. Set plot boundaries to be zoomed to the Australian region and display mean sea level pressure (in hPa) again. ga-> set lat -55 -5 ga-> set lon 110 170 ga-> d mslp/100 Notice plot draws onto of previous output
1
GrADS Tutorial
April 2009
7. Clear display 8. Display mean sea level pressure again d mslp Mean sea level pressure will be drawn in Pa This is because we have not defined (allocated to memory) the converted mslp hPa field. 9. Define new mslp (hPa) field ga-> mslp2 = mslp/100 This will memory allocate the new variable mslp2 Note: Variable names must not include special characters (@_#$% etc) 10. Clear and display new mslp2 variable
Task 2
Command-Line
1. Open mslp.nc 2. Define mslp in hPa Define variable as mslp2 By quitting GrADS or using the reinit command, all previously defined variables are cleared, hence mslp2 needs to be defined again 3. Average mslp (in hPa) covering January – December 1979 and display avemslp ga-> avemslp = ave(mslp2,t=1,t=12) ga-> d avemslp 4. Create a title for image ga-> draw title Average MSLP Jan – Dec 1979 5. Create a hardcopy print using printim command ga-> printim ave1979mslp.gif white 6. Create a hardcopy print using enable print command ga-> enable print ave1979mslp2.gmf ga-> print ga-> disable print 7. Quit GrADS 8. On Unix command line convert outputted ave1979mslp2.gmf to ave1979mslp2.eps > gxeps -acR -i ave1979mslp2.gmf -o ave1979mslp2.eps > rm ave1979mslp2.gmf (optional, simply removes .gmf file) 9. Now display ave1979mslp.gif (image created using printim command) > display ave1979mslp.gif 10. Now display ave1979mslp2.eps (image created using gxeps) > display ave1979mslp2.eps You might notice a difference in image quality between the two images. An eps image is a Postscript image. This vector image can be zoomed to high levels with no pixel distortion
Note that display command will not work under atlas!
2
GrADS Tutorial
Task 3
April 2009
Opening two files from Command-Line
1. Execute GrADS ga-> sdfopen u10m.nc ga-> sdfopen v10m.nc 2. Open both u10m.nc and v10m.nc 3. Query opened files ga-> q files This will print to screen each opened file 4. Set plot boundaries to be zoomed to the Australian region Lat = -55 -5 N, Lon = 110 170 E 5. Display uwnd contours from 1 st file ga-> d uwnd 6. Clear 7. Display vwnd contours from 2 nd file ga-> d vwnd.2 8. Clear display and draw wind vectors and colour magnitudes ga-> set gxout vector ga-> d uwnd;vwnd.2;mag(uwnd,vwnd.2) 9. Clear and draw wind barbs ga-> set gxout barb ga-> d uwnd;vwnd.2 10. Clear and draw wind streamlines ga-> set gxout stream ga-> d uwnd;vwnd.2;mag(uwnd,vwnd.2)
Task 4
Scripting (.gs)
1. Open your favourite text editor (Nedit, Pico, vim) > nedit 2. Create a script that plots wind streamlines. Include in your script code to prompt the user to define plot boundaries. Remember that when scripting, GrADS commands need to be enclosed by ' ' Notice that some commands are not enclosed by ' ' (just to make it confusing!) In addition, you can make comments in your script using * which will not affect the progress of the running script Enter script as follows in text editor (copy - paste): 'reinit' 'sdfopen u10m.nc' 'sdfopen v10m.nc' * Remove GrADS logo from plot 'set grads off' 'set map 1 1 3' 'enable print 10m_stream.gmf' * User prompt section prompt 'Enter min latitude: ' pull minlat prompt 'Enter max latitude: ' pull maxlat
3
GrADS Tutorial
3. 4.
April 2009
prompt 'Enter min longitude: ' pull minlon prompt 'Enter max longitude: ' pull maxlon 'set lat 'minlat%' '%maxlat 'set lon 'minlon%' '%maxlon * Plotting section 'set gxout stream' 'd uwnd;vwnd.2;mag(uwnd,vwnd.2)' 'draw title 10m Streamlines' 'print' 'disable print' * External gxeps plotting '!gxeps -acR -i 10m_stream.gmf -o 10m_stream.eps' '!rm 10m_stream.gmf' 'quit' Save GrADS script as 10m_stream.gs Run 10m_stream.gs > grads -lc 10m_stream.gs (run from Unix terminal) or ga-> run 10m_stream.gs (run from inside GrADS) Open 10m_stream.eps plot as before Plotting via script (.gs) enables you to repeat your data run as many times as needed, without having to type commands on the command line (ga->). This also allows you to record how you plot your images throughout the year. Repeat 4 – 5 for different latitude and longitude ranges
5. 6.
7.
Task 5
Colour Maps (.gs)
In this task we are going to construct a colour map instruction script. This script will define the colours and levels that GrADS will use to create your plots 1. Open text editor 2. Open sample_col.gs 3. You will see the following: 'set rgb 16 255 255 255' 'set rgb 17 31 61 250' 'set rgb 18 0 199 199' 'set rgb 19 0 220 0' 'set rgb 20 161 230 51' 'set rgb 21 230 220 51' 'set rgb 22 240 130 41' 'set rgb 23 240 0 0' 'set clevs 5 10 15 20 25 30 35' 'set ccols 16 17 18 19 20 21 22 23'
cbarn.gs
22 23
21
19
20
18
17 16
4
GrADS Tutorial
April 2009
There are a couple of things to notice: a) Set colours are defined starting at 16 through to 23. This is the “identification number” of that colour which you call in t he set ccols command. Note: numbers 0 – 15 are GrADS default colours b) Each colour is mapped via RGB triplets, where i. White: 255 255 255 ii. Black: 0 0 0 iii. Red: 255 0 0 iv. Green: 0 255 0 v. Blue: 0 0 255 vi. Combinations between these triplets are therefore mixed colours c) set clevs defines contour levels of 5 – 35 in steps of 5. Hint: If using the cbarn.gs colour bar the number of arguments in set clevs must be one less than the number for set ccols. This will insure that the top triangle will have an assigned colour and not default grey. 4. When happy with set colours and levels, save file. 5. Create a GrADS script called mslp_u_v.gs and write the following 6. Open gl2008011800.nc 7. Enable the print as mslp_u_v.gmf 8. set GrADS logo off 9. set the boundaries to the Aust/South Pac [eg: -55 -5 N, 110 200 E] 10. set the level of the plot to 250hPa (set command as in lon and lat) 11. set the graphic output to shaded 12. call new colour map 'run sample_col.gs' 13. display the magnitude of both u and v winds notice that this time u and v parameters are called zonal_wnd and merid_wnd respectively. You would know this if you typed “q file” on the GrADS commandline (ga->) 14. call cbarn.gs to draw colour bar 'run cbarn.gs' 15. set the graphics to contours 16. display mslp Remember you must set the level to 1000hPa before displaying mslp as this is the level in which mslp is defined. If you don't the message “Entire Grid Undefined” will be displayed in your plot 17. Set the title of the plot to: MSLP & 250hPa winds (magnitude) 18. Create hardcopy file (remember print, disable print) 19. Set up external image format conversion '!gxeps -acR -i mslp_u_v.gmf -o mslp_u_v.eps' '!convert +antialias -density 300 -geometry 800x600 mslp_u_v.eps mslp_u_v.gif' 20. Save mslp_u_v.gs and exit text editor 21. Run GrADS script > grads -lc mslp_u_v.gs or ga-> run mslp_u_v.gs 22. Quit GrADS 23. Display mslp_u_v.gif which has been created in your directory > display mslp_u_v.gif
5
GrADS Tutorial
Task 6
April 2009
Command-Line and Data Average Computation
1. Open GrADS 2. Open gl2008011800.nc 3. Set latitude limits for averaging ga-> set lat -40 -20 (set lat -90 90 would be globally) 4. Set fixed longitude ga-> set lon 90 (for display purposes this needs to be fixed) 5. Set level from 1000 – 70 hPa ga-> set lev 1000 70 6. Define azonal as the global average of zonal wind between -40 -20 ga-> azonal = ave(zonal_wnd,lon=0,lon=360) or ga-> azonal = ave(zonal_wnd,lon=0,lon=360,-b) 7. Display new azonal variable. This will show a cross section view of the average global zonal wind 8. clear and repeat 3 - 6 a couple of times, experimenting with latitude and longitude boundaries and also the -b flagging option 9. Clear and now set plot boundaries are follows: Level = 1000 Latitude = -55 -5 N Longitude = 110 170 E 10. Now calculate and area average of mslp between about lat and lon boundaries define aamslp = aave(mslp,lon=110,lon=170,lat=-55,lat=-5) 11. Display aamslp Notice numerical output to GrADS terminal (“Result value = 1008.58” ) 12. Repeat 9 - 11 a couple of times, experimenting with latitude and longitude boundaries
Task 7
Calling files via Data Descriptor files (.ctl)
1. Open sample_desc.ctl in text editor 2. You will see the following, however the ___ spaces require your input. DSET gl%y4%m2%d2%h2.nc DTYPE netcdf OPTIONS template yrev TITLE netCDF file UNDEF -9.e36 XDEF ___ LINEAR ___ 1.5 YDEF ___ LINEAR ___ ___ ZDEF 31 LEVELS ___ 995 ___ 985 ___ ___ ___ ___ TDEF 2 LINEAR 00:00z17JAN2008 1dy VARS ___ sfc_pres=>ps 1 t,y,x Surface pressure sfc_temp=>pt 1 t,y,x Surface temperature mslp 1 t,y,x a geop_ht=>h 31 t,z,y,x a air_temp=>t 31 t,z,y,x a mix_rto=>qv 31 t,z,y,x a zonal_wnd=>u 31 t,z,y,x a merid_wnd=>v 31 t,z,y,x a ENDVARS 6
GrADS Tutorial
April 2009
3. Using ncdump find out the dimensions/resolution of the file > ncdump gl2008011800.nc | more Space bar or Enter moves through the output and q quits Correct answers shown in PowerPoint presentation Data Descriptor file must be 100% correct or files will be opened incorrectly 4. Save sample_desc.ctl and exit text editor 5. Open GrADS 6. Open Data Descriptor File open sample_desc.ctl 7. Display mslp for t = 1 and t = 2 ga-> set t 1 ga-> d mslp ga-> set t 2 ga-> d mslp 8. Now you have a way of opening multiple netCDF files. This is very useful as data provided from the BoM is normally in single time steps. Data Descriptor files are unique to the data you are plotting.
Task 8
Scripting (.gs) and Flow Control
You may recognise this event as the Sydney to Hobart 1998 Cyclone 1. 2. 3. 4.
Run anim_ex.gs and follow on screen instructions Rerun anim_ex.gs for different time steps Open anim_ex.gs in text editor Try to understand how the script operates and “controls the flow” 'reinit' 'sdfopen syd_hobart.nc' 'set t 1' say ' ' prompt 'Enter number of time steps (max 14): ' pull tstep say ' ' say 'Time steps set to: 'tstep 'set dbuff on' i=1 while (i <= tstep) 'set t 'i say 'time set to: 'i 'set mpdset hires' 'set map 1 1 3' 'set ccolor 1' 'd msl/100' 'set gxout stream' 'd p10u;p10v;mag(p10u,p10v)' 'draw title MSLP & 10m Wind ['i' of 'tstep']' 'swap' * Pause screen output (in sec) 7
GrADS Tutorial
April 2009
'!sleep 1' i=i+1 endwhile say 'DONE: Loop completed!' say 'Type: run anim_ex.gs to run again'
Task 9
Climatology and NCEP Website Composite
1. Run clim.gs 'reinit' 'sdfopen mslp.nc' 'set grads off' 'enable print mslp_map.gmf' 'set lon 0 360' 'set lat -55 60' 'set lev 1000' 'set t 1 12' 'define spclim = ave(mslp,t+0,t=348,12)' 'modify spclim seasonal' 'set t 1' 'set gxout shaded' 'd spclim(t=1)/100' 'run cbarn.gs' 'draw title Mean Jan MSLP' 'print' 'disable print' '!gxeps -acR -i mslp_map.gmf -o mslp_map.eps' '!convert +antialias -density 150 -geometry 1024x768 mslp_map.eps mslp_map.gif' 'quit' 2. Display mslp_map.gif 3. Open Firefox Web Browser 4. Go to www.cdc.noaa.gov/cgi-bin/Composites/printpage.pl 5. Select Sea Level Pressure 6. Select January for both starting and ending months 7. Enter range of years: 1979 – 2007 8. Select shading with contours 9. Change plot size to 150% 10. Select custom for map projection 11. Set latitude: -55 60 N 12. Set longitude: 0 360 E 13. Change projection to Cylindrical Equidistant 14. Create Plot 15. Compare Webpage output to GrADS output mslp_map.gif 16. You could then download this data by clicking “ Get a copy of the netcdf data file used for the plot” and plot this data again using your own colour maps and in eps format.
8
GrADS Tutorial
Task 10
April 2009
Plotting from Cyclone Tracking Scheme
1. Change to task10 directory > cd task10 2. Open trackplot.f in a text editor and edit “User defined criteria for filtering tracks” Change minlength to 2 days Pick latitude range Save and exit 3. As you have changed the FORTRAN code, you need to compile the routine > g77 -o trackplot trackplot.f 4. Now execute the FORTRAN routine on your trkdat file > trackplot trkdat.test 5. This would have created a file called tracks in your working directory > ls -t this will show files in chronological order (last to first) 6. Open trackplot.gs.template in a text editor 7. Here you can change Mark, Line and Projection options 8. Run trackplot.gs.template > grads -lc trackplot.gs.template 9. Display plot output > display plot.png
9
