Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;************************************************************
3 ; sort by latitude in decending order (N->S)
5 ; add year_ob_i and year_ob_f
7 ;************************************************************
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
11 ;************************************************************
12 procedure set_line(lines:string,nline:integer,newlines:string)
14 ; add line to ascci/html file
16 nnewlines = dimsizes(newlines)
17 if(nline+nnewlines-1.ge.dimsizes(lines))
18 print("set_line: bad index, not setting anything.")
21 lines(nline:nline+nnewlines-1) = newlines
22 ; print ("lines = " + lines(nline:nline+nnewlines-1))
23 nline = nline + nnewlines
26 ;*************************************************************
32 ;************************************************
34 ;************************************************
38 model_name = "i01.10" + model
40 ; change to unit of observed (u mol/m2/s)
41 ; Model_units [=] gC/m2/s
42 ; 12. = molecular weight of C
48 ;************************************************
50 ;************************************************
52 station = (/"ARM_Oklahoma" \
53 ,"ARM_Oklahoma_burn" \
54 ,"ARM_Oklahoma_control" \
62 ,"Duke_Forest_Hardwoods" \
63 ,"Duke_Forest_Open_Field" \
67 ,"Flagstaff_Managed" \
68 ,"Flagstaff_Unmanaged" \
69 ,"Flagstaff_Wildfire" \
71 ,"FreemanRanch_mesquite" \
74 ,"HarvardForestHemlock" \
75 ,"HowlandForestMain" \
76 ,"HowlandForestWest" \
78 ,"KendallGrasslands" \
79 ,"KennedySpaceCenterPine" \
80 ,"KennedySpaceCenterScrub" \
84 ,"Mead-irrigated-rotation" \
86 ,"Metolius_2nd_YoungPonderosaPine" \
88 ,"MetoliusIntermediatePine" \
89 ,"MetoliusOldPonderosaPine" \
100 ,"SkyOaks_PostFire" \
102 ,"SylvaniaWilderness" \
128 year_ob = (/"2003-2006" \
356 field = (/"NEE Flux" \
364 field_unit = (/"u mol/m2/s" \
372 nfield = dimsizes(field)
374 ;========================================================================
375 ; find (# of year observed) >=4 and year_ob_i <= 2001
377 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
379 nstation = dimsizes(station)
381 year_station = new ((/nstation/),integer)
385 diro = dir_root + station(n)+"/"
386 filo = year_ob_i(n)+"-"+year_ob_f(n)+"_L4_m.nc"
387 fo = addfile (diro+filo,"r")
390 year_station(n) = dimsizes(year)
397 i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
399 station_long = station(i_long_ob)
400 year_ob_long = year_ob(i_long_ob)
401 year_ob_i_long = year_ob_i(i_long_ob)
402 year_ob_f_long = year_ob_f(i_long_ob)
403 year_station_long = year_station(i_long_ob)
405 nstation = dimsizes(station_long)
408 ; print (year_ob_i(i_long_ob))
409 ; print (station_long)
411 ;========================================================================
412 ; get observed lat and lon
414 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
416 lat_ob = new ((/nstation/),float)
417 lon_ob = new ((/nstation/),float)
421 diro = dir_root + station_long(n)+"/"
422 filo = year_ob_i_long(n)+"-"+year_ob_f_long(n)+"_L4_m.nc"
423 fo = addfile (diro+filo,"r")
435 ;=========================================================
436 ; get model data at observed lat-lon
438 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
439 film = model_name+"_ameriflux_1990-2004_monthly.nc"
440 fm = addfile (dirm+film,"r")
446 date_dim = dimsizes(date)
451 data_mod = new ((/nfield,nyear,nmonth,nstation/),float)
453 ;************************************************************
454 ; interpolate model data into observed station
455 ; note: model is 0-360E, 90S-90N
456 ;************************************************************
458 ; to be able to handle observation at (-89.98,-24.80)
461 if (ENERGY .eq. "old") then
464 data_mod(0,:,:,:) = data(:,:,:) * factor
470 data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)
475 ; data = fm->SENSIBLE
477 data_mod(3,:,:,:) = data(:,:,:)
482 data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:)
489 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
490 data_mod(0,:,:,:) = yy(:,:,:) * factor
493 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
494 data_mod(1,:,:,:) = yy(:,:,:)
497 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
498 data_mod(2,:,:,:) = yy(:,:,:)
500 ; data = fm->SENSIBLE
502 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
503 data_mod(3,:,:,:) = yy(:,:,:)
508 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
509 data_mod(4,:,:,:) = yy(:,:,:) * factor
511 if (model .eq. "cn") then
522 yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
523 data_mod(5,:,:,:) = yy(:,:,:) * factor
528 ;*******************************************************************
529 ; for station line plot
530 ;*******************************************************************
532 ; for x-axis in xyplot
534 mon@long_name = "month"
536 res = True ; plot mods desired
537 res@xyLineThicknesses = (/2.0,2.0/) ; make 2nd lines thicker
538 res@xyLineColors = (/"blue","red"/) ; line color (ob,model)
540 res@tmXBFormat = "f" ; not to add trailing zeros
542 ;-------------------------------------------------------------------------
543 ; Add a boxed legend using the more simple method
545 res@pmLegendDisplayMode = "Always"
546 ; res@pmLegendWidthF = 0.1
547 res@pmLegendWidthF = 0.08
548 res@pmLegendHeightF = 0.06
549 ; res@pmLegendOrthogonalPosF = -1.17
550 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
551 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
553 ; res@pmLegendParallelPosF = 0.18
554 res@pmLegendParallelPosF = 0.23 ;(rightward)
556 ; res@lgPerimOn = False
557 res@lgLabelFontHeightF = 0.015
558 res@xyExplicitLegendLabels = (/"observed",model_name/)
559 ;-------------------------------------------------------------------
561 res@gsnFrame = False ; Do not draw plot
562 res@gsnDraw = False ; Do not advance frame
564 pres = True ; panel plot mods desired
565 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
566 ; indiv. plots in panel
567 pres@gsnMaximize = True ; fill the page
568 ;-------------------------------------------------------------------
570 ;==============================================================
571 ; get ob data at each site
573 dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"
580 nyear = year_station_long(n)
582 if (year_ob_f_long(n).eq. 2006) then
583 year_setback = 2006 -2004
585 if (year_ob_f_long(n).eq. 2005) then
586 year_setback = 2005 -2004
589 ntime = (nyear - year_setback) * nmonth
592 data_ob = new ((/nfield, nyear, nmonth/),float)
594 diro = dir_root + station_long(n)+"/"
595 filo = year_ob_i_long(n)+"-"+year_ob_f_long(n)+"_L4_m.nc"
596 fo = addfile (diro+filo,"r")
604 data = fo->NEE_or_fMDS
605 data_ob(0,:,:) = data(:,:)
609 data_ob(1,:,:) = data(:,:)
613 data_ob(2,:,:) = data(:,:)
617 data_ob(3,:,:) = data(:,:)
620 data = fo->GPP_or_MDS
621 data_ob(4,:,:) = data(:,:)
625 data_ob(5,:,:) = data(:,:)
630 timeI = new((/ntime/),integer)
631 timeF = new((/ntime/),float)
632 timeI = ispan(1,ntime,1)
633 timeF = year_ob_i_long(n) + (timeI-1)/12.
634 timeF@long_name = "year"
636 plot_data = new((/2,ntime/),float)
638 ;----------------------------
641 plot_name = station_long(n)+"_model_vs_ob"
642 title = station_long(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
643 res@tiMainString = title
645 wks = gsn_open_wks (plot_type,plot_name)
646 plot=new(nfield,graphic) ; create graphic array
648 i_year_mod_i = year_ob_i_long(n) - 1990
649 i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
651 i_year_ob_f = nyear - year_setback - 1
659 plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
660 plot_data(1,:) = ndtooned(data_mod(i,i_year_mod_i:i_year_mod_f,:,n))
661 plot_data@long_name = field(i)+" ("+field_unit(i)+")"
662 plot(i)=gsn_csm_xy(wks,timeF,plot_data,res) ; create plot
665 gsn_panel(wks,plot,(/3,2/),pres) ; create panel plot
667 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
668 "rm "+plot_name+"."+plot_type)
680 ;###################################################################
681 ; for the following tables,
682 ; sort by latitude in decending order (N->S)
684 isort = dim_pqsort(lat_ob,-1)
686 station_sort = station_long(isort)
687 year_ob_sort = year_ob_long(isort)
688 lat_ob_sort = lat_ob(isort)
689 lon_ob_sort = lon_ob(isort)
690 M_score_sort = M_score(isort,:)
695 ;*******************************************************************
696 ; html table of site: observed
697 ;*******************************************************************
698 output_html = "line_ob.html"
700 header = (/"<HTML>" \
702 ,"<TITLE>CLAMP metrics</TITLE>" \
704 ,"<H1>Energy at Site: Observation</H1>" \
709 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
711 ," <th bgcolor=DDDDDD >Site Name</th>" \
712 ," <th bgcolor=DDDDDD >Latitude</th>" \
713 ," <th bgcolor=DDDDDD >Longitude</th>" \
714 ," <th bgcolor=DDDDDD >Observed</th>" \
717 table_footer = "</table>"
721 lines = new(50000,string)
724 set_line(lines,nline,header)
725 set_line(lines,nline,table_header)
726 ;-----------------------------------------------
731 set_line(lines,nline,row_header)
733 txt0 = station_sort(n)
734 txt1 = sprintf("%5.2f", lat_ob_sort(n))
735 txt2 = sprintf("%5.2f", lon_ob_sort(n))
736 txt3 = year_ob_sort(n)
738 set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")
739 set_line(lines,nline,"<th>"+txt1+"</th>")
740 set_line(lines,nline,"<th>"+txt2+"</th>")
741 set_line(lines,nline,"<th>"+txt3+"</th>")
743 set_line(lines,nline,row_footer)
745 ;-----------------------------------------------
746 set_line(lines,nline,table_footer)
747 set_line(lines,nline,footer)
749 ; Now write to an HTML file.
750 idx = ind(.not.ismissing(lines))
751 if(.not.any(ismissing(idx))) then
752 asciiwrite(output_html,lines(idx))
758 ;*******************************************************************
759 ; score and line table : model vs observed
760 ;*******************************************************************
761 output_html = "score+line_vs_ob.html"
763 header = (/"<HTML>" \
765 ,"<TITLE>CLAMP metrics</TITLE>" \
767 ,"<H1>Energy at Site: Model "+model_name+"</H1>" \
771 delete (table_header)
773 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
775 ," <th bgcolor=DDDDDD >Site Name</th>" \
776 ," <th bgcolor=DDDDDD >Latitude</th>" \
777 ," <th bgcolor=DDDDDD >Longitude</th>" \
778 ," <th bgcolor=DDDDDD >Observed</th>" \
779 ," <th bgcolor=DDDDDD >CO2 Flux</th>" \
780 ," <th bgcolor=DDDDDD >Net Radiation</th>" \
781 ," <th bgcolor=DDDDDD >Latent Heat</th>" \
782 ," <th bgcolor=DDDDDD >Sensible Heat</th>" \
783 ," <th bgcolor=DDDDDD >GPP Glux</th>" \
784 ," <th bgcolor=DDDDDD >Respiration</th>" \
785 ," <th bgcolor=DDDDDD >Average</th>" \
788 table_footer = "</table>"
792 lines = new(50000,string)
795 set_line(lines,nline,header)
796 set_line(lines,nline,table_header)
797 ;-----------------------------------------------
802 set_line(lines,nline,row_header)
804 txt0 = station_sort(n)
805 txt1 = sprintf("%5.2f", lat_ob_sort(n))
806 txt2 = sprintf("%5.2f", lon_ob_sort(n))
807 txt3 = year_ob_sort(n)
808 txt4 = sprintf("%5.2f", M_score_sort(n,0))
809 txt5 = sprintf("%5.2f", M_score_sort(n,1))
810 txt6 = sprintf("%5.2f", M_score_sort(n,2))
811 txt7 = sprintf("%5.2f", M_score_sort(n,3))
812 txt8 = sprintf("%5.2f", M_score_sort(n,4))
813 txt9 = sprintf("%5.2f", M_score_sort(n,5))
814 txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))
816 set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")
817 set_line(lines,nline,"<th>"+txt1+"</th>")
818 set_line(lines,nline,"<th>"+txt2+"</th>")
819 set_line(lines,nline,"<th>"+txt3+"</th>")
820 set_line(lines,nline,"<th>"+txt4+"</th>")
821 set_line(lines,nline,"<th>"+txt5+"</th>")
822 set_line(lines,nline,"<th>"+txt6+"</th>")
823 set_line(lines,nline,"<th>"+txt7+"</th>")
824 set_line(lines,nline,"<th>"+txt8+"</th>")
825 set_line(lines,nline,"<th>"+txt9+"</th>")
826 set_line(lines,nline,"<th>"+txt10+"</th>")
828 set_line(lines,nline,row_footer)
832 set_line(lines,nline,row_header)
834 txt0 = "All_"+sprintf("%.0f", nstation)
846 set_line(lines,nline,"<th>"+txt0+"</th>")
847 set_line(lines,nline,"<th>"+txt1+"</th>")
848 set_line(lines,nline,"<th>"+txt2+"</th>")
849 set_line(lines,nline,"<th>"+txt3+"</th>")
850 set_line(lines,nline,"<th>"+txt4+"</th>")
851 set_line(lines,nline,"<th>"+txt5+"</th>")
852 set_line(lines,nline,"<th>"+txt6+"</th>")
853 set_line(lines,nline,"<th>"+txt7+"</th>")
854 set_line(lines,nline,"<th>"+txt8+"</th>")
855 set_line(lines,nline,"<th>"+txt9+"</th>")
856 set_line(lines,nline,"<th>"+txt10+"</th>")
858 set_line(lines,nline,row_footer)
859 ;-----------------------------------------------
860 set_line(lines,nline,table_footer)
861 set_line(lines,nline,footer)
863 ; Now write to an HTML file.
864 idx = ind(.not.ismissing(lines))
865 if(.not.any(ismissing(idx))) then
866 asciiwrite(output_html,lines(idx))
872 ;***************************************************************************
874 ;***************************************************************************
875 output_dir = model_name+"/ameriflux"
877 system("mv *.png *.html " + output_dir)
878 ;***************************************************************************