DKRZ NCL spiral example

Example script:

;-----------------------------------------------------------------------
;-- DKRZ NCL example: spiral_Mean_Temp_Change_RCP85_NH.ncl
;--
;-- Description: Plot monthly mean temperature change animation 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
;--
;-- 15.06.17 meier-fleischer(at)dkrz.de
;-----------------------------------------------------------------------
;-- function getDate
;--
;-- Read time value and return in this case the year as a string
;-------------------------------------------------------
undef ("getDate")
function getDate(time)
local utc_date, year, mon, day, hours, mins, str_date
begin
  utc_date = cd_calendar(time, 0) ;-- convert date to UT-referenced date
  year = sprinti("%0.4i",tointeger(utc_date(:,0))) ;-- get year as integer value
  mon = sprinti("%0.2i",tointeger(utc_date(:,1))) ;-- get month as integer value
  day = sprinti("%0.2i",tointeger(utc_date(:,2))) ;-- get day as integer value
  hours = sprinti("%0.2i",tointeger(utc_date(:,3))) ;-- get day as integer value
  mins = sprinti("%0.2i",tointeger(utc_date(:,4))) ;-- get day as integer value
  return(year) ;-- return only year of type string
end

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

;-- set plot output file name
 pltvarname = "rcp85" ;-- RCP name
 plotdir = "./Plots_NH_"+pltvarname+"/" ;-- store plot files in directory
 plotfile = "plot_spiral_Mean_Temp_Change_2005-2100_NH_"+pltvarname ;-- 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"
 inr26 = "$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

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

;-- 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 " + inr26 + " "+fldmmNH)
 end if

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

;-- define minimum,maximum and interval of data values
 valMin = 0.0
 valMax = 3.5
 valInt = 0.05
 vrange = toint(ceil((valMax-valMin)/valInt)) ;-- get smallest integer value larger than input
 levels = fspan(valMin,valMax,vrange) ;-- compute number of 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 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

 print("---> number of levels: " + dimsizes(levels))
 print("---> number of colors: " + dimsizes(col))

 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

;-- define viewport values for blank plot
 xf = 0.05 ;-- viewport x position
 yf = 0.95 ;-- viewport y position
 wf = 0.9 ;-- viewport width
 hf = 0.9 ;-- viewport height

;-- 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 = wf/2 + xf ;-- x center position
 ycenter = hf/2 + (yf-hf) ;-- y center position

;-- define the radii
 radius5 = wf/2 ;-- radius of the biggest black filled circle
 dr = (radius5-0.05)/6.1 ;-- distance between the circles
 radius4 = dr*5 ;-- 4.0 degree circle radius
 radius3 = dr*4 ;-- 3.0 degree circle radius
 radius2 = dr*3 ;-- 2.0 degree circle radius
 radius1 = dr*2 ;-- 1.0 degree circle radius
 radius0 = dr ;-- 0.0 degree circle radius

;-- (xcenter + r*cos(n*Pi/180), ycenter + r*sin(n*Pi/180)), n=0,1,2,3,...,359.
;-- 0.0 degree circle
 xc0 = xcenter + (radius0 * xcos)
 yc0 = ycenter + (radius0 * xsin)
;-- 1.0 degree circle
 xc1 = xcenter + (radius1 * xcos)
 yc1 = ycenter + (radius1 * xsin)
;-- 2.0 degree circle
 xc2 = xcenter + (radius2 * xcos)
 yc2 = ycenter + (radius2 * xsin)
;-- 3.0 degree circle
 xc3 = xcenter + (radius3 * xcos)
 yc3 = ycenter + (radius3 * xsin)
;-- 4.0 degree circle
 xc4 = xcenter + (radius4 * xcos)
 yc4 = ycenter + (radius4 * xsin)
;-- big black filled circle
 xc5 = xcenter + (radius5 * xcos)
 yc5 = ycenter + (radius5 * xsin)

;-- 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.80 ;-- set viewport width
 res@vpHeightF = 0.80 ;-- set viewport height
 res@vpXF = 0.10 ;-- set viewport Y position
 res@vpYF = 0.90 ;-- 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)

;-- draw blank base plots
 draw(plotNH)

;-- draw the annotation circles
;-- set resources for circle and plot it with the red circles
 lnres = True
 lnres@gsLineColor = "red"
 lnres@gsLineThicknessF = 4.0
 lnres@gsFillColor = "black"

 blackcircle2 = gsn_add_polygon(wks, plotNH,xc5,yc5,lnres) ;-- filled black circle
 circle42 = gsn_add_polyline(wks,plotNH,xc4,yc4,lnres) ;-- red circle 4.0 deg
 circle32 = gsn_add_polyline(wks,plotNH,xc3,yc3,lnres) ;-- red circle 3.0 deg
 circle22 = gsn_add_polyline(wks,plotNH,xc2,yc2,lnres) ;-- red circle 2.0 deg
 circle12 = gsn_add_polyline(wks,plotNH,xc1,yc1,lnres) ;-- red circle 1.5 deg
 circle02 = gsn_add_polyline(wks,plotNH,xc0,yc0,lnres) ;-- red circle 0.0 deg

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

 txres1 = True ;-- text resources title string
 txres1@txFontColor = "white" ;-- change to white
 txres1@txFontHeightF = 0.020 ;-- text font size
 txres1@txFontHeightF = 0.016 ;-- text font size
 txres1@txJust = "CenterCenter" ;-- text justification

 txres2 = True ;-- text resources year string
 txres2@txFontColor = "white" ;-- change to white
 txres2@txFontHeightF = 0.020 ;-- text font size
 txres2@txJust = "CenterRight" ;-- text justification
 txres2@txFont = 29 ;-- text font "courier-bold"

 txres3 = True ;-- text resources copyright string
 txres3@txFontColor = "white" ;-- change to white
 txres3@txJust = "BottomRight" ;-- text justification
 txres3@txFontHeightF = 0.008 ;-- make font size smaller
 txres3@txAngleF = -90.0

 txres4 = True ;-- text resources copyright string
 txres4@txFontColor = "white" ;-- change to white
 txres4@txJust = "CenterCenter" ;-- text justification
 txres4@txFontHeightF = 0.008 ;-- make font size smaller

 txres5 = True
 txres5@txFontColor = "red"
 txres5@txFontHeightF = 0.009 ;-- make font size smaller
 txres5@txJust = "CenterRight" ;-- text justification
 txres5@txFont = 29 ;-- text font "courier-bold"

;-- monthly labels around plot; xy-coordinates for the monthly labels on circle
 dtx = fspan(90.,-270.,13)
 xtx = xcenter + ((radius5+0.02) * cos(deg2rad*dtx))
 ytx = ycenter + ((radius5+0.02) * sin(deg2rad*dtx))

 monlab = (/"JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP","OCT","NOV","DEC"/)

;-- add month abbreviations to plot
 do k=0,11
    txres4@txAngleF = -30.0 * k ;-- text rotation angle
    mstr = unique_string("mtext") ;-- create unique strings
    plotNH@$mstr$ = gsn_add_text(wks,plotNH,monlab(k),xtx(k),ytx(k),txres4) ;-- center title string
 end do

;-- annotate the red circles with degree values
 s4NH = gsn_add_text(wks,plotNH,"4.0",xcenter+0.028+((radius4)*cos(deg2rad*90.)),xcenter+0.018+((radius4)*sin(deg2rad*90.)),txres5) ;-- draw degrees annotation
 s3NH = gsn_add_text(wks,plotNH,"3.0",xcenter+0.028+((radius3)*cos(deg2rad*90.)),xcenter+0.018+((radius3)*sin(deg2rad*90.)),txres5) ;-- draw degrees annotation
 s2NH = gsn_add_text(wks,plotNH,"2.0",xcenter+0.028+((radius2)*cos(deg2rad*90.)),xcenter+0.018+((radius2)*sin(deg2rad*90.)),txres5) ;-- draw degrees annotation
 s1NH = gsn_add_text(wks,plotNH,"1.0",xcenter+0.028+((radius1)*cos(deg2rad*90.)),xcenter+0.018+((radius1)*sin(deg2rad*90.)),txres5) ;-- draw degrees annotation
 s0NH = gsn_add_text(wks,plotNH,"0.0",xcenter+0.028+((radius0)*cos(deg2rad*90.)),xcenter+0.018+((radius0)*sin(deg2rad*90.)),txres5) ;-- draw degrees annotation

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

;-- compute the start data and its x and y values for the polyline
 data0NH = varNH(0) - var20NH(0)
 xd0NH = 0.5 + ((dr*data0NH)+dr)*cos(deg2rad*ddtx(0))
 yd0NH = 0.5 + ((dr*data0NH)+dr)*sin(deg2rad*ddtx(0))
 data1NH = varNH(1) - var20NH(1)
 xd1NH = 0.5 + ((dr*data1NH)+dr)*cos(deg2rad*ddtx(1))
 yd1NH = 0.5 + ((dr*data1NH)+dr)*sin(deg2rad*ddtx(1))

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

;-- do all polylines

;ntimes = 24 ;-- reduce the number of plots for testing purposes

 do i=1,ntimes-2
    print("Plot: "+j+" year: "+getDate((time(i))) + " month: " + im + " time step: " + i + " of " + ntimes + " time steps")

;-- draw the data polyline segments
    dres@gsLineColor = gscolorsNH(i) ;-- polyline color
    str = unique_string("poly") ;-- define unique string
    plotNH@$str$ = gsn_add_polyline(wks,plotNH,(/xd0NH,xd1NH/),(/yd0NH,yd1NH/),dres) ;-- draw data segment

;-- compute next data polyline points in circle
    data0NH = varNH(i) - var20NH(im)
    xd0NH = 0.5 + ((dr*data0NH)+dr)*cos(deg2rad*ddtx(im))
    yd0NH = 0.5 + ((dr*data0NH)+dr)*sin(deg2rad*ddtx(im))
    data1NH = varNH(i+1) - var20NH(im1)
    xd1NH = 0.5 + ((dr*data1NH)+dr)*cos(deg2rad*ddtx(im1))
    yd1NH = 0.5 + ((dr*data1NH)+dr)*sin(deg2rad*ddtx(im1))

    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

;-- draw the plot
    draw(plotNH)

;-- add the annotations
    gsn_text_ndc(wks,"Mean Temperature Change MPI-ESM LR",0.5,0.95,txres0) ;-- center title string

    gsn_text_ndc(wks,caption, 0.5, 0.08, txres1) ;-- lower left plot title

    txres0@txFontHeightF = 0.018
    gsn_text_ndc(wks,rcptext, 0.15, 0.85, txres0) ;-- plot year string

    gsn_text_ndc(wks,getDate(time(i)),0.90, 0.85, txres2) ;-- plot year string

    gsn_text_ndc(wks,"~F35~c ~F21~~N~DKRZ / MPI-M", 0.90, 0.12, txres3) ;-- plot copyright info

;-- advance the frame of each plot
    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