import komod
import Nio
from IPython.display import Image
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
Contain mostly set of wrapper functions for map plotting with PyNGL. Can be used with any 2D data, not necessarily MITgcm.
coord2d - Convert 1d coordinates to 2d coordinates
colormap - Define custom colormaps
arctpl - plot contours of 2D, 3D or 4D field in the Arctic region from ndarray
arctpltnc - plot contours of 2D, 3D or 4D field in the Arctic region from netCDF file
globpltnc - plot contours of 2D, 3D or 4D field on the global map from netCDF file. You can also specify custom region.
pltgrd - plot model grid from lat and lon arrays
pltgrdnc - plot model grid from netCDF file
plt_vectors - plot vector data from 2D fields
plt_vectors_scalars - plot vectors and scalar filed (e.g. wind direction and speed)
plt_vectors_colors - plot colored vectors
Plot variable from array in the Arctic region. Let's use netCDF file generated by komod.mitopen.
ff =Nio.open_file('MIT_output_2d.nc')
area = ff.variables['AREAtave'][:]
lat = ff.variables['latitude'][:]
lon = ff.variables['longitude'][:]
komod.arctpl(lon, lat, area[0,:,:])
If you have gv installed, then window with the map should pop up. If not, let's convert ps file to png:
!convert -scale 50% output.ps output.png
Image(filename='output.png')
Let's try to specify more parameters. We can choose another region to plot (by default it's Arctic). List of available regions can be found in the reg function. You can easily add new region inside this function.
komod.arctpl(lon, lat, area[0,:,:], region='FramStAnna')
!convert -scale 50% output.ps output.png
Image(filename='output.png')
If you want to quickly plot some region without all this trouble with defining it in reg function, you can do this simply bu using region="Global" and defining min/max lon/lat:
komod.arctpl(lon, lat, area[0,:,:], region='Global', minLon=10, maxLon=50 , minLat=40 , maxLat=75 )
!convert -scale 50% output.ps output.png
Image(filename='output.png')
However polar regions will look ugly and there is still problem with plotting over 0 meridian. For longitudes only 0-360 format is supported, so something like minLon=-20, maxLon=20 will not work :(
If you don't specify any particular field, like we do with area[0,:,:], arctpl should plot all fields you provide, assuming that your dimensions are [time, level, :, :] or [time/level, :, :]. We try it with NCEP reanalysis.
!wget ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.dailyavgs/surface/air.sig995.2012.nc
!cdo seasmean air.sig995.2012.nc air.sig995.2012_sm.nc
--2013-03-11 20:53:43-- ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.dailyavgs/surface/air.sig995.2012.nc => `air.sig995.2012.nc' Resolving ftp.cdc.noaa.gov (ftp.cdc.noaa.gov)... 140.172.38.117 Connecting to ftp.cdc.noaa.gov (ftp.cdc.noaa.gov)|140.172.38.117|:21... connected. Logging in as anonymous ... Logged in! ==> SYST ... done. ==> PWD ... done. ==> TYPE I ... done. ==> CWD (1) /Datasets/ncep.reanalysis.dailyavgs/surface ... done. ==> SIZE air.sig995.2012.nc ... 7700344 ==> PASV ... done. ==> RETR air.sig995.2012.nc ... done. Length: 7700344 (7.3M) (unauthoritative) 100%[======================================>] 7,700,344 1.00M/s in 9.5s 2013-03-11 20:53:56 (794 KB/s) - `air.sig995.2012.nc' saved [7700344] cdo seasmean: Processed 3847392 values from 1 variable over 366 timesteps. ( 0.10s )
ncep = Nio.open_file('air.sig995.2012_sm.nc', 'r')
lat_ncep = ncep.variables['lat'][:]
lon_ncep = ncep.variables['lon'][:]
ncep_air = ncep.variables['air'][:]
lon_ncep2, lat_ncep2 = np.meshgrid(lon_ncep, lat_ncep)
One has to apply scale_factor and add_offset to the NCEP data, and we convert from K to C. Yes and we add cyclic point, in order to get rid of ugly white line at 0 meridian.
komod.arctpl(lon_ncep2, lat_ncep2, ((ncep_air[:,:,:]*0.01)+512.81)-273.15, add_cyclic=True)
Result should be *.ps file with 5 maps, each for one season. Should look something like this:
Image(filename='5plots.png')
Color scale is intentionally made the same among plots. Sorry to atmosphere guys, no values over land.
Same as arctplt, but takes netCDF file and variable name as input values, and tries to get as much information as it can from attributes. If you lucky you get the plot, otherwise, you have to cook the file a bit before feed it to this function.By default plots every field it can squeeze from the file.
komod.arctpltnc('MIT_output_2d.nc', 'AREAtave', lon='longitude', lat='latitude', showfig=False)
By this point I already get annoyed by popping up windows with ps files, so I switch this option off (showfig = False).
!convert -scale 50% output.ps output.png
Image(filename='output.png')
Here we get the time in the header.
Let's try 3D variable. Warning: THIS CAN BE VERY SLOW (about 3-5 minutes)!
komod.arctpltnc('MIT_output_3d.nc', 'Ttave', lon='longitude', lat='latitude', llevel='z', levon=True)
You will get all levels plotted in the same color scale. Option levon=True add depth of the level to the header.
Image(filename='many_levels.png')
Function that displays model grid. In general requires only lat and lon (2d). By default region is "Global".
komod.pltgrd(lon,lat, region = 'FramStAnna')
!convert -scale 50% -rotate -90 grid.ps grid.png
Image(filename='grid.png')
Grid is very dense even for such a relatively small region. We can make it less dense by specifying every parameter.
komod.pltgrd(lon,lat, region = 'FramStAnna', every=5)
!convert -scale 50% -rotate -90 grid.ps grid.png
Image(filename='grid.png')
There are still problems around the North Pole :(
Same as the function above, but takes netCDF file as input. It assumes that the names of lon and lat variables are 'lon' and 'lat' (you can change that):
komod.pltgrdnc('air.sig995.2012.nc', 'air')
!convert -scale 50% -rotate -90 grid.ps grid.png
Image(filename='grid.png')
Plot model grid from lat and lon arrays and line defined by coordinates of the first and the last point.
komod.pltgrd_line(lon, lat, lon1=30, lat1=80, lon2=30, lat2=85, every=3, region = 'FramStAnna')
(array([ 30. , 30.00000191, 30.00000191, 29.99999809, 29.99999809, 30.00000191, 29.99999809, 30.00000191, 29.99999809, 30. ], dtype=float32), array([ 80. , 80.55554962, 81.11110687, 81.66666412, 82.22222137, 82.77777863, 83.33332825, 83.88889313, 84.44444275, 85. ], dtype=float32))
As you can see it also returns coordinates of points along this line (by default 10). Points are calculated by Ngl.gc_interp function.
!convert -scale 50% grid_line.ps grid_line.png
Image(filename='grid_line.png')
Same as above, but instead of end points coordinates took coordinates of points defining the coordinates of the polyline.
Plot transect. Input data for this function calculated by the get_transect function.
First we have to open some 3d data:
import Nio
f3d = Nio.open_file('MIT_output_3d.nc', 'r')
lat3 = f3d.variables['latitude'][:]
lon3 = f3d.variables['longitude'][:]
temp = f3d.variables['Ttave'][:]
lev = f3d.variables['z'][:] # depths
Then we extract transect from the data by get_transect function. Input parameters are 2d arrays of lat and lon, and end points of the transect. First this function will calculate coordinates of the equally distributed points (default is 10) along the line specified by the coordinates of the beginning and the end of the transect. Then it finds grid points that are closest to this points of the transect and get vertical profiles from each of them (collected in data_prof). There is no interpolation, so its a very straight forward approach, that can lead to very strange results on coarse grids. Also the distances between first point of the transect and the rest of the points are calculated (collected in x_kilometers).
data_prof, x_kilometers, m_grid, n_grid = komod.get_transect(lat3, lon3, temp, lon1=30, lat1=80, lon2=30, lat2=85)
(252, 299, 79.951835632324219, 29.880311965942383) (254, 294, 80.590774536132812, 30.156497955322266) (255, 290, 81.073844909667969, 29.921560287475586) (257, 285, 81.714401245117188, 30.258316040039062) (258, 281, 82.198867797851562, 30.008869171142578) (259, 276, 82.792793273925781, 29.401870727539062) (261, 272, 83.326507568359375, 30.160533905029297) (262, 267, 83.921958923339844, 29.471311569213867) (264, 263, 84.456298828125, 30.409765243530273) (265, 259, 84.943344116210938, 30.076766967773438) [ 0. 71.2383141 125.11606104 196.56490763 250.576502 317.20506821 377.42129338 444.17831094 504.52632612 558.79284924] [ 252. 254. 255. 257. 258. 259. 261. 262. 264. 265.] [ 299. 294. 290. 285. 281. 276. 272. 267. 263. 259.]
Now we can plot the transect:
komod.plt_transect(distances=x_kilometers , levels=lev, vvalues=data_prof)
[ 0. 71.2383141 125.11606104 196.56490763 250.576502 317.20506821 377.42129338 444.17831094 504.52632612 558.79284924]
!convert -scale 50% -rotate -90 output.ps output.png
Image(filename='output.png')
To look at the location of the transect, you can use pltgrd_line function.