Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL ICON Mollweide projection + add latitude annotation example
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

DKRZ NCL ICON Mollweide projection + add latitude annotation example

The example uses Mollweide projection and a function to add latitude labels to the plot.

Example script:

;-------------------------------------------------------------------
;-- Script:       ICON_add_polygon_triangles_Mollweide.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
;--               - projection Mollweide
;--               - add latitude label annotation
;--
;-- Notice:       Run time:
;--		    NCL script:		 1.5s
;--		    PyNGL/PyNIO script:	 4.3s
;--		   Matplotlib script: 	28.1s
;--
;-- NCL Version:  6.3.0
;--
;-- 14.03.16      meier-fleischer(at)dkrz.de
;-------------------------------------------------------------------
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
 
;----------------------------------------------------------------------;
;-- Procedure:  add_Mollweide_labels                                 --;
;--    Add latitude labels to plot                                   --;
;----------------------------------------------------------------------;
undef("add_Mollweide_labels")
function add_Mollweide_labels(wks,map,latspc,lonspc)
local lat_values, lat_labels, nlat, txres, x_in, y_in, x_out, y_out
begin

  getvalues map
     "mpProjection"  :  mapproj
  end getvalues
    
  if (mapproj .ne. 15) then
      print("")
      print("Warning: add_Mollweide_labels function only available for Mollweide projection!")
      print("")
      return map
  else
      print("")
      print("--> Add_Mollweide_labels:     add labels for Mollweide projection")
      print("")
  end if
  
  minlat =  -90
  maxlat =   90

;-- pick some "nice" values for the latitude labels
  lat_values = ispan(minlat,maxlat,latspc) * 1.
  nlat       = dimsizes(lat_values)

;-- create the labels; add space before right labels, and after left labels
  lat_labels = where(lat_values.lt.0,abs(lat_values)+"~S~o~N~S",lat_values+"~S~o~N~N")
  lat_labels = where(lat_values.eq.0,"0~S~o~N~",lat_labels)
  lat_labels = where(abs(lat_values).eq.90,"",lat_labels) ;-- no label at lat=abs(90)

;---------------------
;-- latitude labels
;---------------------
  txres                  =  True             ;-- set text resources
  txres@txFontHeightF    =  0.010
  txres@txFontThicknessF =  2

  y_in  = lat_values                         ;-- plot coordinates
  x_in  = new(dimsizes(y_in),float)          ;-- plot coordinates
  x_out = new(dimsizes(x_in),float)          ;-- for NDC coordinates
  y_out = new(dimsizes(y_in),float)          ;-- for NDC coordinates
  
;-- left latitude labels
  txres@txJust = "CenterRight"
  x_in = -180.0
  datatondc(map,x_in,y_in,x_out,y_out)
  gsn_text_ndc(wks, lat_labels, x_out-0.015, y_out, txres)

;-- right latitude labels
  txres@txJust = "CenterLeft"
  x_in = 180.0
  datatondc(map,x_in,y_in,x_out,y_out)
  gsn_text_ndc(wks, lat_labels, x_out+0.015, y_out, txres)

  return(map)
  
end   ;-- end function add_Mollweide_labels

;----------------------------------------------------------------------;
;-- 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.04                ;-- 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:  Variable 'ta'"                  ;-- title string

;  projection = "Robinson"                            ;-- ERROR: edge artifacts
  projection = "Mollweide"
  
  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_"+projection)  ;-- open a workstation
    
;-------------------------------------------------------------------
;-- set resources
;-------------------------------------------------------------------
  res                      =  True
  res@gsnDraw              =  False                   ;-- don't draw the plot
  res@gsnFrame             =  False                   ;-- don't advance the frame
  res@gsnLeftString        = ""
  res@gsnRightString       = ""
  
  res@vpWidthF             =  0.9
  res@vpHeightF            =  0.8
  res@vpXF                 =  0.05
  res@vpYF                 =  0.95
  
  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@mpProjection         =  projection
  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 !!
  res@mpGridAndLimbOn      =  True                    ;-- turn on lat/lon lines
  res@mpPerimOn            =  False                   ;-- turn off box around plot  
  res@mpGridLatSpacingF    =  30.                     ;-- spacing for lat lines
  res@mpGridLonSpacingF    =  30.                     ;-- spacing for lon lines
  
;-------------------------------------------------------------------
;-- 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)) 

  setvalues plot
     "mpPerimOn" :  False                              ;-- get lost using setColorContourClear
  end setvalues
  
;-------------------------------------------------------------------
;-- 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

  print("--> Number of missing values: "+ num(ismissing(var)))
  print("")
  
;-------------------------------------------------------------------
;-- 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.83                                           ;-- text y position
  gsn_text_ndc(wks,title,xl,yb, tires)                ;-- plot title string on top

;-------------------------------------------------------------------
;-- add latitude labels to plot
;-------------------------------------------------------------------
  dlat = 30
  dlon = 30
  map  = add_Mollweide_labels(wks,plot,dlat,dlon)

;-------------------------------------------------------------------
;-- add variable long_name and units between labelbar and plot
;-------------------------------------------------------------------
  stres                =  True
  stres@txFontHeightF  =  0.010
  gsn_text_ndc(wks, var@long_name, 0.223, 0.30, stres)
  gsn_text_ndc(wks, var@units,     0.837, 0.30, stres)
  
;-------------------------------------------------------------------
;-- do the complete plot
;-------------------------------------------------------------------
  draw(plot)        ;-- draw the map
  frame(wks)        ;-- advance the frame

;-------------------------------------------------------------------
;-- print elapsed CPU time
;-------------------------------------------------------------------
  endt = get_cpu_time()
  print("--> Used CPU time:            "+ (endt-startt) + "s")
  
end

Artikelaktionen