Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL ICON add_polygon triangles example
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

DKRZ NCL ICON add_polygon triangles example

The script creates a 1024x1024 pixel PNG plot of the ICON model data on a global map. The traingles are drawn using the add_polygon function of NCL.

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

Artikelaktionen