You are here: Home / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL Spiral of sea ice fraction

DKRZ NCL Spiral of sea ice fraction

This special plot shows the sea ice fraction (yearly mean) for the Arctic reagion as a spiral plot along time based on the plot 'Arctic Death Spiral' by Andy Lee Robinson, 2017. NCL doen't provide this type of plot so we have to build it from scratch.

Example script:

;------------------------------------------------------------------------
;-- DKRZ NCL example: seaIce_spiral_connect_end_start.ncl
;--
;-- Description:      Plot yearly seaIce fraction similar to the spiral
;--                   plot 'Arctic Death Spiral' of Andy Lee Robinson
;--                   (2017)
;--
;-- NCL Version:      6.4.0
;--
;--   DKRZ
;--   (German Climate Computing Center)
;--   Bundesstrasse 45a
;--   20146 Hamburg
;--   Germany
;--  
;--   14.07.17  meier-fleischer(at)dkrz.de
;------------------------------------------------------------------------
;-- global variables
;------------------------------------------------------------------------
NUM_OF_YEARS     =  50          ;-- how many years should be drawn
MIN_VALUE        =   0          ;-- minimum value
MAX_VALUE        =  25          ;-- maximum value (outer circle)
VALUE_INCREMENT  =   5          ;-- distance value of circles

CONNECT          =  True        ;-- True: connect last and first, False: don't connect

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

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

  diri     = "$HOME/data/CMIP5/seaIce/"
  icefile1 = "sic_OImon_MPI-ESM-LR_rcp45_r1i1p1_200601-210012.nc"
  icefile2 = "sic_OImon_MPI-ESM-LR_rcp45_r1i1p1_210101-230012.nc"
  icefile  = "sic_OImon_MPI-ESM-LR_rcp45_r1i1p1_200601-230012.nc"
  ymfile   = "sic_OImon_MPI-ESM-LR_rcp45_r1i1p1_200612-230012_arctic_fldmean_"+NUM_OF_YEARS+"y.nc"
 
  varname  = "sic"
 
  title    = "Sea Ice Fraction ~Z70~(yearly mean)"
  title2   = "~Z55~MPI-ESM LR"                          ;-- title string, top
  subtitle = "Arctic ~Z75~(lat: 50-90~S~o~N~)"          ;-- sub-title string, bottom center
  rcptext  = "RCP 4.5"                                  ;-- rcp string, upper left

;-- merge all data files
  if(.not. fileexists(diri+icefile)) then
     system("cdo -r -mergetime "+diri+icefile1+" "+diri+icefile2+" "+diri+icefile)
  else
     print("-->  sic data file 200612-230012 already exist")
  end if

;-- select 40 years arctic region (lat: 50-90)
  if(.not. fileexists(diri+ymfile)) then
     NUM_OF_MONTHS = NUM_OF_YEARS*12
     system("cdo -r --reduce_dim -fldmean -sellonlatbox,-180.0,180.0,50.0,90.0 -seltimestep,1/"+NUM_OF_MONTHS+" "+diri+icefile+" "+diri+ymfile)
  else
     print("-->  sic "+NUM_OF_YEARS+" years field mean arctic already exist")
  end if

;------------------------------------------------------------------------
;-- open the data file and read variable time
;------------------------------------------------------------------------
  f      = addfile(diri+ymfile,"r")
  sic    = f->sic
  time   = f->time
  ntimes = dimsizes(time)

;------------------------------------------------------------------------
;-- select the data of each month
;------------------------------------------------------------------------
  sic_jan = sic( 0:ntimes-1:12)
  sic_feb = sic( 1:ntimes-1:12)
  sic_mar = sic( 2:ntimes-1:12)
  sic_apr = sic( 3:ntimes-1:12)
  sic_may = sic( 4:ntimes-1:12)
  sic_jun = sic( 5:ntimes-1:12)
  sic_jul = sic( 6:ntimes-1:12)
  sic_aug = sic( 7:ntimes-1:12)
  sic_sep = sic( 8:ntimes-1:12)
  sic_oct = sic( 9:ntimes-1:12)
  sic_nov = sic(10:ntimes-1:12)
  sic_dec = sic(11:ntimes-1:12)

  sic_all = (/sic_jan,sic_feb,sic_mar,sic_apr,sic_may,sic_jun,sic_jul,sic_aug,sic_sep,sic_oct,sic_nov,sic_dec/)

;-- minimum/maximum/increment of data
  valMin    =  MIN_VALUE
  valMax    =  MAX_VALUE
  valInt    =  VALUE_INCREMENT
  circleinc =  valInt               ;-- circle increment
  ncircles  =  valMax/circleinc     ;-- number of circles

;-- define 12 colors
  colors    = (/"red", "orange", "gold", "limegreen", "green","lightblue", \
                "cyan","blue","black","purple","pink","magenta"/)
  ncols     =  dimsizes(colors)

;-- read years
  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
  years    :=  years(::12)             ;-- get years only once
  nyears    =  dimsizes(years)
  print("")
  print("--> Number of years:  "+nyears)
  print("--> Number of colors: "+ncols)
  print("")
 
;------------------------------------------------------------------------
;-- graphics
;------------------------------------------------------------------------
;-- set workstation resources

  window_width  = 1200
  window_height = 1200
  if(window_width .eq. 1200) then
     line_width  =  8
  else
     line_width  = 3
  end if
  end if
 
  wks_type                  = "png"                     ;-- plot output type
  wks_type@wkWidth          =  window_width             ;-- for presentations
  wks_type@wkHeight         =  window_height            ;-- for presentations
  wks = gsn_open_wks(wks_type,"plot_seaIce_arctic_spiral_connect")

;-- define viewport values for blank plot    
  xf    =  0.01                                         ;-- viewport x position
  yf    =  0.95                                         ;-- viewport y position
  wf    =  0.88                                         ;-- viewport width
  hf    =  0.88                                         ;-- viewport height
  scale =  wf/2
 
;------------------------------------------------------------------------
;-- set resources for the blank plot
;------------------------------------------------------------------------
  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.07                        ;-- set viewport Y position
  res@vpYF               =  0.87                         ;-- 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
  plot = gsn_csm_blank_plot(wks,res)

;-----------------------------------------------------------------------------------
;-- draw the circles and labels
;-----------------------------------------------------------------------------------
;-- set resources for circles
  lnres                  =  True
  lnres@gsLineColor      = "gray"
  lnres@gsLineThicknessF =  4.0

;-- text resources for value annotation of circles
  txvals                 =  True     
  txvals@txFontColor     = "black"     
  txvals@txFontHeightF   =  0.014                       ;-- make font size smaller
  txvals@txJust          = "CenterRight"                ;-- text justification
  txvals@txFont          =  21                          ;-- text font "courier-bold"

;-- 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
  ycenter =  0.5
 
;------------------------------------------------------------------------
;-- define the radii, circles and labels dynamical    
;------------------------------------------------------------------------
  dr      =  scale/(ncircles)                    ;-- distance between the circles
  radius  =  True                                       ;-- object contains the radii
  idpl    =  True                                       ;-- object contains the polylines
  idtx    =  True                                       ;-- object contains zero circle
 
  do icic=0,ncircles-1
     rname          = "radius"+icic
     radius@$rname$ =  dr * (icic+1)                    ;-- add radius to object
     
     xname          = "xc"+icic
     radius@$xname$ =  xcenter + (radius@$rname$ * xcos);-- add x-array to object
     
     yname          = "yc"+icic
     radius@$yname$ =  ycenter + (radius@$rname$ * xsin);-- add y-array to object
     
     cicid          =  unique_string("cicid")           ;-- create unique ids for the circles
     idpl@$cicid$   =  gsn_add_polyline(wks, plot,radius@$xname$,radius@$yname$,lnres)   ;-- add gray circles to object
     
     labcic         =  unique_string("labcic")          ;-- create unique labels of the circles
     idtx@$labcic$  =  gsn_add_text(wks,plot,(toint(icic)+1)*circleinc,\
                                    xcenter-0.01+((radius@$rname$)*cos(deg2rad*90.)),\
                                    xcenter     +((radius@$rname$)*sin(deg2rad*90.)),\
                                    txvals)             ;-- add circles to object
     if(icic .eq. 0) then
        idtx@zerocirc =  gsn_add_text(wks,plot,toint(icic),xcenter-0.01,xcenter,txvals)             ;-- add zero circle 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@txFontHeightF  =  0.022                       ;-- text font size
  txtitle@txJust         = "CenterCenter"               ;-- text justification
  txtitle@txFont         =  22                          ;-- text font "courier-bold"
    
  txrcp                  =  txtitle                     ;-- text resource RCP string
  txrcp@txFontHeightF    =  0.015                       ;-- text font size
   
  txsubt                 =  txtitle                     ;-- text resources subtitle string
  txsubt@txFontHeightF   =  0.018                       ;-- text font size
  txsubt@txJust          = "CenterLeft"                 ;-- text justification
      
  txsubt1                =  txsubt                      ;-- text resources subtitle string
  txsubt1@txFontHeightF  =  0.015                       ;-- text font size
      
  txcopy                 =  txtitle                     ;-- text resources copyright string
  txcopy@txFontHeightF   =  0.008                       ;-- make font size smaller
  txcopy@txJust          = "CenterRight"                ;-- text justification

  txyears                =  txtitle                     ;-- text resources years string
  txyears@txFont         =  21                          ;-- text font "courier-bold"
  txyears@txFontColor    = "blue"

;-- resize the font size depending on number of years
  if(NUM_OF_YEARS .ge. 70) then
     txyears@txFontHeightF = 0.009
  else if(NUM_OF_YEARS .ge. 60) then
     txyears@txFontHeightF = 0.008
  else if(NUM_OF_YEARS .ge. 50) then
     txyears@txFontHeightF = 0.009
  else if(NUM_OF_YEARS .ge. 40) then
     txyears@txFontHeightF = 0.010
  else if(NUM_OF_YEARS .lt. 40) then
     txyears@txFontHeightF = 0.012     
  else
     txyears@txFontHeightF = 0.012
  end if
  end if
  end if
  end if
  end if

;-- yearly labels around plot; xy-coordinates for the yearly labels on circle
  dtx = fspan(90.,-270.,nyears+3)
 
  xtx := xcenter + ((outer_radius+0.05) * cos(deg2rad*dtx)) ;-- x-array for years labels
  ytx := ycenter + ((outer_radius+0.05) * sin(deg2rad*dtx)) ;-- y-array for years labels
 
  xlx := xcenter + ((outer_radius+0.02) * cos(deg2rad*dtx)) ;-- x-array for the lines from center to outer circle
  ylx := xcenter + ((outer_radius+0.02) * sin(deg2rad*dtx)) ;-- y-array for the lines from center to outer circle
 
  do k=0,nyears-1
     if(NUM_OF_YEARS .ge. 70) then
        txyears@txAngleF = dtx(k)                       ;-- rotate the year string
     end if
     ystr = unique_string("nyears")                     ;-- create unique strings
     plot@$ystr$ = gsn_add_text(wks,plot,years(k),xtx(k),ytx(k),txyears)  ;-- years string
  end do

;-- add lines from center to max radius to plot
  dres                  =  True
  dres@gsLineColor      =  (/0.7,0.7,0.7,1.0/)          ;-- polyline color
  dres@gsLineThicknessF =  5.0                          ;-- plolyline thickness
 
  do k=0,nyears-1
     lstr = unique_string("nline")                      ;-- create unique strings
     plot@$lstr$ = gsn_add_polyline(wks,plot,(/xcenter,xlx(k)/),(/ycenter,ylx(k)/),dres)  ;-- center title string
  end do

;-----------------------------------------------------------------------------------
;-- data polyline resource settings
;-----------------------------------------------------------------------------------
  vres                  =  True
  vres@gsLineThicknessF =  line_width                       ;-- plolyline thickness
 
;-- create all data polylines
  do j=0,11
     sic_month = sic_all(j,:)

     do i=0,nyears-2
        print("Plot:  month "+j+"   year: "+years(i))

        data0  = ((sic_month(i) * scale) / MAX_VALUE)
        xd0    =  xcenter + (data0 * cos(deg2rad*dtx(i)))
        yd0    =  ycenter + (data0 * sin(deg2rad*dtx(i)))
   
        data1  = ((sic_month(i+1) * scale) / MAX_VALUE)
        xd1    =  xcenter + (data1 * cos(deg2rad*dtx(i+1)))
        yd1    =  ycenter + (data1 * sin(deg2rad*dtx(i+1)))

        vres@gsLineColor = colors(j)                        ;-- polyline color
        vstr = unique_string("vline")                       ;-- define unique string
        plot@$vstr$ = gsn_add_polyline(wks,plot,(/xd0,xd1/),(/yd0,yd1/),vres)   ;-- draw data segment

        ;-- connect last year with first year
        if(i .eq. 0) then
           firstx0 = xd0
           firsty0 = yd0
        end if
        
        if(CONNECT .and. (i .eq. (nyears-2))) then
           lastx1  =  xd1
           lasty1  =  yd1
           vres1   =  True
           vres1@gsLineColor        = "gray60"              ;-- polyline color
           vres1@gsLineThicknessF   =  line_width-2        ;-- plolyline thickness
           vres1@gsLineDashPattern  =  2                    ;-- line type
           vstr1 = unique_string("vline1")                  ;-- define unique string
           plot@$vstr1$ = gsn_add_polyline(wks,plot,(/firstx0,lastx1/),(/firsty0,lasty1/),vres1)   ;-- draw data segment
        end if
        
     end do      ;-- end loop over all years
  end do      ;-- end loop over all months

;-----------------------------------------------------------------------------------
;-- add the annotations
;-----------------------------------------------------------------------------------
  lx = 0.02
  gsn_text_ndc(wks,title,    0.5-lx,  0.97, txtitle)        ;-- top title string
  gsn_text_ndc(wks,title2,   0.5-lx,  0.93, txtitle)        ;-- top title string
  gsn_text_ndc(wks,subtitle, 0.13-lx,  0.88, txsubt)         ;-- subtitle
  gsn_text_ndc(wks,rcptext,  0.82-lx,  0.88, txrcp)          ;-- plot RCP string

;-----------------------------------------------------------------------------------
;-- create legend from scratch
;-----------------------------------------------------------------------------------
  lineres                  =  True
  lineres@gsLineThicknessF =  line_width+5                  ;-- line thickness

  labelres                 =  True
  labelres@txFontHeightF   =  0.012                         ;-- label font height
  labelres@txJust          = "CenterLeft"                   ;-- label font justification

  lgx  = 0.90
  lgy  = 0.59
  lgdx = 0.03
  lgdy = 0.02
 
  monlabels = (/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)

  do ii=0,11
     lineres@gsLineColor  = colors(ii)
     gsn_polyline_ndc(wks, (/lgx,lgx+lgdx/), (/lgy-(lgdy*ii),lgy-(lgdy*ii)/), lineres)
     gsn_text_ndc(wks, monlabels(ii), lgx+lgdx+0.02, lgy-(lgdy*ii), labelres)
  end do

  lgy_down = lgy-(lgdy*11)-0.01
  xn = 0.09
  yn = 0.1
  lgbox_x = (/lgx-0.01, lgx+xn,  lgx+xn,  lgx-0.01,lgx-0.01/)
  lgbox_y = (/lgy+0.01, lgy+0.01, lgy_down, lgy_down,lgy+0.01/)
 
  lineres@gsLineThicknessF =  2                             ;-- line thickness
  lineres@gsLineColor      = "black"                        ;-- line color

  gsn_polyline_ndc(wks, lgbox_x, lgbox_y, lineres)

;-----------------------------------------------------------------------------------
;-- add copyright texts
;-----------------------------------------------------------------------------------
  copyright = "~F35~c ~F21~~N~DKRZ / MPI-M"
  gsn_text_ndc(wks, copyright, 0.93, 0.105, txcopy)         ;-- plot copyright info

  txcopy@txFontHeightF   =  0.007                           ;-- make font size smaller
  txcopy@txFontColor     = "gray20"                         ;-- font color

  original = "~F35~c ~F21~~N~based on 'Arctic Death Spiral',Andy Lee Robinson,2017"
  gsn_text_ndc(wks, original, 0.93, 0.085, txcopy)          ;-- plot copyright info


;-----------------------------------------------------------------------------------
;-- draw the plot and advance the frame
;-----------------------------------------------------------------------------------
  draw(plot)
  frame(wks)

;------------------------------------------------------------------------
;-- 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

Document Actions