Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;************************************************************
2 ; required command line input parameters:
3 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
5 ; using gsn_table for all
6 ; output: line plot for each site (4 fields)
8 ;************************************************************
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
12 ;************************************************************
13 procedure set_line(lines:string,nline:integer,newlines:string)
15 ; add line to ascci/html file
17 nnewlines = dimsizes(newlines)
18 if(nline+nnewlines-1.ge.dimsizes(lines))
19 print("set_line: bad index, not setting anything.")
22 lines(nline:nline+nnewlines-1) = newlines
23 ; print ("lines = " + lines(nline:nline+nnewlines-1))
24 nline = nline + nnewlines
27 ;*************************************************************
33 ; for 6 fields, 12-monthly
37 ;************************************************
39 ;************************************************
43 model_name = "i01.10" + model
45 dirm = "/fis/cgd/cseg/people/jeff/surface_data/"
47 fm = addfile(dirm+film,"r")
51 ;------------------------------------------------
55 data_mod0 = new ((/nfield,nmon,nlat,nlon/),float)
57 ; change to unit of observed (u mol/m2/s)
58 ; Model_units [=] gC/m2/s
59 ; 12. = molecular weight of C
65 ;************************************************
67 ;************************************************
69 station = (/"ARM_Oklahoma" \
70 ,"ARM_Oklahoma_burn" \
71 ,"ARM_Oklahoma_control" \
79 ,"Duke_Forest_Hardwoods" \
80 ,"Duke_Forest_Open_Field" \
84 ,"Flagstaff_Managed" \
85 ,"Flagstaff_Unmanaged" \
86 ,"Flagstaff_Wildfire" \
88 ,"FreemanRanch_mesquite" \
91 ,"HarvardForestHemlock" \
92 ,"HowlandForestMain" \
93 ,"HowlandForestWest" \
95 ,"KendallGrasslands" \
96 ,"KennedySpaceCenterPine" \
97 ,"KennedySpaceCenterScrub" \
101 ,"Mead-irrigated-rotation" \
103 ,"Metolius_2nd_YoungPonderosaPine" \
105 ,"MetoliusIntermediatePine" \
106 ,"MetoliusOldPonderosaPine" \
111 ,"NorthCarolina_cc" \
112 ,"NorthCarolina_lp" \
117 ,"SkyOaks_PostFire" \
119 ,"SylvaniaWilderness" \
145 year_ob = (/"2003-2006" \
221 field = (/"CO2 Flux" \
229 nstation = dimsizes(station)
231 data_mod = new ((/nstation, nfield, nmon/),float)
232 data_ob = new ((/nstation, nfield, nmon/),float)
233 lat_ob = new ((/nstation/),float)
234 lon_ob = new ((/nstation/),float)
236 diro_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
237 dirm_root = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
241 ;-------------------------------------------------
244 diro = diro_root + station(n)+"/"
245 filo = year_ob(n)+"_L4_m.nc"
246 fo = addfile (diro+filo,"r")
253 data = fo->NEE_or_fMDS
254 data_ob(n,0,:) = dim_avg(data(month|:,year|:))
258 data_ob(n,1,:) = dim_avg(data(month|:,year|:))
262 data_ob(n,2,:) = dim_avg(data(month|:,year|:))
266 data_ob(n,3,:) = dim_avg(data(month|:,year|:))
269 data = fo->GPP_or_MDS
270 data_ob(n,4,:) = dim_avg(data(month|:,year|:))
274 data_ob(n,5,:) = dim_avg(data(month|:,year|:))
278 ;---------------------------------------------------
281 ; film = model_name+"_"+ year_ob(n)+"_MONS_climo.nc"
282 film = model_name+"_"+"1990-2004"+"_MONS_climo.nc"
283 fm = addfile (dirm_root+film,"r")
287 if (ENERGY .eq. "old") then
290 data_mod0(0,:,:,:) = data(:,:,:) * factor
297 data_mod0(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
302 ; data = fm->SENSIBLE
304 data_mod0(3,:,:,:) = data(:,:,:)
310 data_mod0(1,:,:,:) = data1(:,:,:)-data2(:,:,:)-data_mod0(2,:,:,:)-data_mod0(3,:,:,:)
317 data_mod0(0,:,:,:) = data(:,:,:) * factor
321 data_mod0(1,:,:,:) = data(:,:,:)
325 data_mod0(2,:,:,:) = data(:,:,:)
328 ; data = fm->SENSIBLE
330 data_mod0(3,:,:,:) = data(:,:,:)
335 data_mod0(4,:,:,:) = data(:,:,:) * factor
338 if (model .eq. "cn") then
349 data_mod0(5,:,:,:) = data(:,:,:) * factor
354 ;************************************************************
355 ; interpolate model data into observed station
356 ; note: model is 0-360E, 90S-90N
357 ;************************************************************
359 ; to be able to handle observation at (-89.98,-24.80)
362 yy = linint2_points_Wrap(xm,ym,data_mod0,True,lon_ob(n),lat_ob(n),0)
366 data_mod(n,:,:) = yy(:,:,0)
371 ;************************************************************
372 ; compute correlation coef and M score
373 ;************************************************************
377 ccr = new ((/nstation, nfield/),float)
378 M_score = new ((/nstation, nfield/),float)
382 ccr(n,m) = esccr(data_ob(n,m,:),data_mod(n,m,:),0)
383 bias = sum(abs(data_mod(n,m,:)-data_ob(n,m,:))/(abs(data_mod(n,m,:))+abs(data_ob(n,m,:))))
384 M_score(n,m) = (1. -(bias/nmon)) * score_max
388 M_co2 = avg(M_score(:,0))
389 M_rad = avg(M_score(:,1))
390 M_lh = avg(M_score(:,2))
391 M_sh = avg(M_score(:,3))
392 M_gpp = avg(M_score(:,4))
393 M_er = avg(M_score(:,5))
394 M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er
396 M_energy_co2 = sprintf("%.2f", M_co2)
397 M_energy_rad = sprintf("%.2f", M_rad)
398 M_energy_lh = sprintf("%.2f", M_lh )
399 M_energy_sh = sprintf("%.2f", M_sh )
400 M_energy_gpp = sprintf("%.2f", M_gpp)
401 M_energy_er = sprintf("%.2f", M_er )
402 M_energy_all = sprintf("%.2f", M_all)
404 ;*******************************************************************
405 ; for station line plot
406 ;*******************************************************************
408 ; for x-axis in xyplot
410 mon@long_name = "month"
412 res = True ; plot mods desired
413 res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
414 res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
415 ;-------------------------------------------------------------------------
416 ; Add a boxed legend using the more simple method
418 res@pmLegendDisplayMode = "Always"
419 ; res@pmLegendWidthF = 0.1
420 res@pmLegendWidthF = 0.08
421 res@pmLegendHeightF = 0.06
422 ; res@pmLegendOrthogonalPosF = -1.17
423 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
424 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
426 ; res@pmLegendParallelPosF = 0.18
427 res@pmLegendParallelPosF = 0.23 ;(rightward)
429 ; res@lgPerimOn = False
430 res@lgLabelFontHeightF = 0.015
431 res@xyExplicitLegendLabels = (/"observed",model_name/)
432 ;-------------------------------------------------------------------
434 res@gsnFrame = False ; Do not draw plot
435 res@gsnDraw = False ; Do not advance frame
437 pres = True ; panel plot mods desired
438 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
439 ; indiv. plots in panel
440 pres@gsnMaximize = True ; fill the page
441 ;-------------------------------------------------------------------
443 plot_data = new((/2,12/),float)
445 plot_data!1 = "month"
448 ;----------------------------
451 plot_name = station(n)+"_ob"
452 title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
453 res@tiMainString = title
455 wks = gsn_open_wks (plot_type,plot_name)
456 plot=new(nfield,graphic) ; create graphic array
458 plot_data(0,:) = (/data_ob (n,0,:)/)
459 plot_data@long_name = field(0)
460 plot(0)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot 1
462 plot_data(0,:) = (/data_ob (n,1,:)/)
463 plot_data@long_name = field(1)
464 plot(1)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot 2
466 plot_data(0,:) = (/data_ob (n,2,:)/)
467 plot_data@long_name = field(2)
468 plot(2)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot 3
470 plot_data(0,:) = (/data_ob (n,3,:)/)
471 plot_data@long_name = field(3)
472 plot(3)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot 4
474 plot_data(0,:) = (/data_ob (n,4,:)/)
475 plot_data@long_name = field(4)
476 plot(4)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot 5
478 plot_data(0,:) = (/data_ob (n,5,:)/)
479 plot_data@long_name = field(5)
480 plot(5)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot 6
482 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
484 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
485 "rm "+plot_name+"."+plot_type)
489 ;----------------------------
492 plot_name = station(n)+"_model_vs_ob"
493 title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
494 res@tiMainString = title
496 wks = gsn_open_wks (plot_type,plot_name)
497 plot=new(nfield,graphic) ; create graphic array
499 plot_data(0,:) = (/data_ob (n,0,:)/)
500 plot_data(1,:) = (/data_mod(n,0,:)/)
501 plot_data@long_name = field(0)
502 plot(0)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot 1
504 plot_data(0,:) = (/data_ob (n,1,:)/)
505 plot_data(1,:) = (/data_mod(n,1,:)/)
506 plot_data@long_name = field(1)
507 plot(1)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot 2
509 plot_data(0,:) = (/data_ob (n,2,:)/)
510 plot_data(1,:) = (/data_mod(n,2,:)/)
511 plot_data@long_name = field(2)
512 plot(2)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot 3
514 plot_data(0,:) = (/data_ob (n,3,:)/)
515 plot_data(1,:) = (/data_mod(n,3,:)/)
516 plot_data@long_name = field(3)
517 plot(3)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot 4
520 plot_data(0,:) = (/data_ob (n,4,:)/)
521 plot_data(1,:) = (/data_mod(n,4,:)/)
522 plot_data@long_name = field(4)
523 plot(4)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot 5
526 plot_data(0,:) = (/data_ob (n,5,:)/)
527 plot_data(1,:) = (/data_mod(n,5,:)/)
528 plot_data@long_name = field(5)
529 plot(5)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot 6
531 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
533 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
534 "rm "+plot_name+"."+plot_type)
540 ;*******************************************************************
541 ; html table of site: observed
542 ;*******************************************************************
543 output_html = "line_ob.html"
545 header = (/"<HTML>" \
547 ,"<TITLE>CLAMP metrics</TITLE>" \
549 ,"<H1>Energy at Site: Observation</H1>" \
554 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
556 ," <th bgcolor=DDDDDD >Site Name</th>" \
557 ," <th bgcolor=DDDDDD >Latitude</th>" \
558 ," <th bgcolor=DDDDDD >Longitude</th>" \
559 ," <th bgcolor=DDDDDD >Observed</th>" \
562 table_footer = "</table>"
566 lines = new(50000,string)
569 set_line(lines,nline,header)
570 set_line(lines,nline,table_header)
571 ;-----------------------------------------------
575 set_line(lines,nline,row_header)
578 txt1 = sprintf("%5.2f", lat_ob(n))
579 txt2 = sprintf("%5.2f", lon_ob(n))
582 set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")
583 set_line(lines,nline,"<th>"+txt1+"</th>")
584 set_line(lines,nline,"<th>"+txt2+"</th>")
585 set_line(lines,nline,"<th>"+txt3+"</th>")
587 set_line(lines,nline,row_footer)
589 ;-----------------------------------------------
590 set_line(lines,nline,table_footer)
591 set_line(lines,nline,footer)
593 ; Now write to an HTML file.
594 idx = ind(.not.ismissing(lines))
595 if(.not.any(ismissing(idx))) then
596 asciiwrite(output_html,lines(idx))
602 ;*******************************************************************
603 ; score and line table : model vs observed
604 ;*******************************************************************
605 output_html = "score+line_vs_ob.html"
607 header = (/"<HTML>" \
609 ,"<TITLE>CLAMP metrics</TITLE>" \
611 ,"<H1>Energy at Site: Model "+model_name+"</H1>" \
615 delete (table_header)
617 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
619 ," <th bgcolor=DDDDDD >Site Name</th>" \
620 ," <th bgcolor=DDDDDD >Latitude</th>" \
621 ," <th bgcolor=DDDDDD >Longitude</th>" \
622 ," <th bgcolor=DDDDDD >Observed</th>" \
623 ," <th bgcolor=DDDDDD >CO2 Flux</th>" \
624 ," <th bgcolor=DDDDDD >Net Radiation</th>" \
625 ," <th bgcolor=DDDDDD >Latent Heat</th>" \
626 ," <th bgcolor=DDDDDD >Sensible Heat</th>" \
627 ," <th bgcolor=DDDDDD >GPP Glux</th>" \
628 ," <th bgcolor=DDDDDD >Respiration</th>" \
629 ," <th bgcolor=DDDDDD >Average</th>" \
632 table_footer = "</table>"
636 lines = new(50000,string)
639 set_line(lines,nline,header)
640 set_line(lines,nline,table_header)
641 ;-----------------------------------------------
645 set_line(lines,nline,row_header)
648 txt1 = sprintf("%5.2f", lat_ob(n))
649 txt2 = sprintf("%5.2f", lon_ob(n))
651 txt4 = sprintf("%5.2f", M_score(n,0))
652 txt5 = sprintf("%5.2f", M_score(n,1))
653 txt6 = sprintf("%5.2f", M_score(n,2))
654 txt7 = sprintf("%5.2f", M_score(n,3))
655 txt8 = sprintf("%5.2f", M_score(n,4))
656 txt9 = sprintf("%5.2f", M_score(n,5))
657 txt10 = sprintf("%5.2f", avg(M_score(n,:)))
659 set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")
660 set_line(lines,nline,"<th>"+txt1+"</th>")
661 set_line(lines,nline,"<th>"+txt2+"</th>")
662 set_line(lines,nline,"<th>"+txt3+"</th>")
663 set_line(lines,nline,"<th>"+txt4+"</th>")
664 set_line(lines,nline,"<th>"+txt5+"</th>")
665 set_line(lines,nline,"<th>"+txt6+"</th>")
666 set_line(lines,nline,"<th>"+txt7+"</th>")
667 set_line(lines,nline,"<th>"+txt8+"</th>")
668 set_line(lines,nline,"<th>"+txt9+"</th>")
669 set_line(lines,nline,"<th>"+txt10+"</th>")
671 set_line(lines,nline,row_footer)
675 set_line(lines,nline,row_header)
677 txt0 = "All_"+sprintf("%.0f", nstation)
689 set_line(lines,nline,"<th>"+txt0+"</th>")
690 set_line(lines,nline,"<th>"+txt1+"</th>")
691 set_line(lines,nline,"<th>"+txt2+"</th>")
692 set_line(lines,nline,"<th>"+txt3+"</th>")
693 set_line(lines,nline,"<th>"+txt4+"</th>")
694 set_line(lines,nline,"<th>"+txt5+"</th>")
695 set_line(lines,nline,"<th>"+txt6+"</th>")
696 set_line(lines,nline,"<th>"+txt7+"</th>")
697 set_line(lines,nline,"<th>"+txt8+"</th>")
698 set_line(lines,nline,"<th>"+txt9+"</th>")
699 set_line(lines,nline,"<th>"+txt10+"</th>")
701 set_line(lines,nline,row_footer)
702 ;-----------------------------------------------
703 set_line(lines,nline,table_footer)
704 set_line(lines,nline,footer)
706 ; Now write to an HTML file.
707 idx = ind(.not.ismissing(lines))
708 if(.not.any(ismissing(idx))) then
709 asciiwrite(output_html,lines(idx))
715 ;***************************************************************************
717 ;***************************************************************************
718 output_dir = model_name+"/ameriflux"
720 system("mv *.png *.html " + output_dir)
721 ;***************************************************************************