Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;************************************************************
2 ; sort by latitude in decending order (N->S)
4 ; add year_ob_i and year_ob_f
6 ;************************************************************
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
10 ;************************************************************
11 procedure set_line(lines:string,nline:integer,newlines:string)
13 ; add line to ascci/html file
15 nnewlines = dimsizes(newlines)
16 if(nline+nnewlines-1.ge.dimsizes(lines))
17 print("set_line: bad index, not setting anything.")
20 lines(nline:nline+nnewlines-1) = newlines
21 ; print ("lines = " + lines(nline:nline+nnewlines-1))
22 nline = nline + nnewlines
25 ;*************************************************************
31 ; for 6 fields, 12-monthly
35 ;************************************************
37 ;************************************************
41 model_name = "i01.10" + model
43 dirm = "/fis/cgd/cseg/people/jeff/surface_data/"
45 fm = addfile(dirm+film,"r")
49 ;------------------------------------------------
53 data_mod0 = new ((/nfield,nmon,nlat,nlon/),float)
55 ; change to unit of observed (u mol/m2/s)
56 ; Model_units [=] gC/m2/s
57 ; 12. = molecular weight of C
63 ;************************************************
65 ;************************************************
67 station = (/"ARM_Oklahoma" \
68 ,"ARM_Oklahoma_burn" \
69 ,"ARM_Oklahoma_control" \
77 ,"Duke_Forest_Hardwoods" \
78 ,"Duke_Forest_Open_Field" \
82 ,"Flagstaff_Managed" \
83 ,"Flagstaff_Unmanaged" \
84 ,"Flagstaff_Wildfire" \
86 ,"FreemanRanch_mesquite" \
89 ,"HarvardForestHemlock" \
90 ,"HowlandForestMain" \
91 ,"HowlandForestWest" \
93 ,"KendallGrasslands" \
94 ,"KennedySpaceCenterPine" \
95 ,"KennedySpaceCenterScrub" \
99 ,"Mead-irrigated-rotation" \
101 ,"Metolius_2nd_YoungPonderosaPine" \
103 ,"MetoliusIntermediatePine" \
104 ,"MetoliusOldPonderosaPine" \
109 ,"NorthCarolina_cc" \
110 ,"NorthCarolina_lp" \
115 ,"SkyOaks_PostFire" \
117 ,"SylvaniaWilderness" \
143 year_ob = (/"2003-2006" \
371 field = (/"NEE Flux" \
379 field_unit = (/"u mol/m2/s" \
387 ;========================================================================
388 ; find # year observed >=4 and year_ob_i <= 2001
390 i_long_ob = ind((year_ob_f - year_ob_i) .ge. 4 .and. year_ob_i .le. 2001)
392 station_long = station(i_long_ob)
393 year_ob_long = year_ob(i_long_ob)
394 year_ob_i_long = year_ob_i(i_long_ob)
395 year_ob_f_long = year_ob_f(i_long_ob)
398 ;print (year_ob_i(i_long_ob))
399 ;print (station_long)
400 ;========================================================================
402 nstation = dimsizes(station_long)
404 data_mod = new ((/nstation, nfield, nmon/),float)
405 data_ob = new ((/nstation, nfield, nmon/),float)
406 lat_ob = new ((/nstation/),float)
407 lon_ob = new ((/nstation/),float)
408 unit_ob = new ((/nfield/),string)
410 diro_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
411 dirm_root = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
415 ;-------------------------------------------------
418 diro = diro_root + station_long(n)+"/"
419 filo = year_ob_i_long(n)+"-"+year_ob_f_long(n)+"_L4_m.nc"
420 fo = addfile (diro+filo,"r")
427 data = fo->NEE_or_fMDS
428 data_ob(n,0,:) = dim_avg(data(month|:,year|:))
432 data_ob(n,1,:) = dim_avg(data(month|:,year|:))
436 data_ob(n,2,:) = dim_avg(data(month|:,year|:))
440 data_ob(n,3,:) = dim_avg(data(month|:,year|:))
443 data = fo->GPP_or_MDS
444 data_ob(n,4,:) = dim_avg(data(month|:,year|:))
448 data_ob(n,5,:) = dim_avg(data(month|:,year|:))
452 ;---------------------------------------------------
455 ; film = model_name+"_"+ year_ob(n)+"_MONS_climo.nc"
456 film = model_name+"_"+"1990-2004"+"_MONS_climo.nc"
457 fm = addfile (dirm_root+film,"r")
461 if (ENERGY .eq. "old") then
464 data_mod0(0,:,:,:) = data(:,:,:) * factor
471 data_mod0(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
476 ; data = fm->SENSIBLE
478 data_mod0(3,:,:,:) = data(:,:,:)
484 data_mod0(1,:,:,:) = data1(:,:,:)-data2(:,:,:)-data_mod0(2,:,:,:)-data_mod0(3,:,:,:)
491 data_mod0(0,:,:,:) = data(:,:,:) * factor
495 data_mod0(1,:,:,:) = data(:,:,:)
499 data_mod0(2,:,:,:) = data(:,:,:)
502 ; data = fm->SENSIBLE
504 data_mod0(3,:,:,:) = data(:,:,:)
509 data_mod0(4,:,:,:) = data(:,:,:) * factor
512 if (model .eq. "cn") then
523 data_mod0(5,:,:,:) = data(:,:,:) * factor
528 ;************************************************************
529 ; interpolate model data into observed station
530 ; note: model is 0-360E, 90S-90N
531 ;************************************************************
533 ; to be able to handle observation at (-89.98,-24.80)
536 yy = linint2_points_Wrap(xm,ym,data_mod0,True,lon_ob(n),lat_ob(n),0)
540 data_mod(n,:,:) = yy(:,:,0)
546 ;************************************************************
547 ; compute correlation coef and M score
548 ;************************************************************
552 ccr = new ((/nstation, nfield/),float)
553 M_score = new ((/nstation, nfield/),float)
557 ccr(n,m) = esccr(data_ob(n,m,:),data_mod(n,m,:),0)
558 bias = sum(abs(data_mod(n,m,:)-data_ob(n,m,:))/(abs(data_mod(n,m,:))+abs(data_ob(n,m,:))))
559 M_score(n,m) = (1. -(bias/nmon)) * score_max
563 M_co2 = avg(M_score(:,0))
564 M_rad = avg(M_score(:,1))
565 M_lh = avg(M_score(:,2))
566 M_sh = avg(M_score(:,3))
567 M_gpp = avg(M_score(:,4))
568 M_er = avg(M_score(:,5))
569 M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er
571 M_energy_co2 = sprintf("%.2f", M_co2)
572 M_energy_rad = sprintf("%.2f", M_rad)
573 M_energy_lh = sprintf("%.2f", M_lh )
574 M_energy_sh = sprintf("%.2f", M_sh )
575 M_energy_gpp = sprintf("%.2f", M_gpp)
576 M_energy_er = sprintf("%.2f", M_er )
577 M_energy_all = sprintf("%.2f", M_all)
579 ;*******************************************************************
580 ; for station line plot
581 ;*******************************************************************
583 ; for x-axis in xyplot
585 mon@long_name = "month"
587 res = True ; plot mods desired
588 res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
589 res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
590 ;-------------------------------------------------------------------------
591 ; Add a boxed legend using the more simple method
593 res@pmLegendDisplayMode = "Always"
594 ; res@pmLegendWidthF = 0.1
595 res@pmLegendWidthF = 0.08
596 res@pmLegendHeightF = 0.06
597 ; res@pmLegendOrthogonalPosF = -1.17
598 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
599 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
601 ; res@pmLegendParallelPosF = 0.18
602 res@pmLegendParallelPosF = 0.23 ;(rightward)
604 ; res@lgPerimOn = False
605 res@lgLabelFontHeightF = 0.015
606 res@xyExplicitLegendLabels = (/"observed",model_name/)
607 ;-------------------------------------------------------------------
609 res@gsnFrame = False ; Do not draw plot
610 res@gsnDraw = False ; Do not advance frame
612 pres = True ; panel plot mods desired
613 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
614 ; indiv. plots in panel
615 pres@gsnMaximize = True ; fill the page
616 ;-------------------------------------------------------------------
618 plot_data = new((/2,12/),float)
620 plot_data!1 = "month"
624 ;----------------------------
627 plot_name = station_long(n)+"_ob"
628 title = station_long(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
629 res@tiMainString = title
631 wks = gsn_open_wks (plot_type,plot_name)
632 plot=new(nfield,graphic) ; create graphic array
635 plot_data(0,:) = (/data_ob (n,i,:)/)
636 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
637 plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res) ; create plot
640 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
642 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
643 "rm "+plot_name+"."+plot_type)
647 ;----------------------------
650 plot_name = station_long(n)+"_model_vs_ob"
651 title = station_long(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
652 res@tiMainString = title
654 wks = gsn_open_wks (plot_type,plot_name)
655 plot=new(nfield,graphic) ; create graphic array
658 plot_data(0,:) = (/data_ob (n,i,:)/)
659 plot_data(1,:) = (/data_mod(n,i,:)/)
660 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
661 plot(i)=gsn_csm_xy(wks,mon,plot_data,res) ; create plot
664 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
666 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
667 "rm "+plot_name+"."+plot_type)
673 ;###################################################################
674 ; for the following tables,
675 ; sort by latitude in decending order (N->S)
677 isort = dim_pqsort(lat_ob,-1)
679 station_sort = station_long(isort)
680 year_ob_sort = year_ob_long(isort)
681 lat_ob_sort = lat_ob(isort)
682 lon_ob_sort = lon_ob(isort)
683 M_score_sort = M_score(isort,:)
688 ;*******************************************************************
689 ; html table of site: observed
690 ;*******************************************************************
691 output_html = "line_ob.html"
693 header = (/"<HTML>" \
695 ,"<TITLE>CLAMP metrics</TITLE>" \
697 ,"<H1>Energy at Site: Observation</H1>" \
702 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
704 ," <th bgcolor=DDDDDD >Site Name</th>" \
705 ," <th bgcolor=DDDDDD >Latitude</th>" \
706 ," <th bgcolor=DDDDDD >Longitude</th>" \
707 ," <th bgcolor=DDDDDD >Observed</th>" \
710 table_footer = "</table>"
714 lines = new(50000,string)
717 set_line(lines,nline,header)
718 set_line(lines,nline,table_header)
719 ;-----------------------------------------------
724 set_line(lines,nline,row_header)
726 txt0 = station_sort(n)
727 txt1 = sprintf("%5.2f", lat_ob_sort(n))
728 txt2 = sprintf("%5.2f", lon_ob_sort(n))
729 txt3 = year_ob_sort(n)
731 set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")
732 set_line(lines,nline,"<th>"+txt1+"</th>")
733 set_line(lines,nline,"<th>"+txt2+"</th>")
734 set_line(lines,nline,"<th>"+txt3+"</th>")
736 set_line(lines,nline,row_footer)
738 ;-----------------------------------------------
739 set_line(lines,nline,table_footer)
740 set_line(lines,nline,footer)
742 ; Now write to an HTML file.
743 idx = ind(.not.ismissing(lines))
744 if(.not.any(ismissing(idx))) then
745 asciiwrite(output_html,lines(idx))
751 ;*******************************************************************
752 ; score and line table : model vs observed
753 ;*******************************************************************
754 output_html = "score+line_vs_ob.html"
756 header = (/"<HTML>" \
758 ,"<TITLE>CLAMP metrics</TITLE>" \
760 ,"<H1>Energy at Site: Model "+model_name+"</H1>" \
764 delete (table_header)
766 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
768 ," <th bgcolor=DDDDDD >Site Name</th>" \
769 ," <th bgcolor=DDDDDD >Latitude</th>" \
770 ," <th bgcolor=DDDDDD >Longitude</th>" \
771 ," <th bgcolor=DDDDDD >Observed</th>" \
772 ," <th bgcolor=DDDDDD >CO2 Flux</th>" \
773 ," <th bgcolor=DDDDDD >Net Radiation</th>" \
774 ," <th bgcolor=DDDDDD >Latent Heat</th>" \
775 ," <th bgcolor=DDDDDD >Sensible Heat</th>" \
776 ," <th bgcolor=DDDDDD >GPP Glux</th>" \
777 ," <th bgcolor=DDDDDD >Respiration</th>" \
778 ," <th bgcolor=DDDDDD >Average</th>" \
781 table_footer = "</table>"
785 lines = new(50000,string)
788 set_line(lines,nline,header)
789 set_line(lines,nline,table_header)
790 ;-----------------------------------------------
795 set_line(lines,nline,row_header)
797 txt0 = station_sort(n)
798 txt1 = sprintf("%5.2f", lat_ob_sort(n))
799 txt2 = sprintf("%5.2f", lon_ob_sort(n))
800 txt3 = year_ob_sort(n)
801 txt4 = sprintf("%5.2f", M_score_sort(n,0))
802 txt5 = sprintf("%5.2f", M_score_sort(n,1))
803 txt6 = sprintf("%5.2f", M_score_sort(n,2))
804 txt7 = sprintf("%5.2f", M_score_sort(n,3))
805 txt8 = sprintf("%5.2f", M_score_sort(n,4))
806 txt9 = sprintf("%5.2f", M_score_sort(n,5))
807 txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))
809 set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")
810 set_line(lines,nline,"<th>"+txt1+"</th>")
811 set_line(lines,nline,"<th>"+txt2+"</th>")
812 set_line(lines,nline,"<th>"+txt3+"</th>")
813 set_line(lines,nline,"<th>"+txt4+"</th>")
814 set_line(lines,nline,"<th>"+txt5+"</th>")
815 set_line(lines,nline,"<th>"+txt6+"</th>")
816 set_line(lines,nline,"<th>"+txt7+"</th>")
817 set_line(lines,nline,"<th>"+txt8+"</th>")
818 set_line(lines,nline,"<th>"+txt9+"</th>")
819 set_line(lines,nline,"<th>"+txt10+"</th>")
821 set_line(lines,nline,row_footer)
825 set_line(lines,nline,row_header)
827 txt0 = "All_"+sprintf("%.0f", nstation)
839 set_line(lines,nline,"<th>"+txt0+"</th>")
840 set_line(lines,nline,"<th>"+txt1+"</th>")
841 set_line(lines,nline,"<th>"+txt2+"</th>")
842 set_line(lines,nline,"<th>"+txt3+"</th>")
843 set_line(lines,nline,"<th>"+txt4+"</th>")
844 set_line(lines,nline,"<th>"+txt5+"</th>")
845 set_line(lines,nline,"<th>"+txt6+"</th>")
846 set_line(lines,nline,"<th>"+txt7+"</th>")
847 set_line(lines,nline,"<th>"+txt8+"</th>")
848 set_line(lines,nline,"<th>"+txt9+"</th>")
849 set_line(lines,nline,"<th>"+txt10+"</th>")
851 set_line(lines,nline,row_footer)
852 ;-----------------------------------------------
853 set_line(lines,nline,table_footer)
854 set_line(lines,nline,footer)
856 ; Now write to an HTML file.
857 idx = ind(.not.ismissing(lines))
858 if(.not.any(ismissing(idx))) then
859 asciiwrite(output_html,lines(idx))
865 ;***************************************************************************
867 ;***************************************************************************
868 output_dir = model_name+"/ameriflux"
870 system("mv *.png *.html " + output_dir)
871 ;***************************************************************************