DKRZ NCL ICON add_polygon triangles example#
DKRZ NCL script:
;-------------------------------------------------------------------
;-- Script: ICON_add_polygon_triangles.ncl
;--
;-- Description: - read ICON data and plot the triangles using
;-- contour plot levels and colors informations
;-- - global plot
;-- - missing values plotted with "gray"
;-- - check if all data values is missing
;--
;-- Settings: - cnFillMode = "CellFill"
;-- - sfXCellBounds = vlon
;-- - sfYCellBounds = vlat
;-- - sfXArray = x
;-- - sfYArray = y
;--
;-- Notice: Run time of comparable scripts:
;-- NCL script: 1.5s
;-- PyNGL/PyNIO script: 4.3s
;-- Matplotlib script: 28.1s
;--
;-- NCL Version: 6.2.0
;--
;-- 29.02.16 meier-fleischer(at)dkrz.de
;-------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;----------------------------------------------------------------------;
;-- Procedure: add_lab --;
;-- Create a new and add labelbar --;
;----------------------------------------------------------------------;
undef("add_lab")
procedure add_lab(wks,map,cmap,levels)
begin
getvalues map
"vpWidthF" : vpw ;-- retrieve the viewport width
"vpHeightF" : vph ;-- retrieve the viewport height
"vpXF" : vpx ;-- retrieve the viewport x pos
"vpYF" : vpy ;-- retrieve the viewport y pos
end getvalues
;-- set the labelbar labels, width, height and position
labs = levels+"" ;-- labelbar labels
nlevs = dimsizes(labs) ;-- number of labels
lbwidth = 0.7 ;-- labelbar width
lbheight = 0.08 ;-- labelbar height
lbx = ((1.0-(lbwidth+vpx))/2) + vpx/2 ;-- labelbar x position
lby = ((1.0-(vpy-vph))/2) - 0.13 ;-- labelbar y position
;-- set global labelbar resources
lbres = True
lbres@gsnFrame = False ;-- don't advance frame
lbres@vpWidthF = lbwidth ;-- width of labelbar
lbres@vpHeightF = lbheight ;-- height of labelbar
lbres@vpXF = lbx ;-- labelbar x-position
lbres@vpYF = lby ;-- labelbar y-position
lbres@lbPerimOn = False ;-- no label bar box
lbres@lbOrientation = "Horizontal" ;-- orientation
lbres@lbLabelFontHeightF = 0.010 ;-- label font height
lbres@lbLabelAlignment = "InteriorEdges" ;-- where to label
lbres@lbMonoFillPattern = True ;-- fill solid
lbres@lbFillColors = cmap ;-- use colors
gsn_labelbar_ndc(wks,nlevs+1,labs,lbx,lby,lbres) ;-- draw transparent labelbar
end
;=======================================================================
;-- MAIN PROGRAM
;=======================================================================
begin
startt = get_cpu_time()
DataFileName = "$HOME/data/ICON/ta_ps_850.nc" ;-- data file
gridinfofile = "$HOME/data/ICON/grids/r2b4_amip.nc" ;-- grid info file
VarName = "ta" ;-- variable name
File = addfile( DataFileName, "r" ) ;-- add data file
GridInfoFile = addfile( gridinfofile, "r" ) ;-- add grid file (not contained in data file!!!)
var = File->$VarName$(time|0,lev|0,ncells|:) ;-- set variable with dims: (time,ncell)
var = var - 273.15 ;-- convert to degrees Celsius
title = "ICON: Surface temperature ["+var@units+"]" ;-- title string
varMin = -32. ;-- data minimum
varMax = 24. ;-- data maximum
varInt = 4. ;-- data increment
lon_min = -180.0 ;-- minimum longitude
lon_max = 180.0 ;-- maximum longitude
lat_min = -90.0 ;-- minimum latitude
lat_max = 90.0 ;-- maximum latitude
;-------------------------------------------------------------------
;-- define the x-, y-values and the polygon points
;-------------------------------------------------------------------
rad2deg = 45./atan(1.) ;-- radians to degrees
x = GridInfoFile->clon *rad2deg ;-- cell center, lon
y = GridInfoFile->clat *rad2deg ;-- cell center, lat
x!0 = "lon" ;-- set named dimension lon
y!0 = "lat" ;-- set named dimension lat
x@units = "degrees_east" ;-- set lon units
y@units = "degrees_north" ;-- set lat units
vlon = GridInfoFile->clon_vertices * rad2deg ;-- cell longitude vertices
vlon = where(vlon.lt.0, vlon + 360, vlon) ;-- longitude: 0-360
vlat = GridInfoFile->clat_vertices * rad2deg ;-- cell lattitude vertices
nv = dimsizes(vlon(0,:)) ;-- number of points in polygon
;-- print some information to stdout
print("")
print("Data Min: " + min(var) + " Max: " + max(var))
print("Data longitude Min: " + min(vlon) + " Max: " + max(vlon))
print("Cell points: " + nv)
print("")
;-------------------------------------------------------------------
;-- open workstation
;-------------------------------------------------------------------
wtype = "png" ;-- plot output type
wtype@wkWidth = 1024 ;-- set workstation width in pixel
wtype@wkHeight = 1024 ;-- set workstation height in pixel
wks = gsn_open_wks(wtype,"plot_ICON_add_polygon_triangles") ;-- open a workstation
;-------------------------------------------------------------------
;-- set resources
;-------------------------------------------------------------------
res = True
res@gsnDraw = False ;-- don't draw the plot
res@gsnFrame = False ;-- don't advance the frame
res@gsnMaximize = True ;-- maximize plot output
res@cnFillPalette = "WhiteBlueGreenYellowRed" ;-- choose colormap - no white
res@cnLinesOn = False ;-- don't draw contour lines
res@cnInfoLabelOn = False ;-- switch off contour info label
res@cnFillOn = True ;-- contour fill on
res@cnMonoFillColor = False
res@cnFillColors = (/20, 35, 50, 60, 70, 80, 90, 115, 125, 132, 145, 155, 165, 175, 190, 205/)
res@cnLevelSelectionMode = "ManualLevels" ;-- use same values for both plots
res@cnMinLevelValF = varMin ;-- minimum contour level
res@cnMaxLevelValF = varMax ;-- maximum contour level
res@cnLevelSpacingF = varInt ;-- contour level spacing
res@cnFillMode = "CellFill" ;-- set fill mode
res@sfXArray = x ;-- transform x to mesh scalar field
res@sfYArray = y ;-- transform y to mesh scalar field
res@sfXCellBounds = vlon ;-- needed if set cnFillMode = "CellFill"
res@sfYCellBounds = vlat ;-- needed if set cnFillMode = "CellFill"
res@lbLabelBarOn = False ;-- don't draw a labelbar yet
res@lbLabelStride = 2 ;-- label every other
res@mpFillOn = False ;-- fill map grey
res@mpDataBaseVersion = "MediumRes" ;-- map resolution
res@mpMinLonF = lon_min ;-- sub-region minimum longitude
res@mpMaxLonF = lon_max ;-- sub-region maximum longitude
res@mpMinLatF = lat_min ;-- sub-region minimum latitude
res@mpMaxLatF = lat_max ;-- sub-region maximum latitude
res@mpGreatCircleLinesOn = False ;-- important: v6.2.0 False !!
;-------------------------------------------------------------------
;-- create the contour plot, but don't draw it (we need to get
;-- the colors of the data values)
;-------------------------------------------------------------------
plot = gsn_csm_contour_map(wks,var,res)
;-------------------------------------------------------------------
;-- to plot the filled polygons we need the data levels and their colors
;-------------------------------------------------------------------
getvalues plot@contour
"cnLevels" : levels ;-- retrieve the levels from plot@contour
"cnFillColors" : colors ;-- retrieve the colors from plot@contour
end getvalues
nlevels = dimsizes(levels) ;-- number of levels
;-------------------------------------------------------------------
;-- clear the plot, but keep all the information
;-------------------------------------------------------------------
plot = setColorContourClear(plot,min(var),max(var))
;-------------------------------------------------------------------
;-- add a labelbar
;-------------------------------------------------------------------
add_lab(wks,plot,colors,levels) ;-- add labelbar
;-------------------------------------------------------------------
;-- create color array for triangles
;-------------------------------------------------------------------
ntri = dimsizes(y) ;-- number of triangles
gscolors = new(ntri,string) ;-- create color array (type of colors integer!)
gscolors = "gray" ;-- initialize to black (for missing colors)
;-- set resources for the triangles (polygons)
pres = True
pres@gsEdgesOn = True ;-- turn on edges
pres@gsEdgeThicknessF = 1.
pres@gsFillIndex = 0 ;-- solid fill
;-- set color for data less than given minimum value var_min
vlow = ind(var .lt. levels(0)) ;-- get the indices of values less levels(0)
if (.not. all(ismissing(vlow))) then
gscolors(vlow) = colors(0) ;-- choose color
ntri_calc = dimsizes(vlow) ;-- number of triangles
end if
;-- set colors for all cells in between var_min and var_max
do i = 1, dimsizes(levels) - 1
vind := ind(var .ge. levels(i-1) .and. var .lt. levels(i)) ;-- get the indices of 'middle' values
if (.not. all(ismissing(vind))) then
gscolors(vind) = colors(i) ;-- choose the colors
ntri_calc = ntri_calc + dimsizes(vind) ;-- number of triangles
end if
end do
;-- set color for data greater than given maximum var_max
nc=dimsizes(colors)-1 ;-- get the number of colors minus one
nl=dimsizes(levels)-1 ;-- get the number of levels minues one
vhgh := ind(var .gt. levels(nl)) ;-- get indices of values greater levels(nl)
if (.not. all(ismissing(vhgh))) then
gscolors(vhgh) = colors(nc) ;-- choose color
ntri_calc = ntri_calc + dimsizes(vhgh) ;-- number of triangles
end if
print("--> triangles calculated: "+ ntri_calc)
;-- attach all the triangles using the list of colors
pres@gsColors = gscolors ;-- set colors for polygons
pres@gsSegments = ispan(0,dimsizes(var) * 3,3) ;-- assign segments array
polygon = gsn_add_polygon(wks,plot,ndtooned(vlon),ndtooned(vlat),pres) ;-- draw all triangles
;-------------------------------------------------------------------
;-- write title on top of the plot
;-------------------------------------------------------------------
tires = True ;-- text resources
tires@txFontHeightF = 0.015 ;-- text font size
tires@txJust = "CenterCenter" ;-- text justification
tires@txFont = "courier" ;-- text font
xl = 0.5 ;-- text x position
yb = 0.85 ;-- text y position
gsn_text_ndc(wks,title,xl,yb, tires) ;-- plot title string on top
;-------------------------------------------------------------------
;-- do the complete plot
;-------------------------------------------------------------------
draw(plot) ;-- draw the map
frame(wks) ;-- advance the frame
print("--> Number of missing values: "+ num(ismissing(var)))
print("")
;-------------------------------------------------------------------
;-- print elapsed CPU time
;-------------------------------------------------------------------
endt = get_cpu_time()
print("--> Used CPU time: "+ (endt-startt) + "s")
end
Result:
