Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL example relative difference plot
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

DKRZ NCL example relative difference plot

This example shows how to compute the relative difference using on one side the difference of the variables and on the other side the absolute value of the difference.

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  meier-fleischer(at)dkrz.de
;-------------------------------------------------------------------
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")     ;-- open a workstation
 
;-------------------------------------------------------------------
;-- 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 third 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

Artikelaktionen