DKRZ NCL Spiral of sea ice fraction

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

Result: