Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / PyNGL / DKRZ PyNGL example Box-Whisker plot
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

DKRZ PyNGL example Box-Whisker plot

This example provides the function boxplot_mod which creates a Box-Whisker plot. Optional the boxes can be colored and a label can be added to the boxes. The temperature variables of two different model runs are used to compute the statistics and create the plot of the detrended data. The Box-Whisker plot is added to the right of the detrended plot showing two colored boxes.

Example script:

#  File:
#    boxplot_function.py
#
#  Synopsis:
#    Read ASCII data and draw countries of global map with colors depending on their values.
#
#  Category:
#    xy plot
#    statistics
#    polylines, polygons
#    functions
#    text
#
#  Based on DKRZ's NCL example:
#    test_boxplot_mod_data.ncl
#
#  Author:
#    Karin Meier-Fleischer
#  
#  Date of initial publication:
#    December, 2018
#
#  Description:
#    Create a boxplot plot with modufied functions.
#
#  This boxplot_mod function is a based on Dennis Shea's NCL boxplot function from
#  the shea_util.ncl library.
#
#  Added resources to boxOpts:
#
#      boxWidth         --> change the width of the boxes, default: 1.0
#      fillColorOn      --> turn on/off fill color
#      fillColors       --> array of colors ['green','blue']; default 'gray70'
#      labelPosition    --> set position: 'right' or 'left'
#      boxLabels        --> array of box labels e.g. ["tas 1","tas 2"]
#      LabelFontHeightF --> change the box label font size, default 0.010
#
#  Effects illustrated:
#    o  Create an xy-plot
#    o  compute statistics
#    o  Define function boxplot_mod
#    o  Create color filled boxplot
#    o  Add polylines and polygons
#    o  Add text
#
#
#  Output:
#     One visualization is produced.     
#
'''
  PyNGL Example:     boxplot_function.py

  -  Create an xy-plot
  -  Compute statistics
  -  Define function boxplot_mod
  -  Create color filled boxplot
  -  Add polylines and polygons
  -  Add text

'''
from __future__ import print_function
import sys, os
import netCDF4
import numpy as np
from scipy import stats,signal
import Ngl,Nio

#---------------------------------------------------------------
# Function boxplot_mod
#
#            Create the Box-Whisker plot from the given statistics.
#---------------------------------------------------------------
def boxplot_mod(wks,y,boxOpts,plotres,lineres):

#-- set box label font size
   if(hasattr(boxOpts,'LabelFontHeightF')):
      textsize = boxOpts.LabelFontHeightF
   else:
#      boxOpts.LabelFontHeightF = 0.010
      textsize = 0.010
   
#-- set box labels
   if(hasattr(boxOpts,'boxLabels')):
      x = np.arange(0,len(boxOpts.boxLabels),1)
   else:
      boxOpts.boxLabels = ['default' for x in range(0,y.shape[0])]
      
#-- set data array, number of boxes to be drawn and their width
   y         = np.array(y)
   numbox    = len(boxOpts.boxLabels)
   boxWidths = np.zeros(numbox,dtype='float')

   xAxis     = np.arange(0,numbox+2,1)              #-- x-axis values
   dx        = len(xAxis)/numbox
   
   labarr = xAxis.astype('str')                     #-- define x-axis labels
   labarr[0] = ''
   labarr[numbox+1] = ''
   for i in range(0,len(plotres.tmXBLabels)):
      labarr[i+1] = plotres.tmXBLabels[i]

#-- set box width
   if(boxOpts):
      if(hasattr(boxOpts,'boxWidth')):
         if(len(boxOpts.boxWidth) != 1 and len(boxOpts.boxWidth) != numbox):
            print("boxplot: Number of input box widths must either equal 1 or the number of boxes (%4i). Using first specified box width only." % numbox)
            boxWidths[:] = boxOpts.boxWidth[0]
         else:
            boxWidths = boxOpts.boxWidth
      else:
         boxWidths[:] = dx*.2
   else:
      boxWidths[:] = dx*.2
   
#-- create a blank plot   
   plotres.trXMinF    =  xAxis.min()
   plotres.trXMaxF    =  xAxis.max()
   plotres.tmXBMode   = 'Explicit'                  #-- set x axis tickmark labels
   plotres.tmXBLabels =  list(labarr)
   plotres.tmXBValues =  xAxis
   
   plot = Ngl.blank_plot(wks,plotres)

#-- set resources for the polylines
   polyres                   =  Ngl.Resources()     #-- set up defaults
   polyres.gsLineColor       = 'black'              #-- color of lines
   polyres.gsLineThicknessF  =  1.5                 #-- thickness of lines
   polyres.gsLineDashPattern =  0

   fsatts = list(lineres.__dict__.keys())           #-- retrieve the resource attributes
   
   if(lineres):
      for ty in range(0,len(fsatts)):
         elements = getattr(lineres,fsatts[ty])
         setattr(polyres,fsatts[ty],elements)
      
   if(boxOpts):
      if(not hasattr(boxOpts,'fillColorOn')):
         boxOpts.fillColorOn = False

      if(hasattr(boxOpts,'boxColors')):
         boxcolor = boxOpts.boxColors
         if(len(boxcolor) == 1 or len(boxcolor) != numbox):
            if(len(boxcolor) != numbox):
               print("boxplot: warning: Number of input colors must either equal 1 or the number of boxes (%2i). Using first specified color only." % numbox)
            polyres.gsLineColor = boxcolor[0]
            cflag = 1
         else:
            cflag = 2
      else:
         cflag = 1
   else:
      cflag = 1         
      
#-- set label text resources
   textres                    =  Ngl.Resources()
   textres.txFontColor        = 'black'             #-- set text color to black
   textres.txFontHeightF      =  textsize           #-- decrease font size
   textres.txJust             = 'CenterCenter'
   textres.txAngleF           =  90.                #-- rotate text counterclockwise 90 deg.

#-- assign graphic array for polylines (outlines); add graphic arrays for filled polygons and labels
   dum   = []
   dumf  = []
   dumtx = []

#-- draw filled boxes and label them if turned on  
   for im in range(0,numbox):
      ff = xAxis[im+1]      
      if (hasattr(boxOpts,'fillColorOn')):
         if(boxOpts.fillColorOn == True):
            if(hasattr(boxOpts,'fillColors')):
               polyres.gsLineColor = 'black'        #-- set to black for the outlines
               polyres.gsFillColor =  boxOpts.fillColors[im]
               cflag = 1
            else:
               polyres.gsFillColor =  'gray70'      #-- default fill color
               
            pxl = (ff-(boxWidths[im]/2.))
            pxr = (ff+(boxWidths[im]/2.))
            px  = [pxl,pxr,pxr,pxl,pxl]
            py  = [y[im,3],y[im,3],y[im,1],y[im,1],y[im,3]]
            dumf.append(Ngl.add_polygon(wks,plot,px,py,polyres))
   
      if(hasattr(boxOpts,'boxLabels')):
         if(hasattr(boxOpts,'labelPosition')):
            if(boxOpts.labelPosition == 'right'):
               txx = (ff+(boxWidths[im]/2.)) + (boxWidths[im]/2.) + 0.03
            else:
               txx = (ff-(boxWidths[im]/2.)) - (boxWidths[im]/2.) - 0.03
         else:
            txx = (ff+(boxWidths[im]/2.)) + (boxWidths[im]/2.) + 0.03 #-- default position 'right'
            
         label = boxOpts.boxLabels[im]
         txy   = y[im,3] - (y[im,3]-y[im,1])/2.
         dumtx.append(Ngl.add_text(wks,plot,label,txx,txy,textres))
   
#-- create the polylines for the box-whisker plot
   for gg in range(0,numbox):
      ff = xAxis[gg+1]
      
      if(cflag == 2):
         polyres.gsLineColor = boxcolor[gg]

#-- y(n,0) = bottom_value,
#-- y(n,1) = bottom_value_of_box,
#-- y(n,2) = mid-value_of_box,
#-- y(n,3) = top_value_of_box,
#-- y(n,4) = top_value.

#-- short top horizontal line
      yy = np.array([y[gg,4],y[gg,4]])
      xx = np.array([(ff-(boxWidths[gg]/4.)),(ff+(boxWidths[gg]/4.))])
      if(not((np.isnan(xx).any() or np.isnan(yy).any()))):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))
#-- upper vertical line
      yy = [y[gg,3],y[gg,4]]
      xx = [ff,ff]
      polyres.gsLineDashPattern = 1
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))

#-- upper horizontal box line
      yy = [y[gg,3],y[gg,3]]
      xx = [(ff-(boxWidths[gg]/2.)),(ff+(boxWidths[gg]/2.))]
      polyres.gsLineDashPattern = 0
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))
#-- lower horizontal box line
      yy = [y[gg,1],y[gg,3]]
      xx = [(ff-(boxWidths[gg]/2.)),(ff-(boxWidths[gg]/2.))]
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))
#-- median line of box
      yy = [y[gg,2],y[gg,2]]
      xx = [(ff-(boxWidths[gg]/2.)),(ff+(boxWidths[gg]/2.))]
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))
#-- right vertical line of box
      yy = [y[gg,1],y[gg,3]]
      xx = [(ff+(boxWidths[gg]/2.)),(ff+(boxWidths[gg]/2.))]
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))
#-- left vertical line of box
      yy = [y[gg,1],y[gg,1]]
      xx = [(ff-(boxWidths[gg]/2.)),(ff+(boxWidths[gg]/2.))]
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))

#-- lower vertical line
      yy = [y[gg,0],y[gg,1]]
      xx = [ff,ff]
      polyres.gsLineDashPattern = 1
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))
#-- lower horizontal line
      yy = [y[gg,0],y[gg,0]]
      xx = [(ff-(boxWidths[gg]/4.)),(ff+(boxWidths[gg]/4.))]
      polyres.gsLineDashPattern = 0
      if(not (np.isnan(xx).any() or np.isnan(yy).any())):
         dum.append(Ngl.add_polyline(wks,plot,xx,yy,polyres))

   return(plot)


#-------------------------------------------------------
#--              MAIN
#-------------------------------------------------------
#-- data path and file names
home  = os.environ.get('HOME')
fname1 = home + '/data/CORDEX/tas_AFR-44_HadGEM2-ES_rcp45_r1i1p1_CCLM_4-8-17_ym_20060101-20981231-tser.nc'
fname2 = home + '/data/CORDEX/tas_AFR-44_MPI-ESM-LR_rcp45_r1i1p1_CCLM_4-8-17_ym_20060101-20981231-tser.nc'

#-- open files
f1   = Nio.open_file(fname1,"r")                #-- open data file
f2   = Nio.open_file(fname2,"r")                #-- open data file

#-- read time and convert to UTC date
t        = f1.variables['time']
attrs    = t.attributes
units    = attrs['units']
calendar = attrs['calendar']
utc_date = netCDF4.num2date(t, units, calendar)

year     = [date.year for date in utc_date]
month    = [date.month for date in utc_date]
day      = [date.day for date in utc_date]
hour     = [date.hour for date in utc_date]
minute   = [date.minute for date in utc_date]
year     = np.array(year)
tmin     = min(year)
tmax     = max(year)

#-- select variables of timesteps
t_ind  = Ngl.ind(year >= tmin)
nyear  = len(t_ind)
ntimes = np.arange(0,nyear-1,1)
tas1   = f1.variables['tas'][t_ind[0]:t_ind[-1],0,0,0]
tas2   = f2.variables['tas'][t_ind[0]:t_ind[-1],0,0,0]

#-- detrend data
slope1,intercept1,r_value1,p_value1,std_err1 = stats.linregress(ntimes,tas1)
trend1     = intercept1 + (slope1*ntimes)
detrended1 = signal.detrend(tas1) + intercept1

slope2,intercept2,r_value2,p_value2,std_err2 = stats.linregress(ntimes,tas2)
trend2     = intercept2 + (slope2*ntimes)
detrended2 = signal.detrend(tas2) + intercept2

#-- assign array to hold the basic statistics list for Box-Whisker plot
#--    yb(n,0)=bottom_value,
#--    yb(n,1)=bottom_value_of_box,
#--    yb(n,2)=mid-value_of_box,
#--    yb(n,3)=top_value_of_box,
#--    yb(n,4)=top_value.
yb = np.ndarray([2,5],dtype='float')

#-- retrieve statistics data for first variable
yb[0,0] = np.amin(detrended1)                           #-- min
yb[0,1] = np.percentile(detrended1, 25)                 #-- low percentile (25%)
yb[0,2] = np.median(detrended1)                         #-- median
yb[0,3] = np.percentile(detrended1, 75)                 #-- high percentile (75%)
yb[0,4] = np.amax(detrended1)                           #-- max                       
   
#-- retrieve statistics data for second variable   
yb[1,0] = np.amin(detrended2)                           #-- min
yb[1,1] = np.percentile(detrended2, 25)                 #-- low percentile (25%)
yb[1,2] = np.median(detrended2)                         #-- median
yb[1,3] = np.percentile(detrended2, 75)                 #-- high percentile (75%)
yb[1,4] = np.amax(detrended2)                           #-- max

#-- open graphic output
wks =  Ngl.open_wks('png','plot_boxplot_mod')

#-- font size of titles and axis labels
font_size = 0.012
 
#-- timeseries resources - xy-plot
res                               =  Ngl.Resources()
res.nglDraw                       =  False              #-- don't draw plot
res.nglFrame                      =  False              #-- don't advance the frame
res.nglMaximize                   =  False              #-- don't maximize the plot
               
res.vpXF                          =  0.08               #-- set plot x-position
res.vpYF                          =  0.70               #-- set plot y-position
res.vpWidthF                      =  0.62               #-- set plot width
res.vpHeightF                     =  0.30               #-- set plot height
              
res.tiXAxisString                 = "Year"              #-- x-axis titles
res.tiYAxisString                 = "temperature [degC]"#-- y-axis titles
res.tiXAxisFontHeightF            =  font_size
res.tiYAxisFontHeightF            =  font_size
        
res.tmXBLabelFontHeightF          =  font_size
res.tmYLLabelFontHeightF          =  font_size - 0.002

res.tmXMajorGrid                  =  True               #-- draw x-grid lines
res.tmXMajorGridThicknessF        =  1.0
res.tmXMajorGridLineDashPattern   =  2

res.tmYMajorGrid                  =  True               #-- draw y-grid lines
res.tmYMajorGridThicknessF        =  1.0
res.tmYMajorGridLineDashPattern   =  2

res.trXMinF                       =  tmin               #-- x-axis minimum
res.trXMaxF                       =  tmax               #-- x-axis maximum
              
res.trYMinF                       =  min(min(detrended1),min(detrended2))-2 #-- x-axis minimum
res.trYMaxF                       =  max(max(detrended1),max(detrended2))+2 #-- y-axis maximum

res.xyLineThicknessF              =  3.0                #-- line width
res.xyDashPattern                 =  0                  #-- 0: solid line

#-------------------------------------------------------------------
#-- left plot: detrended time series
#-------------------------------------------------------------------
res.xyLineColor = "darkgreen"
plot0 = Ngl.xy(wks, year[t_ind[0]:t_ind[-1]], detrended1, res)  #-- variable 1

res.xyLineColor = "red"
plot00 = Ngl.xy(wks, year[t_ind[0]:t_ind[-1]], detrended2, res) #-- variable 2

#-- overlay plot00 on plot0 and draw plot0
Ngl.overlay(plot0,plot00)
Ngl.draw(plot0)

#-------------------------------------------------------------------
#-- right plot: Box-Whisker plot
#-------------------------------------------------------------------
vpx = Ngl.get_float(plot0,"vpXF")                       #-- retrieve plot x-position
vpy = Ngl.get_float(plot0,"vpYF")                       #-- retrieve plot y-position
vpw = Ngl.get_float(plot0,"vpWidthF")                   #-- retrieve plot width
vph = Ngl.get_float(plot0,"vpHeightF")                  #-- retrieve plot height

#-- resources for the second plot
resBW                       =  Ngl.Resources()
resBW.vpXF                  =  vpx+vpw+0.08             #-- set plot position
resBW.vpWidthF              =  0.2                      #-- set plot width
resBW.vpHeightF             =  0.6                      #-- set plot height
resBW.vpYF                  =  vpy-(vph/2)+(resBW.vpHeightF/2) #-- set plot position

resBW.trYMinF               =  np.floor(yb.min())       #-- y-axis minimum value
resBW.trYMaxF               =  np.ceil(yb.max())        #-- y-axis maximum value

resBW.tmXTOn                =  False                    #-- don't draw top x-axis
resBW.tmXTBorderOn          =  False
resBW.tmYROn                =  False                    #-- don't draw right y-axis
resBW.tmYRBorderOn          =  False

resBW.tiMainString          = "Box-Whisker Plot ~C~of detrended tas"
resBW.tiMainFontHeightF     =  font_size

resBW.tmXBLabels            = ["tas 1","tas 2"]         #-- set labels for boxes

resBW.tmXBLabelFontHeightF  =  font_size
resBW.tmYLLabelFontHeightF  =  font_size

#-- line resources
lresBW                      =  Ngl.Resources()            
lresBW.gsLineThicknessF     =  2.5                      #-- polyline width

#-- special resources for the Box-Whisker plot
OptBW                       =  Ngl.Resources()            
OptBW.fillColorOn           =  True                     #-- turn fill color on, default False
OptBW.boxLabels             = ["HadGEM2","MPI-ESM"]     #-- set box labels, default 'default'
OptBW.fillColors            =  ["darkgreen","red"]      #-- set fill colors, default gray70
OptBW.boxWidth              =  [.2,.2]                  #-- set box width, default 1.0
OptBW.labelPosition         = "right"                   #-- label position "right" or "left", default right
OptBW.LabelFontHeightF      =  0.009                    #-- box label font size, default 0.01
              
#-- create the right plot
plot1 = boxplot_mod(wks, yb, OptBW, resBW, lresBW)

#-- draw plot1 and advance the frame
Ngl.draw(plot1)

#-- write title to right plot
txres               =  Ngl.Resources()
txres.txJust        = 'CenterLeft'
txres.txFontHeightF =  font_size

Ngl.text_ndc(wks,"Detrended tas from "+str(tmin)+"-"+str(tmax),vpx,vpy +0.03,txres)

#-- advance the frame
Ngl.frame(wks)

Result:

Artikelaktionen