Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ; ***********************************************
2 ; using gsn_table for all
3 ; combine 19.metric_plot.ncl and 24.lines.ncl
4 ; ***********************************************
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
8 ;************************************************
9 procedure set_line(lines:string,nline:integer,newlines:string)
11 ; add line to ascci/html file
13 nnewlines = dimsizes(newlines)
14 if(nline+nnewlines-1.ge.dimsizes(lines))
15 print("set_line: bad index, not setting anything.")
18 lines(nline:nline+nnewlines-1) = newlines
19 ; print ("lines = " + lines(nline:nline+nnewlines-1))
20 nline = nline + nnewlines
23 ;****************************************************************************
30 ;************************************************
32 ;************************************************
34 ; from command line inputs
36 ;--------------------------------------------------
37 ; edit table.html of current model for movel1_vs_model2
39 if (isvar("compare")) then
40 html_name2 = compare+"/table.html"
41 html_new2 = html_name2 +".new"
43 ; the following only needs to execute once
44 ; system("sed 1,/model_nameA/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \
45 ; "mv -f "+html_new2+" "+html_name2+";"+ \
46 ; "sed 1,/model_nameB/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \
47 ; "mv -f "+html_new2+" "+html_name2+";"+ \
48 ; "sed s#"+modeln+"#"+model_name+"# "+html_name2+" > "+html_new2+";"+ \
49 ; "mv -f "+html_new2+" "+html_name2)
52 ;-------------------------------------
53 ; edit table.html for current model
55 html_name = model_name+"/table.html"
56 html_new = html_name +".new"
58 ; the following only needs to execute once
59 ; system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
60 ; "mv -f "+html_new+" "+html_name)
61 ;------------------------------------------------
62 fm = addfile(dirm+film3,"r")
72 ; get co2 at the lowest level
75 ; change to unit of observed (u mol/mol)
76 ; Model_units [=] kgCO2 / kgDryAir
77 ; 28.966 = molecular weight of dry air
78 ; 44. = molecular weight of CO2
81 factor = (28.966/44.) * 1e6
86 ;************************************************
88 ;************************************************
89 diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/ob/"
90 fili = "co2_globalView_98.nc"
91 g = addfile (diri+fili,"r")
95 sta = chartostring(g->STATION)
99 ;**************************************************************
100 ; get only the lowest level at each station
101 ;**************************************************************
103 lat_tmp@_FillValue = 1.e+36
106 if (.not. ismissing(lat_tmp(n))) then
107 indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
108 if (dimsizes(indexes) .gt. 1) then
109 lat_tmp(indexes(1:)) = lat_tmp@_FillValue
115 indexes = ind(.not. ismissing(lat_tmp))
117 lat_ob = lat(indexes)
118 lon_ob = lon(indexes)
119 val_ob = val(indexes,:)
120 ;************************************************************
121 ; interpolate model data into observed station
122 ; note: model is 0-360E, 90S-90N
123 ;************************************************************
124 ; to be able to handle observation at (-89.98,-24.80)
127 i = ind(lon_ob .lt. 0.)
128 lon_ob(i) = lon_ob(i) + 360.
130 yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
132 val_model = yo(pts|:,time|:)
133 val_model_0 = val_model
134 ;************************************************************
136 ;************************************************************
137 val_model = val_model - conform(val_model,dim_avg(val_model),0)
139 ;*******************************************************************
140 ; res for station line plot
141 ;*******************************************************************
142 ; for x-axis in xyplot
144 mon@long_name = "month"
146 res = True ; plot mods desired
147 res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker
148 res@xyLineColors = (/"red","black"/) ; change line color
150 ; Add a boxed legend using the more simple method
151 res@pmLegendDisplayMode = "Always"
152 ; res@pmLegendWidthF = 0.1
153 res@pmLegendWidthF = 0.08
154 res@pmLegendHeightF = 0.06
155 ; res@pmLegendOrthogonalPosF = -1.17
156 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
157 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
159 ; res@pmLegendParallelPosF = 0.18
160 res@pmLegendParallelPosF = 0.23 ;(rightward)
162 ; res@lgPerimOn = False
163 res@lgLabelFontHeightF = 0.015
164 res@xyExplicitLegendLabels = (/model_name,"observed"/)
165 ;************************************************************
166 ; number of latitude zone
167 ;************************************************************
170 ; saving data for zone
171 ; number of rows for zone table (with data)
174 ; number of columns for zone table
177 text4 = new((/nrow_zone,ncol_zone/),string)
184 ind_z = ind(lat_ob .ge. 60.)
190 ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
196 ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
202 ind_z = ind(lat_ob .lt. 0. )
205 npts = dimsizes(ind_z)
207 ;------------------------------------------------------
208 ; for metric table computation
209 amp_ob = new((/npts/),float)
210 amp_model = new((/npts/),float)
212 amp_ratio_sta = new((/npts/),float)
213 ccr_sta = new((/npts/),float)
214 M_sta = new((/npts/),float)
215 score_sta = new((/npts/),float)
216 ;-----------------------------------------------------
217 ; for station line plot
221 plot_data = new((/2,12,npts/),float)
222 plot_data_0 = new((/12,npts/),float)
225 plot_data!1 = "month"
227 plot_data@long_name = "CO2 Seasonal"
229 plot_data_0!0 = "month"
230 plot_data_0!1 = "pts"
231 plot_data_0@long_name = "CO2"
232 ;--------------------------------------------------------------------------
235 amp_ob(n) = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))
236 amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
238 amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
239 ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
240 M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
241 score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
243 ;----------------------------------------------------------------------
244 ; for station line plot
246 plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
247 plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
249 plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
251 plot_name = sta(ind_z(n))
252 title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"
254 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
255 ;------------------------------------------
258 plot=new(2,graphic) ; create graphic array
259 res@gsnFrame = False ; Do not draw plot
260 res@gsnDraw = False ; Do not advance frame
262 pres = True ; panel plot mods desired
263 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
264 ; indiv. plots in panel
265 pres@gsnMaximize = True ; fill the page
266 ;------------------------------------------
267 res@tiMainString = title ; add title
269 plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res) ; create plot 1
271 plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
273 gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
275 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
276 system("rm "+plot_name+"."+plot_type)
279 ;---------------------------------------------------------------------------
281 ;-------------------------------------------------------------------------
282 ; for line plot in a zone
284 plot_name = "All_"+npts_str
285 title = plot_name + " in "+ zone
287 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
288 ;-----------------------------------------
291 plot=new(2,graphic) ; create graphic array
292 res@gsnFrame = False ; Do not draw plot
293 res@gsnDraw = False ; Do not advance frame
295 pres = True ; panel plot mods desired
296 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
297 ; indiv. plots in panel
298 pres@gsnMaximize = True ; fill the page
299 ;-----------------------------------------
300 res@tiMainString = title ; add title
302 plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res) ; create plot 1
304 plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
306 gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
308 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
309 system("rm "+plot_name+"."+plot_type)
315 ;---------------------------------------------------------------------------
316 ; values saved for zone table
318 amp_ratio_zone = avg(amp_ratio_sta)
319 ccr_zone = avg(ccr_sta)
320 M_zone = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)
321 score_zone = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
324 text4(z,1) = sprintf("%.0f", npts)
325 text4(z,2) = sprintf("%.2f", amp_ratio_zone)
326 text4(z,3) = sprintf("%.2f", ccr_zone)
327 text4(z,4) = sprintf("%.2f", M_zone)
328 text4(z,5) = sprintf("%.2f", score_zone)
330 ;*******************************************************************
331 ; html table -- station
332 ;*******************************************************************
333 output_html = "score+line_"+zone+".html"
335 header = (/"<HTML>" \
337 ,"<TITLE>CLAMP metrics</TITLE>" \
339 ,"<H1>CO2 in Zone "+zone+": Model "+model_name+"</H1>" \
344 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
346 ," <th bgcolor=DDDDDD >Site Name</th>" \
347 ," <th bgcolor=DDDDDD >Latitude</th>" \
348 ," <th bgcolor=DDDDDD >Longitude</th>" \
349 ," <th bgcolor=DDDDDD >Amplitude Ratio</th>" \
350 ," <th bgcolor=DDDDDD >Correleration Coef</th>" \
351 ," <th bgcolor=DDDDDD >M Score</th>" \
352 ," <th bgcolor=DDDDDD >Combined Score</th>" \
355 table_footer = "</table>"
359 lines = new(50000,string)
362 set_line(lines,nline,header)
363 set_line(lines,nline,table_header)
364 ;-----------------------------------------------
368 set_line(lines,nline,row_header)
371 txt1 = sprintf("%5.2f", (/lat(ind_z(n))/))
372 txt2 = sprintf("%5.2f", (/lon(ind_z(n))/))
373 txt3 = sprintf("%5.2f", (/amp_ratio_sta(n)/))
374 txt4 = sprintf("%5.2f", (/ccr_sta(n)/))
375 txt5 = sprintf("%5.2f", (/M_sta(n)/))
376 txt6 = sprintf("%5.2f", (/score_sta(n)/))
378 set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")
379 set_line(lines,nline,"<th>"+txt1+"</th>")
380 set_line(lines,nline,"<th>"+txt2+"</th>")
381 set_line(lines,nline,"<th>"+txt3+"</th>")
382 set_line(lines,nline,"<th>"+txt4+"</th>")
383 set_line(lines,nline,"<th>"+txt5+"</th>")
384 set_line(lines,nline,"<th>"+txt6+"</th>")
386 set_line(lines,nline,row_footer)
390 set_line(lines,nline,row_header)
392 txt0 = "All_"+sprintf("%.0f", (/npts/))
395 txt3 = sprintf("%5.2f", (/amp_ratio_zone/))
396 txt4 = sprintf("%5.2f", (/ccr_zone/))
397 txt5 = sprintf("%5.2f", (/M_zone/))
398 txt6 = sprintf("%5.2f", (/score_zone/))
400 set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")
401 set_line(lines,nline,"<th>"+txt1+"</th>")
402 set_line(lines,nline,"<th>"+txt2+"</th>")
403 set_line(lines,nline,"<th>"+txt3+"</th>")
404 set_line(lines,nline,"<th>"+txt4+"</th>")
405 set_line(lines,nline,"<th>"+txt5+"</th>")
406 set_line(lines,nline,"<th>"+txt6+"</th>")
408 set_line(lines,nline,row_footer)
409 ;-----------------------------------------------
410 set_line(lines,nline,table_footer)
411 set_line(lines,nline,footer)
413 ; Now write to an HTML file.
414 idx = ind(.not.ismissing(lines))
415 if(.not.any(ismissing(idx))) then
416 asciiwrite(output_html,lines(idx))
421 ;-----------------------------------------------------------------
426 delete (amp_ratio_sta)
433 ;*******************************************************************
435 ;*******************************************************************
436 output_html = "score+line_vs_ob.html"
438 header = (/"<HTML>" \
440 ,"<TITLE>CLAMP metrics</TITLE>" \
442 ,"<H1>CO2 in Latitude Zone: Model "+model_name+"</H1>" \
446 delete (table_header)
448 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
450 ," <th bgcolor=DDDDDD >Zone</th>" \
451 ," <th bgcolor=DDDDDD >Number of Site</th>" \
452 ," <th bgcolor=DDDDDD >Amplitide Ratio</th>" \
453 ," <th bgcolor=DDDDDD >Correleration Coef</th>" \
454 ," <th bgcolor=DDDDDD >M Score</th>" \
455 ," <th bgcolor=DDDDDD >Combined Score</th>" \
458 table_footer = "</table>"
462 lines = new(50000,string)
465 set_line(lines,nline,header)
466 set_line(lines,nline,table_header)
467 ;-----------------------------------------------
471 set_line(lines,nline,row_header)
473 set_line(lines,nline,"<th><a href=score+line_"+text4(n,0)+".html>"+text4(n,0)+"</th>")
474 set_line(lines,nline,"<th>"+text4(n,1)+"</th>")
475 set_line(lines,nline,"<th>"+text4(n,2)+"</th>")
476 set_line(lines,nline,"<th>"+text4(n,3)+"</th>")
477 set_line(lines,nline,"<th>"+text4(n,4)+"</th>")
478 set_line(lines,nline,"<th>"+text4(n,5)+"</th>")
480 set_line(lines,nline,row_footer)
486 txt1 = sum(stringtofloat(text4(0:3,1)))
490 txt5 = sum(stringtofloat(text4(0:3,5)))
492 set_line(lines,nline,row_header)
494 set_line(lines,nline,"<th>"+txt0+"</th>")
495 set_line(lines,nline,"<th>"+txt1+"</th>")
496 set_line(lines,nline,"<th>"+txt2+"</th>")
497 set_line(lines,nline,"<th>"+txt3+"</th>")
498 set_line(lines,nline,"<th>"+txt4+"</th>")
499 set_line(lines,nline,"<th>"+txt5+"</th>")
501 set_line(lines,nline,row_footer)
502 ;-----------------------------------------------
503 set_line(lines,nline,table_footer)
504 set_line(lines,nline,footer)
506 ; Now write to an HTML file.
507 idx = ind(.not.ismissing(lines))
508 if(.not.any(ismissing(idx))) then
509 asciiwrite(output_html,lines(idx))
514 ;--------------------------------------------------------------------------
518 if (isvar("compare")) then
519 system("sed 1,/M_co2/s//"+M_co2+"/ "+html_name2+" > "+html_new2+";"+ \
520 "mv -f "+html_new2+" "+html_name2)
523 system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \
524 "mv -f "+html_new+" "+html_name)
525 ;***************************************************************************
527 ;***************************************************************************
528 output_dir = model_name+"/co2"
530 system("mv *.png *.html " + output_dir)
531 ;***************************************************************************