EarthPy

Plot grid and transect with PyNGL and komod

Task:

Draw a grid and a transect with PyNGL module

Solution:

komod module

Notebook file

This is a short follow up of the previous post about komod module, that is essentially a set of wrapper functions for PyNGL module. Here I am going to show how to plot a grid of your model and how to draw a transect. I am going to use the same data set af before: mean temperature from the World Ocean Atlas 2009 (5 deg. resolution).

Import modules:

In [2]:

import komod
import Nio
from IPython.display import Image
import numpy as np

Get variables from the netCDF file:

In [3]:

ff =Nio.open_file('temperature_annual_5deg.nc')
temp = ff.variables['t_mn'][:]
lat  = ff.variables['lat'][:]
lon  = ff.variables['lon'][:]
lev  = ff.variables['depth'][:]

Create lat/lon mesh grid:

In [4]:

lon_m, lat_m = np.meshgrid(lon,lat)

Now we can plot our grid. This function use centers of the grid boxes, not corners, so result might not be a perfectly accurate representation of your grid, but rather rough estimate of its shape and density. Only necessary input parameters are 2d arrays of lats and lons:

In [5]:

komod.pltgrd(lon_m,lat_m, region = 'Global', add_cyclic=True)

In [6]:

!convert -scale 70% -rotate -90 grid.ps grid.png
Image(filename='grid.png')

Out[6]:

In this case everything is very simple and rectangular. Let's have a look at more interesting case. We will use a grid from one of the MPIOM configurations, and you can download this sample file from here.

In [7]:

g = Nio.open_file('GR15s.nc')
lat_g = g.variables['grid_center_lat'][:]
lon_g = g.variables['grid_center_lon'][:]

In [8]:

komod.pltgrd(lon_g,lat_g, region = 'Global', add_cyclic=True)
!convert -scale 70% -rotate -90 grid.ps grid.png
Image(filename='grid.png')

Out[8]:

Very dense grid, you can hardly see any details. We can show only every n$^{th}$ element of grid by specifying every parameter:

In [9]:

komod.pltgrd(lon_g,lat_g, region = 'Global', every=3, add_cyclic=True)
!convert -scale 70% -rotate -90 grid.ps grid.png
Image(filename='grid.png')

Out[9]:

Yes I know there is an ugly white in the middle, but fixing it would require a bit more effort than I want to invest in this function :) At the end it do the thing it was created for - show how the grid looks like.

You can look at the regions as well:

In [10]:

komod.pltgrd(lon_g,lat_g, region = 'Arctic', every=1, add_cyclic=True)
!convert -scale 70% -rotate -90 grid.ps grid.png
Image(filename='grid.png')

Out[10]:

You can see that the Fram Strait have much better resolution than the Bering Strait. North pole is still a problem :(

Now we are almost ready to plot a transect. But first let's have a look at its position on the map and how its related to the grid. For this we use komod.pltgrd_line function. We have to provide lat/lon arrays with coordinates and also start and end points of your transect:

In [11]:

komod.pltgrd_line(lon_m, lat_m, lon1=335, lat1=50, lon2=335, lat2=-50, every=1, region = 'Global', add_cyclic=True )

Out[11]:

(array([ 335.,  335.,  335.,  335.,  335.,  335.,  335.,  335.,  335.,  335.], dtype=float32),
 array([ 50.        ,  38.88888931,  27.77777672,  16.66666412,
         5.55555248,  -5.55555916, -16.66667175, -27.77777672,
       -38.88889313, -50.        ], dtype=float32))

In [12]:

!convert -scale 70% -rotate -90 grid_line.ps grid_line.png
Image(filename='grid_line.png')

Out[12]:

It's maybe a bit hard to see, but there is a line now at the center of the Atlantic Ocean. This function uses Ngl.gc_interp to calculate positions of the points between the end points. Number of intermediate points is 10 by default, and can be changed (npoints parameter).

After we find a good position for our transect, we have to get the data to plot it. This is done with komod.get_transect and you have to provide lat/lon arrays, 3D array with your data and start/end points of the transect. Two main variables that this function return is data_prof with vertical profiles at the points along your line (10 point by default) and x_kilometers - the distance (in kilometers) between your start point and every point of your transect. This variables will be used as an input for plotting function.

In [14]:

data_prof, x_kilometers, m_grid, n_grid = komod.get_transect(lat_m, lon_m, temp[0,0:20,:,:], lon1=335, lat1=50, lon2=335, lat2=-50)

Finally we can plot our transect:

In [15]:

komod.plt_transect(distances=x_kilometers , levels=lev[0:20], vvalues=data_prof, datamin = -2., datamax = 30., datastep = 1)

In [16]:

!convert -scale 50% -rotate -90 output.ps output.png
Image(filename='output.png')

Out[16]:

With only 10 points between start and end the transect look kind of rough. Let's increase number of points:

In [17]:

data_prof, x_kilometers, m_grid, n_grid = komod.get_transect(lat_m, lon_m, temp[0,0:20,:,:], 
                                                             lon1=335, lat1=50, lon2=335, lat2=-50, npoints=30)
komod.plt_transect(distances=x_kilometers , levels=lev[0:20], vvalues=data_prof, datamin = -2., datamax = 30., datastep = 1)
!convert -scale 50% -rotate -90 output.ps output.png
Image(filename='output.png')

Out[17]:

Now it look squared :) The reason is that for some points on the line the same grid boxes are chosen, so we have repetitions. You can see it in the x_kilometers variable - some distances occur twice:

In [18]:

x_kilometers

Out[18]:

array([     0.        ,    555.99391722,   1111.98783444,   1667.98175166,
         1667.98175166,   2223.97566888,   2779.9695861 ,   2779.9695861 ,
         3335.96350332,   3891.95742055,   3891.95742055,   4447.95133777,
         5003.94525499,   5003.94525499,   5559.93917221,   6115.93308943,
         6671.92700665,   6671.92700665,   7227.92092387,   7783.91484109,
         7783.91484109,   8339.90875831,   8895.90267553,   8895.90267553,
         9451.89659275,  10007.89050997,  10007.89050997,  10563.88442719,
        11119.87834441,  11675.87226164])

We can filter out this repetition by norep parameter:

In [19]:

data_prof, x_kilometers, m_grid, n_grid = komod.get_transect(lat_m, lon_m, temp[0,0:20,:,:], norep=True,
                                                             lon1=335, lat1=50, lon2=335, lat2=-50, npoints=30)
komod.plt_transect(distances=x_kilometers , levels=lev[0:20], vvalues=data_prof, datamin = -2., datamax = 30., datastep = 1)
!convert -scale 50% -rotate -90 output.ps output.png
Image(filename='output.png')

Out[19]:

Better :)

In [20]:

x_kilometers

Out[20]:

array([     0.        ,    555.99391722,   1111.98783444,   1667.98175166,
         2223.97566888,   2779.9695861 ,   3335.96350332,   3891.95742055,
         4447.95133777,   5003.94525499,   5559.93917221,   6115.93308943,
         6671.92700665,   7227.92092387,   7783.91484109,   8339.90875831,
         8895.90267553,   9451.89659275,  10007.89050997,  10563.88442719,
        11119.87834441,  11675.87226164])

One last thing to mention is that you can change aspect ratio of your plot by changing values of vpWidthF and vpHeightF parameters:

In [21]:

komod.plt_transect(distances=x_kilometers , levels=lev[0:20], 
                   vvalues=data_prof, datamin = -2., datamax = 30., datastep = 2,
                   vpWidthF=0.5, vpHeightF=0.3)
!convert -scale 50% -rotate -90 output.ps output.png
Image(filename='output.png')

Out[21]:

There are comments.

EarthPy

Other articles

Plot maps with PyNGL and komod

Other articles

links

social