source: trunk/UTIL/PYTHON/simple_scripts/2dplot.py @ 1242

Last change on this file since 1242 was 943, checked in by aslmd, 12 years ago

clean and organized UTIL/PYTHON folder

  • Property svn:executable set to *
File size: 1.7 KB
Line 
1#! /usr/bin/env python
2
3from netCDF4 import Dataset
4import matplotlib.pyplot as mpl
5import numpy as np
6from myplot import reducefield,getfield,getcoorddef,calculate_bounds,bounds,fmtvar,ptitle,makeplotres
7from matplotlib.pyplot import contourf,colorbar,show,xlabel,ylabel
8from matplotlib.cm import get_cmap
9
10name = "wrfout_d01_9999-09-09_09:00:00_z"
11itime = 12
12#itime = 1
13ndiv = 10
14zey = 0
15var = "W"
16vmin = -1.
17vmax = 1.
18title = "Vertical velocity"
19var = "Um"
20vmin = -2.
21vmax = 18.
22title = "Horizontal velocity"
23#var = "tk"
24#vmin = 150.
25#vmax = 170.
26#title = "Atmospheric temperature"
27
28name = "wrfout_d01_2024-03-22_01:00:00_z"
29itime = 12
30ndiv = 10
31zey = 120
32var = "VMR_ICE"
33title = "Volume mixing ratio of water ice [ppm]"
34vmin = -0.5
35vmax = 400.
36
37nc = Dataset(name)
38
39what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d2=zey )
40
41
42y = nc.variables["vert"][:]
43
44horinp = len(what_I_plot[0,:])
45x = np.linspace(0.,horinp*500.,horinp) / 1000.
46xlabel("Horizontal coordinate (km)")
47
48horinp = len(what_I_plot[0,:])
49x = np.linspace(0.,horinp,horinp) 
50xlabel("x grid points")
51
52
53zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
54#if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar))
55#else:                           palette = get_cmap(name=colorb)
56palette = get_cmap(name="jet")
57what_I_plot = bounds(what_I_plot,zevmin,zevmax)
58zelevels = np.linspace(zevmin,zevmax)
59contourf( x, y, what_I_plot, zelevels, cmap = palette )
60colorbar(fraction=0.05,pad=0.03,format=fmtvar(var),\
61                       ticks=np.linspace(zevmin,zevmax,ndiv+1),\
62                       extend='neither',spacing='proportional')
63ptitle(title)
64ylabel("Altitude (m)")
65makeplotres(var+str(itime),res=200.,disp=False)
66#show()
Note: See TracBrowser for help on using the repository browser.