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

DKRZ NCL spiral example

The example script shows how to create an animation of multi-year monthly mean temperature change data of the northern hemisphere similar to Ed Hawkins's spiral animation (https://www.climate-lab-book.ac.uk/spirals/).

Example script:

;-----------------------------------------------------------------------
;-- DKRZ NCL example: NCL_spiral_Mean_Temp_Change_RCP85_NH.ncl
;--
;-- Description:      Plot monthly mean temperature change data similar to
;--                      Ed Hawkins's spiral animation:
;--
;--                   https://www.climate-lab-book.ac.uk/spirals/
;--
;--                   Difference of RCP 8.5 data and 20 years global
;--                   monthly mean on the northern hemisphere.
;--
;-- NCL Version:      6.4.0
;--
;--   DKRZ
;--   (German Climate Computing Center)
;--   Bundesstrasse 45a
;--   20146 Hamburg
;--   Germany
;--  
;--   13.07.17  meier-fleischer(at)dkrz.de
;-----------------------------------------------------------------------

;------------------------------------------------------------------------
;-- global variables
;------------------------------------------------------------------------

MIN_VALUE         =  0.0
MAX_VALUE         =  6.0
VALUE_INCREMENT   =  0.5
LABELS_CIRCLE     = (/0.0,1.0,2.0,3.0,4.0,5.0,6.0/)

;------------------------------------------------------------------------

begin

  start_date = toint(systemfunc("date +%s"))            ;-- computing start time

;-- set plot output file name
  rcpname    = "rcp85"                                  ;-- RCP name

  plotdir    = "./Plots_NH_"+rcpname+"/"             ;-- store plot files in directory
  plotfile   = "plot_spiral_Mean_Temp_Change_2005-2100_NH_"+rcpname  ;-- plot file name without extension
  system("mkdir -p "+plotdir)                           ;-- make directory if it doesn't exist
 
  rcptext    = "RCP 8.5"                                ;-- RCP annotation string
  caption    = "northern hemisphere"                    ;-- plot caption

;------------------------------------------------------------------------
;-- compute the multi-year monthly mean and select the northern hemisphere
;------------------------------------------------------------------------
;-- input monthly mean data
  inhis  = "$HOME/data/CMIP5/atmos/hist_mean1-3_LR_temp2_mm_1850-2005.nc"
  inr85  = "$HOME/data/CMIP5/atmos/rcp85_mean1-3_LR_temp2_mm_2005-2100.nc"

;-- create the 20 years monthly means from historical data
  ymonmeanNH = "$HOME/data/CMIP5/atmos/hist_mean1-3_LR_temp2_ymm_1986-2005_NH.nc"
  if(.not. fileexists(ymonmeanNH)) then
     system("cdo -fldmean -ymonmean -sellonlatbox,-180.0,180.0,0.0,90.0 -selyear,1986/2005 "+ inhis +" "+ymonmeanNH)
  end if

;-- compute the field means of the northern hemisphere
  fldmmNH  = "$HOME/data/CMIP5/atmos/rcp85_mean1-3_LR_temp2_fldmm_2005-2100_NH.nc"
  if(.not. fileexists(fldmmNH)) then
     system("cdo -fldmean -sellonlatbox,-180.0,180.0,0.0,90.0  " + inr85 + " "+fldmmNH)
  end if

;------------------------------------------------------------------------
;-- open files and read data
;------------------------------------------------------------------------
  f20NH    =  addfile(ymonmeanNH,"r")
  var20NH  =  f20NH->temp2(:,0,0)                       ;-- temperature data 1st time step 20 years mean

;-- open file and read variable temp2 and time
  frcpNH =  addfile(fldmmNH,"r")
  varNH  =  frcpNH->temp2(:,0,0)                        ;-- temperature data 1st time step rcp85
  time   =  frcpNH->time                                ;-- time values

  utc_date = cd_calendar(time, 0)                       ;-- convert date to UT-referenced date
  years    = sprinti("%0.4i",tointeger(utc_date(:,0)))  ;-- get year as integer value
   
;-- define minimum/maximum and interval of data values
  vrange   = toint(ceil((MAX_VALUE-MIN_VALUE)/VALUE_INCREMENT))          ;-- get smallest integer value larger than input
  levels   = fspan(MIN_VALUE,MAX_VALUE,vrange+1)
  ncircles = dimsizes(levels)

;------------------------------------------------------------------------
;-- graphics
;------------------------------------------------------------------------
;-- set workstation resources
  wks_type                   = "png"                    ;-- plot output type
  wks_type@wkBackgroundColor = "gray18"                 ;-- set workstation background to black (or grey18)
  wks_type@wkForegroundColor = "white"                  ;-- set workstation background to black (or grey18)
  wks_type@wkWidth           =  2500                    ;-- for presentations
  wks_type@wkHeight          =  2500                    ;-- for presentations
  wks = gsn_open_wks(wks_type,plotdir+plotfile)

;-- set resources
  res                      =  True                      ;-- set resources for plot
  res@gsnDraw              =  False                     ;-- don't draw plot yet
  res@gsnFrame             =  False                     ;-- don't advance frame
      
  res@tiMainFontHeightF    =  0.02                      ;-- main title font size
  res@tiMainOffsetYF       =  0.06                      ;-- move title upward
  res@tiMainFontColor      = "white"                    ;-- set font color to white
    
  res@vpWidthF             =  0.8                       ;-- set viewport width
  res@vpHeightF            =  0.8                       ;-- set viewport height
  res@vpXF                 =  0.1                       ;-- set viewport Y position
  res@vpYF                 =  0.9                       ;-- set viewport Y position
    
  res@trXMinF              =  0.                        ;-- x-axis minimum
  res@trXMaxF              =  1.                        ;-- x-axis maximum
  res@trYMinF              =  0.                        ;-- y-axis minimum
  res@trYMaxF              =  1.                        ;-- y-axis maximum
      
  res@tmXBBorderOn         =  False                     ;-- no x-axis bottom line
  res@tmXBOn               =  False                     ;-- don't draw bottom x-axis
  res@tmXTBorderOn         =  False    
  res@tmXTOn               =  False                     ;-- don't draw top x-axis
  res@tmYLBorderOn         =  False    
  res@tmYLOn               =  False                     ;-- don't draw left y-axis
  res@tmYRBorderOn         =  False    
  res@tmYROn               =  False                     ;-- don't draw right y-axis

;-- create blank base plots
  plotNH = gsn_csm_blank_plot(wks,res)

;------------------------------------------------------------------------
;-- set color map and retrieve the colors to change them
;------------------------------------------------------------------------
  gsn_define_colormap(wks,"NCV_jet")                    ;-- define (256 colors) color map
  colors     =  gsn_retrieve_colormap(wks)              ;-- retrieve color map for common labelbar
  ncols      =  dimsizes(colors(:,0))                   ;-- number of colors in color map
  colinc     = (ncols/vrange)                           ;-- color index increment
  col        =  ispan(2,ncols-1,colinc)                 ;-- number of colors we need
  nlines     =  dimsizes(time)                          ;-- number of triangles

  gscolorsNH =  new(nlines,integer)                     ;-- create color array (type of colors integer!)
  gscolorsNH =  col(0)                                  ;-- initialize to black (for missing colors)

;------------------------------------------------------------------------
;-- create differences for color allocation
;-- 20 years monthly means contains only 12 months Jan-Dec
;------------------------------------------------------------------------
  vNH = varNH                                           ;-- copy variable and remain all attributes

  il = 0
  do l=0,dimsizes(vNH)-1
     vNH(l) = varNH(l) - var20NH(il)
     if(il .eq. 11) then
        il = 0
     else
        il = il + 1
     end if
  end do
  il = 0

;-- set color for data less than given minimum value var_min
  vlow = ind(vNH .lt. levels(0))                        ;-- get the indices of values less levels(0)
  if(.not. all(ismissing(vlow))) then
     gscolorsNH(vlow) = col(0)                          ;-- choose color
  end if

;-- set colors for all cells in between var_min and var_max
  do i = 1, dimsizes(levels) - 1
     vind := ind(vNH .ge. levels(i-1) .and. vNH .lt. levels(i))  ;-- get the indices of 'middle' values
     if(.not. all(ismissing(vind))) then
        gscolorsNH(vind) = col(i)                       ;-- choose the colors
     end if
  end do

;-- set color for data greater than given maximum var_max
  nc=dimsizes(col)-1                                    ;-- get the number of colors minus one
  nl=dimsizes(levels)-1                                 ;-- get the number of levels minues one
  vhgh := ind(vNH .gt. levels(nl))                      ;-- get indices of values greater levels(nl)
  if(.not. all(ismissing(vhgh))) then    
     gscolorsNH(vhgh) = col(nc)                         ;-- choose color
  end if    

;-----------------------------------------------------------------------------------
;-- line and text resources
;-----------------------------------------------------------------------------------
;-- set resources for circles and plot circle lines in red
  lnres                   =  True
  lnres@gsLineColor       = "gray40"                ;-- polyline color
  lnres@gsLineThicknessF  =  4.0                    ;-- make lines thicker
  lnres@gsLineDashPattern =  2                      ;-- line dash pattern
 
  pgres                   =  True
  pgres@gsFillColor       = "black"                 ;-- polygon fill color

;-- text resources for value annotation of circles
  txvals                  =  True     
  txvals@txFontColor      = "gray30"                ;-- text color darker gray
  txvals@txFontHeightF    =  0.014                  ;-- make font size smaller
  txvals@txJust           = "CenterRight"           ;-- text justification
  txvals@txFont           =  21                     ;-- text font "courier-bold"
 
;-----------------------------------------------------------------------------------
;-- define viewport values for blank plot    
;-----------------------------------------------------------------------------------
  xf      =  0.05                                       ;-- viewport x position
  yf      =  0.95                                       ;-- viewport y position
  wf      =  0.8                                        ;-- viewport width
  hf      =  0.8                                        ;-- viewport height
  scale   = (wf/2)
 
;-- degrees in radians (pi/180)    
  deg2rad = (4.0*atan(1.))/180.                         ;-- degrees to radians
  degrees =  ispan(0,360,1)                             ;-- assign array 0-360 degrees
  xcos    =  cos(deg2rad * degrees)                     ;-- convert degrees to radians
  xsin    =  sin(deg2rad * degrees)                     ;-- convert to radians
  xcenter =  0.5                                        ;-- x center position
  ycenter =  0.5                                        ;-- y center position

;-- define radii, circles and labels dynamical
  dr       =  scale/(ncircles-1)                   ;-- distance between the circles
  radius   =  True                                  ;-- object contains the radii
  idplNH   =  True                                  ;-- object contains the polylines
  idtxNH   =  True                                  ;-- object contains zero circle

;-- draw the filled black base circle
  radiusblack =  dr * (ncircles)                    ;-- add radius to object
  xblack      =  xcenter + (radiusblack * xcos);-- add x-array to object
  yblack      =  ycenter + (radiusblack * xsin);-- add y-array to object
  blackidNH   =  gsn_add_polygon(wks, plotNH, xblack, yblack, pgres)   ;-- add gray circles to object

  do icic=0,ncircles-1
     rname          = "radius"+icic
     radius@$rname$ =  dr * (icic+1)                    ;-- add radius to object
      
     xname          = "xc"+icic
     radius@$xname$ =  0.5 + (radius@$rname$ * xcos);-- add x-array to object
      
     yname          = "yc"+icic
     radius@$yname$ =  0.5 + (radius@$rname$ * xsin);-- add y-array to object
      
     cicid          =  unique_string("cicid")           ;-- create unique ids for the circles
     idplNH@$cicid$ =  gsn_add_polyline(wks, plotNH, radius@$xname$, radius@$yname$, lnres)   ;-- add gray circles to object

     ilab = ind(LABELS_CIRCLE .eq. levels(icic))

     if(.not. ismissing(ilab)) then
        labcic          =  unique_string("labcic")          ;-- create unique labels of the circles
        idtxNH@$labcic$ =  gsn_add_text(wks,plotNH, toint(LABELS_CIRCLE(ilab)),\
                                        0.5+0.018+((radius@$rname$ + dr)*cos(deg2rad*90.)),\
                                        0.5-0.035+((radius@$rname$ + dr)*sin(deg2rad*90.)),txvals) ;-- add circles to object
     end if
     
     outer_radius   = radius@$rname$
  end do

;-- text resources: labelbar units, copyright and year string
  txtitle                   =  True                 ;-- text resources title string
  txtitle@txFontColor       = "white"               ;-- change to white
  txtitle@txFontHeightF     =  0.022                ;-- text font size
  txtitle@txJust            = "CenterCenter"        ;-- text justification

  txcapt                    =  txtitle              ;-- text resources title string
  txcapt@txFontHeightF      =  0.016                ;-- text font size

  txrcp                     =  txtitle              ;-- text resources rcp string
  txrcp@txFontHeightF       =  0.018                ;-- text font size


  txyear                    =  txtitle              ;-- text resources year string
  txyear@txFontHeightF      =  0.020                ;-- text font size
  txyear@txFont             =  29                   ;-- text font "courier-bold"
 
  txcopy                    =  True                 ;-- text resources copyright string
  txcopy@txFontColor        = "white"               ;-- change to white
  txcopy@txJust             = "BottomRight"         ;-- text justification
  txcopy@txFontHeightF      =  0.008                ;-- make font size smaller
  txcopy@txAngleF           =  -90.0
 
  txmonth                   =  txtitle              ;-- text resources copyright string
  txmonth@txFontHeightF     =  0.008                ;-- make font size smaller

;-- draw the monthly labels around the outer circle
  dtx    = fspan(90.,-270.,13)
  xtx    = xcenter + ((outer_radius+0.03) * cos(deg2rad*dtx))
  ytx    = ycenter + ((outer_radius+0.03) * sin(deg2rad*dtx))
  monlab = (/"JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP","OCT","NOV","DEC"/)
 
  do k=0,11
     txmonth@txAngleF = -30.0 * k
     mstr = unique_string("mtext")
     plotNH@$mstr$ = gsn_add_text(wks,plotNH,monlab(k),xtx(k),ytx(k),txmonth)       ;-- center title string

     plineNH = unique_string("plineNH")
     plotNH@$plineNH$ = gsn_add_polyline(wks, plotNH,(/0.5,xtx(k)/),(/0.5,ytx(k)/), True)
  end do
 
;-----------------------------------------------------------------------------------
;-- loop over time
;-----------------------------------------------------------------------------------
  ntimes = dimsizes(time)                               ;-- number of time steps
  j      = 0                                            ;-- init plot number
  im     = 1                                            ;-- start point
  im1    = 2                                            ;-- next point
  ddtx   = dtx(0:dimsizes(dtx)-2)                       ;-- the 12 months

;-- polyline resource settings
  dres                  =  True
  dres@gsLineColor      =  gscolorsNH(0)                ;-- polyline color
  dres@gsLineThicknessF =  5.0                          ;-- plolyline thickness

;-- do all polylines
  do i=0,ntimes-2
     print("Plot:  "+j+"   year: "+years(i) +"  time step: " + i + " of " + ntimes + " time steps")

;-- compute next data polyline points
     data0NH = ((vNH(i) * scale) / MAX_VALUE) + dr
     xd0NH    = xcenter + (data0NH * cos(deg2rad*dtx(im)))
     yd0NH    = ycenter + (data0NH * sin(deg2rad*dtx(im)))

     data1NH = ((vNH(i+1) * scale) / MAX_VALUE) + dr
     xd1NH    = xcenter + (data1NH * cos(deg2rad*dtx(im1)))
     yd1NH    = ycenter + (data1NH * sin(deg2rad*dtx(im1)))

;-- draw the data polyline segments
     str = unique_string("poly")

     dres@gsLineColor = gscolorsNH(i)
     plotNH@$str$ = gsn_add_polyline(wks,plotNH,(/xd0NH,xd1NH/),(/yd0NH,yd1NH/),dres)         ;-- draw data segment

     if(im .eq. 10) then
       im1 = 0
     else
       im1 = im1 + 1
     end if
     if(im .eq. 11) then
       im = 0
     else
       im = im + 1
     end if
 
     j  = j + 1

;-- add the annotations to the workstation (NDC coordinate system)
     gsn_text_ndc(wks,"Mean Temperature Change  MPI-ESM LR",0.5,0.98,txtitle) ;-- center title string
     gsn_text_ndc(wks,rcptext,  0.2, 0.88, txrcp)                   ;-- plot year string
     gsn_text_ndc(wks,years(i), 0.8, 0.88, txyear)                  ;-- plot year string
     gsn_text_ndc(wks,caption,  0.5, 0.08, txcapt)                  ;-- lower left plot title
     gsn_text_ndc(wks,"~F35~c ~F21~~N~DKRZ / MPI-M", 0.97, 0.15, txcopy) ;-- plot copyright info

;-- advance the frame of each plot
     draw(plotNH)
     frame(wks)
   
  end do            ;-- end loop over all polylines
 
;------------------------------------------------------------------------
;-- print some computing time information
;------------------------------------------------------------------------
  end_date = toint(systemfunc("date +%s"))              ;-- computing used time
  print("")
  print("Start Time:                         "+start_date)
  print("End Time:                           "+end_date)
  print("Time for "+(j)+"  time steps: "+(end_date-start_date)+"s")
  print("")
 
end

Artikelaktionen