Benutzerspezifische Werkzeuge

##### Sektionen
Sie sind hier: NCL / 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 takeend`

Artikelaktionen