forrest@0: ; ***********************************************
forrest@0: ; using gsn_table for all
forrest@0: ; combine 19.metric_plot.ncl and 24.lines.ncl
forrest@0: ; ***********************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;************************************************
forrest@0: procedure set_line(lines:string,nline:integer,newlines:string)
forrest@0: begin
forrest@0: ; add line to ascci/html file
forrest@0:
forrest@0: nnewlines = dimsizes(newlines)
forrest@0: if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0: print("set_line: bad index, not setting anything.")
forrest@0: return
forrest@0: end if
forrest@0: lines(nline:nline+nnewlines-1) = newlines
forrest@0: ; print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0: nline = nline + nnewlines
forrest@0: return
forrest@0: end
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: begin
forrest@0:
forrest@0: plot_type = "ps"
forrest@0: plot_type_new = "png"
forrest@0:
forrest@0: ;************************************************
forrest@0: ; read model data
forrest@0: ;************************************************
forrest@0:
forrest@0: ; from command line inputs
forrest@0:
forrest@0: ;--------------------------------------------------
forrest@0: ; edit table.html of current model for movel1_vs_model2
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: html_name2 = compare+"/table.html"
forrest@0: html_new2 = html_name2 +".new"
forrest@0:
forrest@0: ; the following only needs to execute once
forrest@0: ; system("sed 1,/model_nameA/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0: ; "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: ; "sed 1,/model_nameB/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0: ; "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: ; "sed s#"+modeln+"#"+model_name+"# "+html_name2+" > "+html_new2+";"+ \
forrest@0: ; "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: ;-------------------------------------
forrest@0: ; edit table.html for current model
forrest@0:
forrest@0: html_name = model_name+"/table.html"
forrest@0: html_new = html_name +".new"
forrest@0:
forrest@0: ; the following only needs to execute once
forrest@0: ; system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
forrest@0: ; "mv -f "+html_new+" "+html_name)
forrest@0: ;------------------------------------------------
forrest@0: fm = addfile(dirm+film3,"r")
forrest@0:
forrest@0: x = fm->CO2
forrest@0: xi = fm->lon
forrest@0: yi = fm->lat
forrest@0:
forrest@0: xdim = dimsizes(x)
forrest@0: nlev = xdim(1)
forrest@0: y = x(:,0,:,:)
forrest@0:
forrest@0: ; get co2 at the lowest level
forrest@0: y = x(:,nlev-1,:,:)
forrest@0:
forrest@0: ; change to unit of observed (u mol/mol)
forrest@0: ; Model_units [=] kgCO2 / kgDryAir
forrest@0: ; 28.966 = molecular weight of dry air
forrest@0: ; 44. = molecular weight of CO2
forrest@0: ; u mol = 1e-6 mol
forrest@0:
forrest@0: factor = (28.966/44.) * 1e6
forrest@0: y = y * factor
forrest@0:
forrest@0: y@_FillValue = 1.e36
forrest@0: y@units = "u mol/mol"
forrest@0: ;************************************************
forrest@0: ; read data: observed
forrest@0: ;************************************************
forrest@0: diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/ob/"
forrest@0: fili = "co2_globalView_98.nc"
forrest@0: g = addfile (diri+fili,"r")
forrest@0: val = g->CO2_SEAS
forrest@0: lon = g->LON
forrest@0: lat = g->LAT
forrest@0: sta = chartostring(g->STATION)
forrest@0: delete (g)
forrest@0:
forrest@0: ncase = dimsizes(lat)
forrest@0: ;**************************************************************
forrest@0: ; get only the lowest level at each station
forrest@0: ;**************************************************************
forrest@0: lat_tmp = lat
forrest@0: lat_tmp@_FillValue = 1.e+36
forrest@0:
forrest@0: do n = 0,ncase-1
forrest@0: if (.not. ismissing(lat_tmp(n))) then
forrest@0: indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
forrest@0: if (dimsizes(indexes) .gt. 1) then
forrest@0: lat_tmp(indexes(1:)) = lat_tmp@_FillValue
forrest@0: end if
forrest@0: delete (indexes)
forrest@0: end if
forrest@0: end do
forrest@0:
forrest@0: indexes = ind(.not. ismissing(lat_tmp))
forrest@0:
forrest@0: lat_ob = lat(indexes)
forrest@0: lon_ob = lon(indexes)
forrest@0: val_ob = val(indexes,:)
forrest@0: ;************************************************************
forrest@0: ; interpolate model data into observed station
forrest@0: ; note: model is 0-360E, 90S-90N
forrest@0: ;************************************************************
forrest@0: ; to be able to handle observation at (-89.98,-24.80)
forrest@0: yi(0) = -90.
forrest@0:
forrest@0: i = ind(lon_ob .lt. 0.)
forrest@0: lon_ob(i) = lon_ob(i) + 360.
forrest@0:
forrest@0: yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
forrest@0:
forrest@0: val_model = yo(pts|:,time|:)
forrest@0: val_model_0 = val_model
forrest@0: ;************************************************************
forrest@0: ; remove annual mean
forrest@0: ;************************************************************
forrest@0: val_model = val_model - conform(val_model,dim_avg(val_model),0)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; res for station line plot
forrest@0: ;*******************************************************************
forrest@0: ; for x-axis in xyplot
forrest@0: mon = ispan(1,12,1)
forrest@0: mon@long_name = "month"
forrest@0:
forrest@0: res = True ; plot mods desired
forrest@0: res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker
forrest@0: res@xyLineColors = (/"red","black"/) ; change line color
forrest@0:
forrest@0: ; Add a boxed legend using the more simple method
forrest@0: res@pmLegendDisplayMode = "Always"
forrest@0: ; res@pmLegendWidthF = 0.1
forrest@0: res@pmLegendWidthF = 0.08
forrest@0: res@pmLegendHeightF = 0.06
forrest@0: ; res@pmLegendOrthogonalPosF = -1.17
forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0: res@pmLegendOrthogonalPosF = -0.30 ;(downward)
forrest@0:
forrest@0: ; res@pmLegendParallelPosF = 0.18
forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward)
forrest@0:
forrest@0: ; res@lgPerimOn = False
forrest@0: res@lgLabelFontHeightF = 0.015
forrest@0: res@xyExplicitLegendLabels = (/model_name,"observed"/)
forrest@0: ;************************************************************
forrest@0: ; number of latitude zone
forrest@0: ;************************************************************
forrest@0: nzone = 4
forrest@0:
forrest@0: ; saving data for zone
forrest@0: ; number of rows for zone table (with data)
forrest@0: nrow_zone = nzone
forrest@0:
forrest@0: ; number of columns for zone table
forrest@0: ncol_zone = 6
forrest@0:
forrest@0: text4 = new((/nrow_zone,ncol_zone/),string)
forrest@0:
forrest@0: do z = 0,nzone-1
forrest@0:
forrest@0: if (z .eq. 0) then
forrest@0: zone = "60N-90N"
forrest@0: score_max = 5.0
forrest@0: ind_z = ind(lat_ob .ge. 60.)
forrest@0: end if
forrest@0:
forrest@0: if (z .eq. 1) then
forrest@0: zone = "30N-60N"
forrest@0: score_max = 5.0
forrest@0: ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
forrest@0: end if
forrest@0:
forrest@0: if (z .eq. 2) then
forrest@0: zone = "EQ-30N"
forrest@0: score_max = 5.0
forrest@0: ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
forrest@0: end if
forrest@0:
forrest@0: if (z .eq. 3) then
forrest@0: zone = "90S-EQ"
forrest@0: score_max = 5.0
forrest@0: ind_z = ind(lat_ob .lt. 0. )
forrest@0: end if
forrest@0:
forrest@0: npts = dimsizes(ind_z)
forrest@0:
forrest@0: ;------------------------------------------------------
forrest@0: ; for metric table computation
forrest@0: amp_ob = new((/npts/),float)
forrest@0: amp_model = new((/npts/),float)
forrest@0:
forrest@0: amp_ratio_sta = new((/npts/),float)
forrest@0: ccr_sta = new((/npts/),float)
forrest@0: M_sta = new((/npts/),float)
forrest@0: score_sta = new((/npts/),float)
forrest@0: ;-----------------------------------------------------
forrest@0: ; for station line plot
forrest@0: npts_str = ""
forrest@0: npts_str = npts
forrest@0:
forrest@0: plot_data = new((/2,12,npts/),float)
forrest@0: plot_data_0 = new((/12,npts/),float)
forrest@0:
forrest@0: plot_data!0 = "case"
forrest@0: plot_data!1 = "month"
forrest@0: plot_data!2 = "pts"
forrest@0: plot_data@long_name = "CO2 Seasonal"
forrest@0:
forrest@0: plot_data_0!0 = "month"
forrest@0: plot_data_0!1 = "pts"
forrest@0: plot_data_0@long_name = "CO2"
forrest@0: ;--------------------------------------------------------------------------
forrest@0: do n=0,npts-1
forrest@0:
forrest@0: amp_ob(n) = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))
forrest@0: amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
forrest@0:
forrest@0: amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
forrest@0: ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
forrest@0: M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
forrest@0: score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
forrest@0:
forrest@0: ;----------------------------------------------------------------------
forrest@0: ; for station line plot
forrest@0:
forrest@0: plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
forrest@0: plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
forrest@0:
forrest@0: plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
forrest@0:
forrest@0: plot_name = sta(ind_z(n))
forrest@0: title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: ;------------------------------------------
forrest@0: ; for panel plot
forrest@0:
forrest@0: plot=new(2,graphic) ; create graphic array
forrest@0: res@gsnFrame = False ; Do not draw plot
forrest@0: res@gsnDraw = False ; Do not advance frame
forrest@0:
forrest@0: pres = True ; panel plot mods desired
forrest@0: pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0: ; indiv. plots in panel
forrest@0: pres@gsnMaximize = True ; fill the page
forrest@0: ;------------------------------------------
forrest@0: res@tiMainString = title ; add title
forrest@0:
forrest@0: plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res) ; create plot 1
forrest@0:
forrest@0: plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
forrest@0:
forrest@0: gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
forrest@0:
forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0: system("rm "+plot_name+"."+plot_type)
forrest@0:
forrest@0: clear (wks)
forrest@0: ;---------------------------------------------------------------------------
forrest@0: end do
forrest@0: ;-------------------------------------------------------------------------
forrest@0: ; for line plot in a zone
forrest@0:
forrest@0: plot_name = "All_"+npts_str
forrest@0: title = plot_name + " in "+ zone
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: ;-----------------------------------------
forrest@0: ; for panel plot
forrest@0:
forrest@0: plot=new(2,graphic) ; create graphic array
forrest@0: res@gsnFrame = False ; Do not draw plot
forrest@0: res@gsnDraw = False ; Do not advance frame
forrest@0:
forrest@0: pres = True ; panel plot mods desired
forrest@0: pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0: ; indiv. plots in panel
forrest@0: pres@gsnMaximize = True ; fill the page
forrest@0: ;-----------------------------------------
forrest@0: res@tiMainString = title ; add title
forrest@0:
forrest@0: plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res) ; create plot 1
forrest@0:
forrest@0: plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
forrest@0:
forrest@0: gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
forrest@0:
forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0: system("rm "+plot_name+"."+plot_type)
forrest@0:
forrest@0: clear (wks)
forrest@0: ; delete (ind_z)
forrest@0: delete (plot_data)
forrest@0: delete (plot_data_0)
forrest@0: ;---------------------------------------------------------------------------
forrest@0: ; values saved for zone table
forrest@0:
forrest@0: amp_ratio_zone = avg(amp_ratio_sta)
forrest@0: ccr_zone = avg(ccr_sta)
forrest@0: M_zone = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)
forrest@0: score_zone = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
forrest@0:
forrest@0: text4(z,0) = zone
forrest@0: text4(z,1) = sprintf("%.0f", npts)
forrest@0: text4(z,2) = sprintf("%.2f", amp_ratio_zone)
forrest@0: text4(z,3) = sprintf("%.2f", ccr_zone)
forrest@0: text4(z,4) = sprintf("%.2f", M_zone)
forrest@0: text4(z,5) = sprintf("%.2f", score_zone)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; html table -- station
forrest@0: ;*******************************************************************
forrest@0: output_html = "score+line_"+zone+".html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"
" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"CO2 in Zone "+zone+": Model "+model_name+"
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," Site Name | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@0: ," Amplitude Ratio | " \
forrest@0: ," Correleration Coef | " \
forrest@0: ," M Score | " \
forrest@0: ," Combined Score | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ; row of table
forrest@0:
forrest@0: do n = 0,npts-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt0 = sta(ind_z(n))
forrest@0: txt1 = sprintf("%5.2f", (/lat(ind_z(n))/))
forrest@0: txt2 = sprintf("%5.2f", (/lon(ind_z(n))/))
forrest@0: txt3 = sprintf("%5.2f", (/amp_ratio_sta(n)/))
forrest@0: txt4 = sprintf("%5.2f", (/ccr_sta(n)/))
forrest@0: txt5 = sprintf("%5.2f", (/M_sta(n)/))
forrest@0: txt6 = sprintf("%5.2f", (/score_sta(n)/))
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0:
forrest@0: ; last row, summary
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt0 = "All_"+sprintf("%.0f", (/npts/))
forrest@0: txt1 = "-"
forrest@0: txt2 = "-"
forrest@0: txt3 = sprintf("%5.2f", (/amp_ratio_zone/))
forrest@0: txt4 = sprintf("%5.2f", (/ccr_zone/))
forrest@0: txt5 = sprintf("%5.2f", (/M_zone/))
forrest@0: txt6 = sprintf("%5.2f", (/score_zone/))
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0: delete (idx)
forrest@0: ;-----------------------------------------------------------------
forrest@0:
forrest@0: delete (ind_z)
forrest@0: delete (amp_model)
forrest@0: delete (amp_ob)
forrest@0: delete (amp_ratio_sta)
forrest@0: delete (ccr_sta)
forrest@0: delete (M_sta)
forrest@0: delete (score_sta)
forrest@0: clear (wks)
forrest@0: end do
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; html table -- zone
forrest@0: ;*******************************************************************
forrest@0: output_html = "score+line_vs_ob.html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"CO2 in Latitude Zone: Model "+model_name+"
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: delete (table_header)
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," Zone | " \
forrest@0: ," Number of Site | " \
forrest@0: ," Amplitide Ratio | " \
forrest@0: ," Correleration Coef | " \
forrest@0: ," M Score | " \
forrest@0: ," Combined Score | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: do n = 0,nrow_zone-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: set_line(lines,nline,""+text4(n,0)+" | ")
forrest@0: set_line(lines,nline,""+text4(n,1)+" | ")
forrest@0: set_line(lines,nline,""+text4(n,2)+" | ")
forrest@0: set_line(lines,nline,""+text4(n,3)+" | ")
forrest@0: set_line(lines,nline,""+text4(n,4)+" | ")
forrest@0: set_line(lines,nline,""+text4(n,5)+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0:
forrest@0: ; for the last row
forrest@0:
forrest@0: txt0 = "All"
forrest@0: txt1 = sum(stringtofloat(text4(0:3,1)))
forrest@0: txt2 = "-"
forrest@0: txt3 = "-"
forrest@0: txt4 = "-"
forrest@0: txt5 = sum(stringtofloat(text4(0:3,5)))
forrest@0:
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: set_line(lines,nline,""+txt0+" | ")
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0: delete (idx)
forrest@0: ;--------------------------------------------------------------------------
forrest@0:
forrest@0: M_co2 = txt5
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed 1,/M_co2/s//"+M_co2+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0: ;***************************************************************************
forrest@0: ; output plots
forrest@0: ;***************************************************************************
forrest@0: output_dir = model_name+"/co2"
forrest@0:
forrest@0: system("mv *.png *.html " + output_dir)
forrest@0: ;***************************************************************************
forrest@0: end