DKRZ NCL example relative difference plotΒΆ

Left plot:

reldiff = (v1-v2)/max(v1,v2)

Right plot:

reldiff = abs(v1-v2)/max(v1,v2)

Source code:

;-------------------------------------------------------------------
; Script: stat_relative_difference.ncl
;
; Description:  Read data variable tas from model 1 data and
;-------------  compare it with the model 2 data
;
;    Compute the relative difference:
;
;        reldiff    = (v1-v2)/max(v1,v2)
;
;        reldiffabs = abs(v1-v2)/max(v1,v2)
;
;    NCL has not a max(v1,v2) function but you can use
;        where(v1 .gt. v2, v1, v2)
;
; Output file:  plot_stat_rel_diff.png
;-------------
;
; Input files:  tas_model1.nc, tas_model2.nc
;-------------
;
; 09.12.16 kmf
;-------------------------------------------------------------------
begin
  startt = get_cpu_time()                 ;-- returns the CPU time used by NCL

  ;-------------------------------------------------------------------
  ;-- choose variable
  ;-------------------------------------------------------------------
  varname = "tas"                             ;-- tas(lat,lon)

  ;-------------------------------------------------------------------
  ;-- old data file
  ;-------------------------------------------------------------------
  diri = "$NCL_TUT/data/"                     ;-- data directory
  file_1 = "tas_model1.nc"                    ;-- file name 1
  file_2 = "tas_model2.nc"                    ;-- file name 2

  ;-------------------------------------------------------------------
  ;-- open data and grid files
  ;-------------------------------------------------------------------
  f1 = addfile(diri+file_1, "r")              ;-- open old data file 1
  f2 = addfile(diri+file_2, "r")              ;-- open new data file 2

  ;-------------------------------------------------------------------
  ;-- read the variables and lat/lon (same for both)
  ;-------------------------------------------------------------------
  tas_1 = f1->$varname$                       ;-- read variable 1
  tas_2 = f2->$varname$                       ;-- read variable 2

  ;-------------------------------------------------------------------
  ;-- check for missing value (if exists missing value/_FillValue of
  ;-- tas_1 and tas_2 must be the same
  ;-------------------------------------------------------------------
  missval = default_fillvalue(typeof(tas_1))  ;-- define missing value of type tas_1

  if(isatt(tas_1,"missing_value")) then
     if(.not. isatt(tas_1,"_FillValue")) then
        tas_1@_FillValue = tas_1@missing_value
     end if
  else
     if(.not. isatt(tas_1,"_FillValue")) then
        tas_1@_FillValue = missval
     end if
  end if

  if(isatt(tas_2,"missing_value")) then
     if(.not. isatt(tas_2,"_FillValue")) then
        tas_2@_FillValue = tas_2@missing_value
     end if
  else
     if(.not. isatt(tas_1,"_FillValue")) then
        tas_2@_FillValue = missval
     end if
  end if

  ;-------------------------------------------------------------------
  ;-- read latitudes and longitudes (must be the same for tas_1 and tas_2)
  ;-------------------------------------------------------------------
  lon = f1->lon ;-- read lon
  lon@units = "degrees_east"                  ;-- lon units

  lat = f1->lat                               ;-- read lat
  lat@units = "degrees_north"                 ;-- lat units

  ;-------------------------------------------------------------------
  ;-- compute the difference (v1-v2)
  ;-------------------------------------------------------------------
  difference = tas_1                          ;-- copy variable and retain metadata
  difference = tas_1-tas_2                    ;-- compute the difference
  difference@long_name = tas_1@long_name+" - difference" ;-- change long_name attribute

  ;-------------------------------------------------------------------
  ;-- absolute difference abs(v1-v2)
  ;-------------------------------------------------------------------
  absdiff = tas_1                             ;-- copy variable and retain metadata
  absdiff = abs(tas_1-tas_2)                  ;-- compute the absolute difference
  absdiff@long_name = tas_2@long_name+" - abs. difference" ;-- change long_name attribute

  ;-------------------------------------------------------------------
  ;-- element-wise maximum of two arrays
  ;-------------------------------------------------------------------
  emax = where(tas_1 .gt. tas_2, tas_1, tas_2) ;-- max(v1,v2)
  emax@long_name = tas_2@long_name+" - max(v1,v2)"    ;-- change long_name attribute

  ;-------------------------------------------------------------------
  ;-- relative difference 1: reldiff1 = (v1-v2)/emax
  ;-------------------------------------------------------------------
  reldiff = difference/emax                   ;-- compute the relative difference
  copy_VarMeta(tas_1,reldiff)                 ;-- copy metadata from tas_1 to reldiff
  reldiff@long_name = tas_2@long_name+" - rel. difference" ;-- change long_name attribute
  reldiff@units = ""                          ;-- change units attribute
  reldiff&lat@units = "degrees_north"         ;-- set latitude units of variable
  reldiff&lon@units = "degrees_east"          ;-- set longitude units of variable

  ;-------------------------------------------------------------------
  ;-- relative difference 2: reldiff2 = abs(v1-v2)/emax
  ;-------------------------------------------------------------------
  reldiffabs = absdiff/emax ;-- compute the relative difference
  copy_VarMeta(tas_1,reldiffabs)              ;-- copy metadata from tas_1 to reldiff2
  reldiffabs@long_name = tas_2@long_name+" - rel. difference using absolute difference"
                                              ;-- change long_name attribute
  reldiffabs@units = ""                       ;-- change units attribute
  reldiffabs&lat@units = "degrees_north"      ;-- set latitude units of variable
  reldiffabs&lon@units = "degrees_east"       ;-- set longitude units of variable

  ;-------------------------------------------------------------------
  ;-- define the region
  ;-------------------------------------------------------------------
  lon_min = min(lon) ;-- minimum longitude
  lon_max = max(lon) ;-- maximum longitude
  lat_min = min(lat) ;-- minimum latitude
  lat_max = max(lat) ;-- maximum latitude

  ;-------------------------------------------------------------------
  ;-- print some information
  ;-------------------------------------------------------------------
  print("")
  printMinMax(tas_1,0)
  printMinMax(tas_2,0)
  print("")
  printMinMax(difference,0)
  printMinMax(absdiff,0)
  printMinMax(emax,0)
  printMinMax(reldiff,0)
  printMinMax(reldiffabs,0)
  print("")
  print("MinLon: " + min(lon) + " / MaxLon: " + max(lon)+" / MinLat: "+min(lat)+" / MaxLat: "+max(lat))
  print("")

  ;-------------------------------------------------------------------
  ;-- open workstation, define type, size and background color
  ;-------------------------------------------------------------------
  wks_type = "png" ;-- plot output type
  wks_type@wkWidth = 800 ;-- workstation width
  wks_type@wkHeight = 800 ;-- workstation height
  wks_type@wkBackgroundColor = "grey90"
  wks = gsn_open_wks(wks_type,"plot_stat_rel_diff_panel_diff_abs")

  ;-------------------------------------------------------------------
  ;-- set common resources
  ;-------------------------------------------------------------------
  res = True
  res@gsnDraw = False                         ;-- don't draw plot yet
  res@gsnFrame = False                        ;-- don't advance frame
  res@gsnLeftString = ""                      ;-- don't draw left string
  res@gsnRightString = ""                     ;-- don't draw right string
  res@gsnAddCyclic = False                    ;-- regional data
  res@gsnMaximize = True                      ;-- maximize output

  ;-- contour resources
  res@cnFillOn = True                         ;-- contour fill on
  res@cnFillPalette = "rainbow"               ;-- choose colormap
  res@cnLinesOn = False                       ;-- don't draw contour lines
  res@cnInfoLabelOn = False                   ;-- switch off contour info label
  res@cnFillMode = "RasterFill"               ;-- set fill mode
  res@cnLevelSelectionMode = "ManualLevels"   ;-- use manual settings

  ;-- map resources
  res@mpFillOn = False                        ;-- turn off fill map grey
  res@mpDataBaseVersion = "MediumRes"         ;-- use medium map resolution
  res@mpDataSetName = "Earth..4"              ;-- choose map data set
  res@mpOutlineOn = True                      ;-- draw coastlines
  res@mpLimitMode = "LatLon"                  ;-- choose map limit mode
  res@mpMinLatF = lat_min                     ;-- sub-region minimum latitude
  res@mpMaxLatF = lat_max                     ;-- sub-region maximum latitude
  res@mpMinLonF = lon_min                     ;-- sub-region minimum longitude
  res@mpMaxLonF = lon_max                     ;-- sub-region maximum longitude
  res@mpGridAndLimbOn = True                  ;-- draw grid lines
  res@mpGridLonSpacingF = 5.                  ;-- lon grid line spacing
  res@mpGridLatSpacingF = 5.                  ;-- lat grid line spacing
  res@mpGridLineColor = "black"               ;-- grid line color
  res@pmTickMarkDisplayMode = "Always"        ;-- turn on tickmarks

  ;-- axes resources
  res@tmYLLabelFontHeightF = 0.010            ;-- left y-axis font size
  res@tmXBLabelFontHeightF = 0.010            ;-- bottom x-axis font size

  ;-- labelbar resources
  res@lbOrientation = "vertical"              ;-- vertical labelbar
  res@lbLabelFontHeightF = 0.012              ;-- labelbar labe font size
  res@lbBoxMinorExtentF = 0.12                ;-- decrease the height of the labelbar
  res@cnLabelBarEndStyle = "ExcludeOuterBoxes" ;-- don't draw outer labelbar boxes

  ;-- title resources
  res@tiMainFontHeightF = 0.018                ;-- title font size
  res@tiMainOffsetYF = -0.005                 ;-- move title downward

  ;***********************************
  ; prepare plots
  ;***********************************

  ;-------------------------------------------------------------------
  ;-- define value minimum, maximum and interval
  ;-------------------------------------------------------------------
  dmin = min(reldiff)                         ;-- minimum value
  dmax = max(reldiff)                         ;-- maximum value
  dint = (max(reldiff)-min(reldiff))/20       ;-- value interval

  ;-------------------------------------------------------------------
  ;-- use difference: see how many diff values exist and in
  ;-- which ranges they are and print them
  ;-------------------------------------------------------------------
  print("")
  print("--> use difference to compute the relative difference")
  print("")

  reldiff1d = ndtooned(reldiff) ;-- convert 2d variable to 1d

  val = 0.
  do i=dmin,dmax,dint
     ii = i + dint
     if((i .eq. dmin) .or. (i .eq. dmax)) then
        val := ind(reldiff1d .eq. i) ;-- retrieve number of indices
        print("--> values equal "+sprintf("%5.2f",i)+": "+sprinti("%8i",dimsizes(val)))
     else
        val := ind((reldiff1d .ge. i) .and. (reldiff1d .lt. ii)) ;-- retrieve number of indices where values are in between
        print("--> values <= "+sprintf("%5.2f",i)+" > "+sprintf("%5.2f",ii)+": "+sprinti("%8i",dimsizes(val)))
     end if
  end do
  print("-----------------------------------------------")
  print("")

  ;-- variable specific resources
  title = "Relative difference: (v1-v2)/max(v1,v2)" ;-- title string

  res@cnMinLevelValF = dmin                   ;-- minimum contour value
  res@cnMaxLevelValF = dmax                   ;-- maximum contour value
  res@cnLevelSpacingF = dint                  ;-- contour interval

  res@tiMainString = title                    ;-- title string

  ;-------------------------------------------------------------------
  ;-- create first plot (using differences to compute the relative difference)
  ;-------------------------------------------------------------------
  plot1 = gsn_csm_contour_map(wks, reldiff, res) ;-- create plot

  ;-- change labelbar label format for all plots
  getvalues plot1@contour
     "cnLevels" : levels1                     ;-- retrieve contour level values
  end getvalues
  labels = sprintf("%9.6f",levels1)           ;-- labelbar label format
  res@lbLabelStrings := labels                ;-- format the labels
  plot1 = gsn_csm_contour_map(wks, reldiff, res) ;-- recreate plot

  ;***********************************
  ; prepare 'abs' plots
  ;***********************************

  ;-------------------------------------------------------------------
  ;-- define value minimum, maximum and interval
  ;-------------------------------------------------------------------
  dmin = min(reldiffabs) ;-- minimum value
  dmax = max(reldiffabs) ;-- maximum value
  dint = (max(reldiffabs)-min(reldiffabs))/20 ;-- value interval

  ;-------------------------------------------------------------------
  ;-- use absolute difference: see how many diff values exist and in
  ;-- which ranges they are and print them
  ;-------------------------------------------------------------------
  print("")
  print("--> use difference to compute the relative difference")
  print("")

  reldiffabs1d := ndtooned(reldiffabs)        ;-- convert 2d variable to 1d

  val := 0.
  do i=dmin,dmax,dint
     ii := i + dint
     if((i .eq. dmin) .or. (i .eq. dmax)) then
        val := ind(reldiffabs1d .eq. i)       ;-- retrieve number of indices
        print("--> values equal "+sprintf("%5.2f",i)+": "+sprinti("%8i",dimsizes(val)))
     else
        val := ind((reldiffabs1d .ge. i) .and. (reldiffabs1d .lt. ii)) ;-- retrieve number of indices where values are in between
        print("--> values <= "+sprintf("%5.2f",i)+" > "+sprintf("%5.2f",ii)+": "+sprinti("%8i",dimsizes(val)))
     end if
  end do
  print("")

  ;-- variable specific resources
  title = "Relative difference: abs(v1-v2)/max(v1,v2)"

  res@cnMinLevelValF = dmin                   ;-- minimum contour value
  res@cnMaxLevelValF = dmax                   ;-- maximum contour value
  res@cnLevelSpacingF = dint                  ;-- contour interval
  res@tiMainString = title                    ;-- change title

  ;-------------------------------------------------------------------
  ;-- create second plot
  ;-------------------------------------------------------------------
  plot2 = gsn_csm_contour_map(wks, reldiffabs, res) ;-- create plot

  ;-- change labelbar label format
  getvalues plot2@contour
      "cnLevels" : levels2                    ;-- retrieve contour level values
  end getvalues
  labels = sprintf("%9.6f",levels2)           ;-- labelbar label format
  res@lbLabelStrings := labels                ;-- format the labels
  plot2 = gsn_csm_contour_map(wks, reldiffabs, res) ;-- create plot

  ;-------------------------------------------------------------------
  ;-- create the panel
  ;-------------------------------------------------------------------
  pres = True
  gsn_panel(wks,(/plot1,plot2/),(/1,2/), pres) ;-- create panel plot

  ;-------------------------------------------------------------------
  ;-- print elapsed CPU time
  ;-------------------------------------------------------------------
  endt = get_cpu_time() ;-- returns the CPU time used by NCL
  print("--> Used CPU time: "+ (endt-startt) + "s") ;-- how long did it take

end

Result:

../../../../../../_images/plot_stat_rel_diff_panel_diff_abs_w400.png