PyNGL example ICON triangles with different projections#

Software requirements:

  • Python 2.7.x

  • Numpy 1.9.2

  • PyNGL/PyNIO 1.5.0

Run the ICON triangles projections example script:

python PyNGL_unstructured_ICON_triangles_projections.py

Script PyNGL_unstructured_ICON_triangles_projections.py:

'''
 DKRZ PyNGL script:
                PyNGL_unstructured_ICON_triangles_projections.py

 Grid type:     unstructured grid
 Data/Model:    ICON, 20480 cells

 Settings:      polygons, labelbar, wksColorMap

 Projections:   Robinson, Mollweide, WinkelTripel, CylindricalEquidistant

 Info:          The Ngl.polygon function is used
                because it's faster!

 Plots:
        Mollweide:
                plot_PyNGL_unstructured_ICON_triangles_projections.000001.png
        Winkel Tripel:
                plot_PyNGL_unstructured_ICON_triangles_projections.000002.png
        Cylindrical Equidistant:
                plot_PyNGL_unstructured_ICON_triangles_projections.000003.png

 25.02.2016   meier-fleischer(at)dkrz.de
'''
import numpy as np
import math, time, sys, os
import Nio, Ngl

#-----------------------------------------------------------------------
#-- Function: main
#-----------------------------------------------------------------------
def main():

    t1 = time.time()                                    #-- retrieve start time

    #-- projections
    proj = ["Mollweide","WinkelTripel"]

    #--  define variables
    diri, fname  = '$HOME/data/ICON/', 'ta_ps_850.nc'   #-- data path and file name
    gname        = 'grids/r2b4_amip.nc'                 #-- grid info file
    VarName      = 'ta'                                 #-- variable name

    #--  open file and read variables
    f = Nio.open_file(diri + fname,'r')                 #-- add data file
    g = Nio.open_file(diri + gname,'r')                 #-- add grid file (not contained in data file!!!)

    #-- read a timestep of 'ta'
    variable =  f.variables['ta']                       #-- first time step, lev, ncells
    data     =  variable[0,0,:]                         #-- ta [time,lev,ncells]; miss _FillValue
    var      =  data - 273.15                           #-- convert to degrees Celsius; miss _FillValue

    #-- define _FillValue and missing_value if not existing
    missing = -1e20

    if not hasattr(var,'_FillValue'):
       var._FillValue  =  missing                       #-- set _FillValue
    if not hasattr(var,'missing_value'):
       var.missing_value =  missing                     #-- set missing_value

    varM = np.ma.array(var, mask=np.equal(var,missing)) #-- mask array with missing values
    nummissing = np.count_nonzero(varM.mask)            #-- number of missing values

    #-- set data intervals, levels, labels, color indices
    varMin, varMax, varInt = -32, 28, 4                 #-- set data minimum, maximum, interval

    levels   =  range(varMin,varMax,varInt)             #-- set levels array
    nlevs    =  len(levels)                             #-- number of levels
    labels   = ['{:.2f}'.format(x) for x in levels]     #-- convert list of floats to list of strings
    colors   =  range(2,nlevs+6)                        #-- create color indices

    #-- print info to stdout
    print ''
    print 'min/max:          %.2f' %np.min(varM) + ' /' + ' %.2f' %np.max(varM)
    print ''
    print 'varMin:           %3d' %varMin
    print 'varMax:           %3d' %varMax
    print 'varInt:           %3d' %varInt
    print ''
    print 'missing_value:    ', missing
    print 'missing values:   ', nummissing
    #-------------------------------------------------------------------
    #-- define the x-, y-values and the polygon points
    #-------------------------------------------------------------------
    rad2deg = 45./np.arctan(1.)                         #-- radians to degrees

    x, y       =  g.variables['clon'][:], g.variables['clat'][:]
    vlon, vlat =  g.variables['clon_vertices'][:], g.variables['clat_vertices'][:]

    x, y       =  x*rad2deg,  y*rad2deg                 #-- cell center, lon, lat
    vlat, vlon =  vlat*rad2deg, vlon * rad2deg          #-- cell latitude/longitude vertices
    ncells, nv =  vlon.shape                            #-- ncells: number of cells; nv: number of edges

    #-- print information to stdout
    print ''
    print 'cell points:      ', nv
    print 'cells:            ', str(ncells)
    print ''

    #-- rearrange the longitude values to -180.-180.
    def rearrange(vlon):
        less_than    = vlon < -180.
        greater_than = vlon >  180.
        vlon[less_than]    = vlon[less_than] + 360.
        vlon[greater_than] = vlon[greater_than] - 360.
        return vlon

    vlon = rearrange(vlon)                              #-- set longitude values to -180.-180. degrees

    print 'min/max vlon:     ', np.min(vlon), np.max(vlon)
    print 'min/max vlat:     ', np.min(vlat), np.max(vlat)
    print ''

    #-- open a workstation for second plot:  triangles plot
    wkres       =  Ngl.Resources()
    wkres.wkColorMap =  'WhiteBlueGreenYellowRed'       #-- choose colormap
    wkres.wkWidth, wkres.wkHeight  =  1024, 1024
    wks1_type   = 'png'
    wks1 =  Ngl.open_wks(wks1_type,'plot_PyNGL_unstructured_ICON_triangles_projections',wkres)

    #-- define colormap
    cmap     =  Ngl.retrieve_colormap(wks1)             #-- RGB ! [256,3]
    ncmap    =  cmap.shape[0]                           #-- number of colors
    colormap =  cmap[:ncmap:12,:]                       #-- select every 13th color
    ncol     =  colormap.shape[0]
    colormap[20] = ([1.,1.,1.])                         #-- white for missing values

    print 'colors index:     ', colors
    print ''
    print 'levels:           ', levels
    print 'labels:           ', labels
    print ''
    print 'nlevs:            %3d' %nlevs
    print 'ncols:            %3d' %ncol
    print ''

    #-- overwrite resources of wks1
    setlist                    =  Ngl.Resources()
    setlist.wkColorMap         =  colormap              #-- set color map to new colormap array
    setlist.wkBackgroundColor  = 'white'                #-- has to be set when wkColorMap is set to colormap array
    setlist.wkForegroundColor  = 'black'                #-- has to be set when wkColorMap is set to colormap array
    Ngl.set_values(wks1,setlist)

    #-- set map resources
    mpres                             =  Ngl.Resources()
    mpres.nglDraw                     =  False          #-- turn off plot draw and frame advance. We will
    mpres.nglFrame                    =  False          #-- do it later after adding subtitles.
    mpres.mpGridAndLimbOn             =  False
    mpres.mpGeophysicalLineThicknessF =  2.

    #-- loop over projection list

    for ip in xrange(len(proj)):
       print 'Projection: '+proj[ip]

       if proj[ip] == 'Mollweide' or proj[ip] == 'WinkelTripel':
          mpres.mpGridAndLimbOn     =  True             #-- turn on lat/lon/limb lines
          mpres.mpGridLineColor     = 'transparent'     #-- we don't want lat/lon lines
          mpres.pmTickMarkDisplayMode = 'Never'         #-- don't draw tickmark border (box) around plot
          mpres.mpPerimOn           =  False            #-- don't draw the box around the plot

       mpres.mpProjection       =  proj[ip]             #-- set projection
       mpres.tiMainString       = 'ICON - PyNGL '+proj[ip]+' projection' #-- title string

       #-- create and draw the basic map
       map = Ngl.map(wks1,mpres)
       Ngl.draw(map)

       #-- assign and initialize array which will hold the color indices of the cells
       gscolors = -1*(np.ones((ncells,),dtype=np.int))     #-- assign array containing zeros; init to transparent: -1

       #-- set color index of all cells in between levels
       for m in xrange(0,nlevs):
           vind = []                                       #-- empty list for color indices
           for i in xrange(0,ncells-1):
               if (varM[i] >= levels[m] and varM[i] < levels[m+1]):
                  gscolors[i] = colors[m+1]
                  vind.append(i)
           print 'finished level %3d' % m , ' -- %5d ' % len(vind) , ' polygons considered - gscolors %3d' % colors[m]
           del vind

       gscolors[varM < varMin]         =  colors[0]        #-- set color index for cells less than level[0]
       gscolors[varM >= varMax]        =  colors[(nlevs-1)+2] #-- set color index for cells greater than levels[nlevs-1]
       gscolors[np.nonzero(varM.mask)] =  20               #-- set color index for missing values

       #-- set polygon resources
       pgres                   =  Ngl.Resources()
       pgres.gsEdgesOn         =  True                     #-- draw the edges
       pgres.gsFillIndex       =  0                        #-- solid fill
       pgres.gsLineColor       = 'black'                   #-- edge line color
       pgres.gsLineThicknessF  =  0.7                      #-- line thickness
       pgres.gsColors          =  gscolors                 #-- use color array
       pgres.gsSegments        =  range(0,len(vlon[:,0])*3,3) #-- define segments array

       x1d, y1d = np.ravel(vlon), np.ravel(vlat)           #-- convert to 1D-arrays

       #-- add polygons to map
       polyg  = Ngl.add_polygon(wks1,map,x1d,y1d,pgres)

       #-- add a labelbar
       lbres                   =  Ngl.Resources()
       lbres.vpWidthF          =  0.85
       lbres.vpHeightF         =  0.15
       lbres.lbOrientation     = 'Horizontal'
       lbres.lbFillPattern     = 'SolidFill'
       lbres.lbMonoFillPattern =  21                       #-- must be 21 for color solid fill
       lbres.lbMonoFillColor   =  False                    #-- use multiple colors
       lbres.lbFillColors      =  colors                   #-- indices from loaded colormap
       lbres.lbBoxCount        =  len(colormap[colors,:])
       lbres.lbLabelFontHeightF=  0.014
       lbres.lbLabelAlignment  = 'InteriorEdges'
       lbres.lbLabelStrings    =  labels

       if proj[ip] == 'WinkelTripel':
          print '--> move labelbar downward (WT)'
          lbx, lby  = 0.1, 0.14
       else:
          lbx, lby  = 0.1, 0.24

       lb = Ngl.labelbar_ndc(wks1,nlevs+1,labels,lbx,lby,lbres)

       #-- maximize and draw the plot and advance the frame
       Ngl.maximize_plot(wks1, map)
       Ngl.draw(map)
       Ngl.frame(wks1)

    #-- get wallclock time
    t2 = time.time()
    print ''
    print 'Wallclock time:  %0.3f seconds' % (t2-t1)
    print ''

    #-- done
    Ngl.end()

if __name__ == '__main__':
    main()

Results:

PyNGL unstructured ICON triangles projections w400

PyNGL unstructured ICON triangles 2 projections w400