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:
