-
Notifications
You must be signed in to change notification settings - Fork 35
GC Tutorials
Quick intro to:
- converting to CF compliant NetCDF
- viewing file as CDL (like ncdump)
- Plot tileplot/map figures
- Plot vertical profiles.
- Extract monitor locations.
- Perturb emissions in an existing file.
- Plot a pressure average
Assuming:
- /path/to/ctm.bpch is the path to a bpch file.
- /path/for/output/ctm.nc is the path you want to output to.
- tracerinfo.dat and diaginfo.dat are either in your working directory or colocated with the bpch input file.
Example creates a file that is readable by Panoply or IDV for plotting.
pncgen -f bpch /path/to/ctm.bpch /path/for/output/ctm.nc
Viewing CDL is often useful to know what you have (dimensions, variable, etc)
pncdump -f bpch --header /path/to/ctm.bpch
You can also view all the data as ascii. To avoid overload you might look just at the zonal averages at the surface.
pncdump -f bpch -r time,mean -s layer47,0 -r longitude,mean /path/to/ctm.bpch
pncmap.py -f bpch -r time,mean -s layer47,0 /path/to/ctm.bpch /path/for/output/ctm
Example averages time and space to make a vertical profile. The user then select the variable and chooses to make a profile. The mean could be replaced by max, std, var, median, etc.
pncvertprofile.py -f bpch -r time,mean -r latitude,mean -r longitude,mean /path/to/ctm.bpch /path/for/output/ctm
Example extracts groundlayer Ox at two locations and displays it on the screen or makes an output file
pncdump -f bpch --extract=-74,24/-90,30 -s layer47,0 -v Ox /path/to/ctm.bpch
pncgen -f bpch --extract=-74,24/-90,30 -s layer47,0 -v Ox /path/to/ctm.bpch /path/for/output/ctm.nc
This example increases emissions from just one file by 50%
$pncload -f bpch,noscale=True /path/to/emis.bpch
>>> NOx_0[:] *= 1.50
>>> ifile0.sync()
>>> exit()
The simple command that I emailed you does not, but that is easy to do with PNC processing. To provide pressure integration, you'll need PEDGE-$_PSURF in your output. I'll start by assuming you have the latest PseudoNetCDF (after 2016-11-01) and that you have the PEDGE-$_PSURF in the output at the same time-resolution as the species to be averaged.
In this example, I'll work with a 72 layer 4x5 benchmark and I'll create an pressure averaged value for ozone between 491hPa and 288.9hPa. Each line prefixed with a number and dollar sign, is a bash command and can be run separately, but must be run in this order. Below the commands is a detailed description for each.
First, create a pressure averaging script. The script is a reusable way to pressure average and create some derived variables for documentation. Note that I have asked it to print the nominal pressure range (p0 = 1013.25) and the average pressure range for use later. That way, you can change the slice definitions (in step two) and know what you're working with.
echo "dP=np.diff(PSURF[:], axis = 1)
dP.dimensions=('time','layer72','latitude','longitude')
O3w=(O3*dP).sum(1)/dP.sum(1)
O3w.dimensions=('time','latitude','longitude')
TOPP=PSURF[:, -1]
TOPP.dimensions=('time','latitude','longitude')
BOTP=PSURF[:,0]
BOTP.dimensions=('time','latitude','longitude')
print('Nominal Pressure Range: %4g - %4g' % (etai_pressure[0],
etai_pressure[-1]))
print('Average Pressure Range: %4g - %4g' % (BOTP.mean(), TOPP.mean()))
" > pavg.pncexpr
In the second command, I am doing a bunch of stuff. (a) First letting it know this is a NATIVE grid and telling it to open without groups for IJ-AVG-$ and PEDGE-$. That is necessary for the pressure averager to work. (b) I am telling it to slice layer72 from level 23 to 28 and layer73 from 23 to 29 since layer73 (used for pressure edges) starts at the bottom. (c) The next step (--exprscript) is using the expression from step 1 in the context of the input file.
pncgen.py -O -v O3w,dP,TOPP,BOTP --
--format="bpch2,nogroup=('IJ-AVG-$',
'PEDGE-$'),vertgrid='GEOS-5-NATIVE'" --slice=layer72,23,28
--slice=layer_bounds,23,29 --slice=layer73,23,29
--exprscript=pavg.pncexpr trac_avg.geosfp_4x5_benchmark.201307010000
trac_avg.geosfp_4x5_benchmark.201307010000.nc
Next create a plot with a title derived from the bottom and top nominal pressure level coordinates (see etai_pressure and etam_pressure in output). The output is saved as ctm_layeravg12-24_O3w.png. If this is your first time working with matplotlib on your server, you might get an easily fixed* error about X11 or DISPLAY.
pncmap.py --axes-keywords="title='At Nominal Pressure 491-289hPa
(surface=1013.25)'" --norm="Normalize(30,100)" -v O3w -r time,mean
trac_avg.geosfp_4x5_benchmark.201307010000.nc ctm_layeravg12-24
*You can fix a DISPLAY or X11 error by making sure your ~/.config/matplotlib/matplotlibrc file has "backend: Agg" (without the quotes) in it. If need be, create the file and add just that line. If it already exists, edit the backend line(s).