Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / source codes / PyNGL example ICON triangles with different projections
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

PyNGL example ICON triangles with different projections

In this script the projection is changed from Mollweide to Winkel Tripel within a loop.

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

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()

Plot results:

PyNGL unstructured ICON triangles projections w400 PyNGL unstructured ICON triangles 2 projections w400

 

 

 

 

 

Artikelaktionen