clamp: ameriflux/11.plot.ncl@4be95183fbcd

     1 ;************************************************************

     2 ; required command line input parameters:

     3 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

4 ;

     5 ; using gsn_table for all

     6 ; output: line plot for each site (4 fields)

     7 ;         table for M_score

     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)

    14 begin

    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.")

    20     return

    21   end if

    22   lines(nline:nline+nnewlines-1) = newlines

    23 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;*************************************************************

    28 begin

    30   plot_type     = "ps"

    31   plot_type_new = "png"

    33 ; for 6 fields, 12-monthly

    34   nmon      = 12

    35   nfield    = 6

    37 ;************************************************

    38 ; read model data

    39 ;************************************************

    40 ; model = "cn"

    41   model = "casa"

    43   model_name = "i01.10" + model

    45   dirm  = "/fis/cgd/cseg/people/jeff/surface_data/"

    46   film  = "lnd_T42.nc"

    47   fm    = addfile(dirm+film,"r")

    49   xm    = fm->lon

    50   ym    = fm->lat

    51 ;------------------------------------------------

    52   nlat = dimsizes(ym)

    53   nlon = dimsizes(xm)

    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

    60 ; u mol = 1e-6 mol

    61   factor = 1e6 /12.

    63   ENERGY ="new"

    65 ;************************************************

    66 ; read data: observed

    67 ;************************************************

    69  station = (/"Duke_Forest_Hardwoods" \

    70             ,"Duke_Forest_Open_Field" \

    71             ,"Duke_Forest_Pine" \

    72             ,"HarvardForest" \

    73             /)

    75  year_ob = (/"2003-2005" \

    76             ,"2001-2005" \

    77             ,"2001-2005" \

    78             ,"1991-2004" \

    79             /)

    81  field   = (/"CO2 Flux" \

    82             ,"Net Radiation" \

    83             ,"Latent Heat" \

    84             ,"Sensible Heat" \

    85             ,"GPP Flux" \

    86             ,"Respiration" \

    87             /)

    89  nstation  = dimsizes(station)

    91  data_mod  = new ((/nstation, nfield, nmon/),float)

    92  data_ob   = new ((/nstation, nfield, nmon/),float)

    93  lat_ob    = new ((/nstation/),float)

    94  lon_ob    = new ((/nstation/),float)

    96  diro_root  = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/"

    97  dirm_root  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    99  do n = 0,nstation-1

   101 ;-------------------------------------------------

   102 ;   get ob data

   104     diro = diro_root + station(n)+"/"

   105     filo = year_ob(n)+"_L4_m.nc"

   106     fo   = addfile (diro+filo,"r")

   108     print (filo)

   110     lon_ob(n) = fo->lon

   111     lat_ob(n) = fo->lat

   113     data      = fo->NEE_or_fMDS

   114     data_ob(n,0,:) = dim_avg(data(month|:,year|:))

   115     delete (data)

   117     data      = fo->Rg_f

   118     data_ob(n,1,:) = dim_avg(data(month|:,year|:))

   119     delete (data)

   121     data      = fo->LE_f

   122     data_ob(n,2,:) = dim_avg(data(month|:,year|:))

   123     delete (data)

   125     data      = fo->H_f

   126     data_ob(n,3,:) = dim_avg(data(month|:,year|:))

   127     delete (data)

   129     data      = fo->GPP_or_MDS

   130     data_ob(n,4,:) = dim_avg(data(month|:,year|:))

   131     delete (data)

   133     data      = fo->Reco_or

   134     data_ob(n,5,:) = dim_avg(data(month|:,year|:))

   135     delete (data)

   137     delete (fo)

   138 ;---------------------------------------------------

   139 ;   get model data

   141     film = model_name+"_"+year_ob(n)+"_MONS_climo.nc"

   142     fm   = addfile (dirm_root+film,"r")

   144     print (film)

   146 if (ENERGY .eq. "old") then

   148   data = fm->NEE

   149   data_mod0(0,:,:,:) = data(:,:,:) * factor

   150   delete (data)

   152 ; data  = fm->LATENT

   153   data1 = fm->FCEV

   154   data2 = fm->FCTR

   155   data3 = fm->FGEV

   156   data_mod0(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:)

   157   delete (data1)

   158   delete (data2)

   159   delete (data3)

   161 ; data = fm->SENSIBLE

   162   data  = fm->FSH

   163   data_mod0(3,:,:,:) = data(:,:,:)

   164   delete (data)

   166 ; data  = fm->NETRAD

   167   data1 = fm->FSA

   168   data2 = fm->FIRA

   169   data_mod0(1,:,:,:) = data1(:,:,:)-data2(:,:,:)-data_mod0(2,:,:,:)-data_mod0(3,:,:,:)

   170   delete (data1)

   171   delete (data2)

   173 else

   175   data = fm->NEE

   176   data_mod0(0,:,:,:) = data(:,:,:) * factor

   177   delete (data)

   179   data = fm->NETRAD

   180   data_mod0(1,:,:,:) = data(:,:,:)

   181   delete (data)

   183   data = fm->LATENT

   184   data_mod0(2,:,:,:) = data(:,:,:)

   185   delete (data)

   187 ; data = fm->SENSIBLE

   188   data = fm->FSH

   189   data_mod0(3,:,:,:) = data(:,:,:)

   190   delete (data)

   191 end if

   193   data = fm->GPP

   194   data_mod0(4,:,:,:) = data(:,:,:) * factor

   195   delete (data)

   197   if (model .eq. "cn") then

   198      data = fm->ER

   199   else

   200      data1 = fm->AR

   201      data2 = fm->HR

   202      data  = data1 + data2

   204      delete (data1)

   205      delete (data2)

   206   end if

   208   data_mod0(5,:,:,:) = data(:,:,:) * factor

   209   delete (data)

   211   delete (fm)

   213 ;************************************************************

   214 ; interpolate model data into observed station

   215 ; note: model is 0-360E, 90S-90N

   216 ;************************************************************

   218 ; to be able to handle observation at (-89.98,-24.80)

   219   ym(0) = -90.

   221   yy = linint2_points_Wrap(xm,ym,data_mod0,True,lon_ob(n),lat_ob(n),0)

   223 ; print (yy)

   225   data_mod(n,:,:) = yy(:,:,0)

   227   delete (yy)

   229  end do

   230 ;************************************************************

   231 ; compute correlation coef and M score

   232 ;************************************************************

   234  score_max = 5.

   236  ccr     = new ((/nstation, nfield/),float)

   237  M_score = new ((/nstation, nfield/),float)

   239  do n=0,nstation-1

   240  do m=0,nfield-1

   241     ccr(n,m) = esccr(data_ob(n,m,:),data_mod(n,m,:),0)

   242     bias = sum(abs(data_mod(n,m,:)-data_ob(n,m,:))/(abs(data_mod(n,m,:))+abs(data_ob(n,m,:))))

   243     M_score(n,m) = (1. -(bias/nmon)) * score_max

   244  end do

   245  end do

   247  M_co2 = avg(M_score(:,0))

   248  M_rad = avg(M_score(:,1))

   249  M_lh  = avg(M_score(:,2))

   250  M_sh  = avg(M_score(:,3))

   251  M_gpp = avg(M_score(:,4))

   252  M_er  = avg(M_score(:,5))

   253  M_all = M_co2+ M_rad +M_lh + M_sh + M_gpp + M_er

   255  M_energy_co2 = sprintf("%.2f", M_co2)

   256  M_energy_rad = sprintf("%.2f", M_rad)

   257  M_energy_lh  = sprintf("%.2f", M_lh )

   258  M_energy_sh  = sprintf("%.2f", M_sh )

   259  M_energy_gpp = sprintf("%.2f", M_gpp)

   260  M_energy_er  = sprintf("%.2f", M_er )

   261  M_energy_all = sprintf("%.2f", M_all)

   263 ;*******************************************************************

   264 ; for station line plot

   265 ;*******************************************************************

   267 ; for x-axis in xyplot

   268   mon = ispan(1,12,1)

   269   mon@long_name = "month"

   271   res                   = True               ; plot mods desired

   272   res@xyLineThicknesses = (/2.0,2.0/)        ; make 2nd lines thicker

   273   res@xyLineColors      = (/"blue","red"/)   ; line color (ob,model)

   274 ;-------------------------------------------------------------------------

   275 ; Add a boxed legend using the more simple method

   277   res@pmLegendDisplayMode    = "Always"

   278 ; res@pmLegendWidthF         = 0.1

   279   res@pmLegendWidthF         = 0.08

   280   res@pmLegendHeightF        = 0.06

   281 ; res@pmLegendOrthogonalPosF = -1.17

   282 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   283   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   285 ; res@pmLegendParallelPosF   =  0.18

   286   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   288 ; res@lgPerimOn             = False

   289   res@lgLabelFontHeightF     = 0.015

   290   res@xyExplicitLegendLabels = (/"observed",model_name/)

   291 ;-------------------------------------------------------------------

   292 ; for panel plot

   293   res@gsnFrame     = False                   ; Do not draw plot

   294   res@gsnDraw      = False                   ; Do not advance frame

   296   pres                            = True     ; panel plot mods desired

   297   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   298                                              ; indiv. plots in panel

   299   pres@gsnMaximize                = True     ; fill the page

   300 ;-------------------------------------------------------------------

   302   plot_data   = new((/2,12/),float)

   303   plot_data!0 = "case"

   304   plot_data!1 = "month"

   306   do n = 0,nstation-1

   307 ;----------------------------

   308 ; for observed

   310     plot_name = station(n)+"_ob"

   311     title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"

   312     res@tiMainString = title

   314     wks = gsn_open_wks (plot_type,plot_name)

   315     plot=new(nfield,graphic)                        ; create graphic array

   317     plot_data(0,:) = (/data_ob (n,0,:)/)

   318     plot_data@long_name = field(0)

   319     plot(0)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 1

   321     plot_data(0,:) = (/data_ob (n,1,:)/)

   322     plot_data@long_name = field(1)

   323     plot(1)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 2

   325     plot_data(0,:) = (/data_ob (n,2,:)/)

   326     plot_data@long_name = field(2)

   327     plot(2)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 3

   329     plot_data(0,:) = (/data_ob (n,3,:)/)

   330     plot_data@long_name = field(3)

   331     plot(3)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 4

   333     plot_data(0,:) = (/data_ob (n,4,:)/)

   334     plot_data@long_name = field(4)

   335     plot(4)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 5

   337     plot_data(0,:) = (/data_ob (n,5,:)/)

   338     plot_data@long_name = field(5)

   339     plot(5)=gsn_csm_xy(wks,mon,plot_data(0,:),res)   ; create plot 6

   341     gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot

   343     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   344            "rm "+plot_name+"."+plot_type)

   346     clear (wks)

   347     delete (plot)

   348 ;----------------------------

   349 ; for model_vs_ob

   351     plot_name = station(n)+"_model_vs_ob"

   352     title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"

   353     res@tiMainString = title

   355     wks = gsn_open_wks (plot_type,plot_name)

   356     plot=new(nfield,graphic)                        ; create graphic array

   358     plot_data(0,:) = (/data_ob (n,0,:)/)

   359     plot_data(1,:) = (/data_mod(n,0,:)/)

   360     plot_data@long_name = field(0)

   361     plot(0)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 1

   363     plot_data(0,:) = (/data_ob (n,1,:)/)

   364     plot_data(1,:) = (/data_mod(n,1,:)/)

   365     plot_data@long_name = field(1)

   366     plot(1)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 2

   368     plot_data(0,:) = (/data_ob (n,2,:)/)

   369     plot_data(1,:) = (/data_mod(n,2,:)/)

   370     plot_data@long_name = field(2)

   371     plot(2)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 3

   373     plot_data(0,:) = (/data_ob (n,3,:)/)

   374     plot_data(1,:) = (/data_mod(n,3,:)/)

   375     plot_data@long_name = field(3)

   376     plot(3)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 4

   379     plot_data(0,:) = (/data_ob (n,4,:)/)

   380     plot_data(1,:) = (/data_mod(n,4,:)/)

   381     plot_data@long_name = field(4)

   382     plot(4)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 5

   385     plot_data(0,:) = (/data_ob (n,5,:)/)

   386     plot_data(1,:) = (/data_mod(n,5,:)/)

   387     plot_data@long_name = field(5)

   388     plot(5)=gsn_csm_xy(wks,mon,plot_data,res)   ; create plot 6

   390     gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot

   392     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   393            "rm "+plot_name+"."+plot_type)

   395     clear (wks)

   396     delete (plot)

   397  end do

   399 ;*******************************************************************

   400 ; html table of site: observed

   401 ;*******************************************************************

   402   output_html = "line_ob.html"

   404   header = (/"<HTML>" \

   405             ,"<HEAD>" \

   406             ,"<TITLE>CLAMP metrics</TITLE>" \

   407             ,"</HEAD>" \

   408             ,"<H1>Energy at Site: Observation</H1>" \

   409             /)

   410   footer = "</HTML>"

   412   table_header = (/ \

   413         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \

   414        ,"<tr>" \

   415        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   416        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   417        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   418        ,"   <th bgcolor=DDDDDD >Observed</th>" \

   419        ,"</tr>" \

   420        /)

   421   table_footer = "</table>"

   422   row_header = "<tr>"

   423   row_footer = "</tr>"

   425   lines = new(50000,string)

   426   nline = 0

   428   set_line(lines,nline,header)

   429   set_line(lines,nline,table_header)

   430 ;-----------------------------------------------

   431 ; row of table

   433   do n = 0,nstation-1

   434      set_line(lines,nline,row_header)

   436      txt0 = station(n)

   437      txt1 = sprintf("%5.2f", lat_ob(n))

   438      txt2 = sprintf("%5.2f", lon_ob(n))

   439      txt3 = year_ob(n)

   441      set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")

   442      set_line(lines,nline,"<th>"+txt1+"</th>")

   443      set_line(lines,nline,"<th>"+txt2+"</th>")

   444      set_line(lines,nline,"<th>"+txt3+"</th>")

   446      set_line(lines,nline,row_footer)

   447   end do

   448 ;-----------------------------------------------

   449   set_line(lines,nline,table_footer)

   450   set_line(lines,nline,footer)

   452 ; Now write to an HTML file.

   453   idx = ind(.not.ismissing(lines))

   454   if(.not.any(ismissing(idx))) then

   455     asciiwrite(output_html,lines(idx))

   456   else

   457    print ("error?")

   458   end if

   459   delete (idx)

   461 ;*******************************************************************

   462 ; score and line table : model vs observed

   463 ;*******************************************************************

   464   output_html = "score+line_vs_ob.html"

   466   header = (/"<HTML>" \

   467             ,"<HEAD>" \

   468             ,"<TITLE>CLAMP metrics</TITLE>" \

   469             ,"</HEAD>" \

   470             ,"<H1>Energy at Site: Model "+model_name+"</H1>" \

   471             /)

   472   footer = "</HTML>"

   474   delete (table_header)

   475   table_header = (/ \

   476         "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \

   477        ,"<tr>" \

   478        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   479        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   480        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   481        ,"   <th bgcolor=DDDDDD >Observed</th>" \

   482        ,"   <th bgcolor=DDDDDD >CO2 Flux</th>" \

   483        ,"   <th bgcolor=DDDDDD >Net Radiation</th>" \

   484        ,"   <th bgcolor=DDDDDD >Latent Heat</th>" \

   485        ,"   <th bgcolor=DDDDDD >Sensible Heat</th>" \

   486        ,"   <th bgcolor=DDDDDD >GPP Glux</th>" \

   487        ,"   <th bgcolor=DDDDDD >Respiration</th>" \

   488        ,"   <th bgcolor=DDDDDD >Average</th>" \

   489        ,"</tr>" \

   490        /)

   491   table_footer = "</table>"

   492   row_header = "<tr>"

   493   row_footer = "</tr>"

   495   lines = new(50000,string)

   496   nline = 0

   498   set_line(lines,nline,header)

   499   set_line(lines,nline,table_header)

   500 ;-----------------------------------------------

   501 ; row of table

   503   do n = 0,nstation-1

   504      set_line(lines,nline,row_header)

   506      txt0  = station(n)

   507      txt1  = sprintf("%5.2f", lat_ob(n))

   508      txt2  = sprintf("%5.2f", lon_ob(n))

   509      txt3  = year_ob(n)

   510      txt4  = sprintf("%5.2f", M_score(n,0))

   511      txt5  = sprintf("%5.2f", M_score(n,1))

   512      txt6  = sprintf("%5.2f", M_score(n,2))

   513      txt7  = sprintf("%5.2f", M_score(n,3))

   514      txt8  = sprintf("%5.2f", M_score(n,4))

   515      txt9  = sprintf("%5.2f", M_score(n,5))

   516      txt10 = sprintf("%5.2f", avg(M_score(n,:)))

   518      set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")

   519      set_line(lines,nline,"<th>"+txt1+"</th>")

   520      set_line(lines,nline,"<th>"+txt2+"</th>")

   521      set_line(lines,nline,"<th>"+txt3+"</th>")

   522      set_line(lines,nline,"<th>"+txt4+"</th>")

   523      set_line(lines,nline,"<th>"+txt5+"</th>")

   524      set_line(lines,nline,"<th>"+txt6+"</th>")

   525      set_line(lines,nline,"<th>"+txt7+"</th>")

   526      set_line(lines,nline,"<th>"+txt8+"</th>")

   527      set_line(lines,nline,"<th>"+txt9+"</th>")

   528      set_line(lines,nline,"<th>"+txt10+"</th>")

   530      set_line(lines,nline,row_footer)

   531   end do

   533 ; last row, summary

   534   set_line(lines,nline,row_header)

   536   txt0  = "All_"+sprintf("%.0f", nstation)

   537   txt1  = "-"

   538   txt2  = "-"

   539   txt3  = "-"

   540   txt4  = M_energy_co2

   541   txt5  = M_energy_rad

   542   txt6  = M_energy_lh

   543   txt7  = M_energy_sh

   544   txt8  = M_energy_gpp

   545   txt9  = M_energy_er

   546   txt10 = M_energy_all

   548   set_line(lines,nline,"<th>"+txt0+"</th>")

   549   set_line(lines,nline,"<th>"+txt1+"</th>")

   550   set_line(lines,nline,"<th>"+txt2+"</th>")

   551   set_line(lines,nline,"<th>"+txt3+"</th>")

   552   set_line(lines,nline,"<th>"+txt4+"</th>")

   553   set_line(lines,nline,"<th>"+txt5+"</th>")

   554   set_line(lines,nline,"<th>"+txt6+"</th>")

   555   set_line(lines,nline,"<th>"+txt7+"</th>")

   556   set_line(lines,nline,"<th>"+txt8+"</th>")

   557   set_line(lines,nline,"<th>"+txt9+"</th>")

   558   set_line(lines,nline,"<th>"+txt10+"</th>")

   560   set_line(lines,nline,row_footer)

   561 ;-----------------------------------------------

   562   set_line(lines,nline,table_footer)

   563   set_line(lines,nline,footer)

   565 ; Now write to an HTML file.

   566   idx = ind(.not.ismissing(lines))

   567   if(.not.any(ismissing(idx))) then

   568     asciiwrite(output_html,lines(idx))

   569   else

   570    print ("error?")

   571   end if

   572   delete (idx)

   574 ;***************************************************************************

   575 ; output plots

   576 ;***************************************************************************

   577   output_dir = model_name+"/ameriflux"

   579   system("mv *.png *.html " + output_dir)

   580 ;***************************************************************************

   581 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--