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:
