clamp: ameriflux/10.plot.ncl@0c6405ab2ff4

     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 = (/"HarvardForest" \

    70             /)

    72  year_ob = (/"1991-2004" \

    73             /)

    75  field   = (/"CO2 Flux" \

    76             ,"Net Radiation" \

    77             ,"Latent Heat" \

    78             ,"Sensible Heat" \

    79             ,"GPP Flux" \

    80             ,"Respiration" \

    81             /)

    83  nstation  = dimsizes(station)

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

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

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

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

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

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

    93  do n = 0,nstation-1

    95 ;-------------------------------------------------

    96 ;   get ob data

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

    99     filo = station(n)+"_"+year_ob(n)+"_monthly.nc"

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

   102     print (filo)

   104     lon_ob(n) = fo->lon

   105     lat_ob(n) = fo->lat

   107     data      = fo->NEE_or_fMDS

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

   109     delete (data)

   111     data      = fo->Rg_f

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

   113     delete (data)

   115     data      = fo->LE_f

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

   117     delete (data)

   119     data      = fo->H_f

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

   121     delete (data)

   123     data      = fo->GPP_or_MDS

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

   125     delete (data)

   127     data      = fo->Reco_or

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

   129     delete (data)

   131     delete (fo)

   132 ;---------------------------------------------------

   133 ;   get model data

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

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

   138     print (film)

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

   142   data = fm->NEE

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

   144   delete (data)

   146 ; data  = fm->LATENT

   147   data1 = fm->FCEV

   148   data2 = fm->FCTR

   149   data3 = fm->FGEV

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

   151   delete (data1)

   152   delete (data2)

   153   delete (data3)

   155 ; data = fm->SENSIBLE

   156   data  = fm->FSH

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

   158   delete (data)

   160 ; data  = fm->NETRAD

   161   data1 = fm->FSA

   162   data2 = fm->FIRA

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

   164   delete (data1)

   165   delete (data2)

   167 else

   169   data = fm->NEE

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

   171   delete (data)

   173   data = fm->NETRAD

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

   175   delete (data)

   177   data = fm->LATENT

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

   179   delete (data)

   181 ; data = fm->SENSIBLE

   182   data = fm->FSH

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

   184   delete (data)

   185 end if

   187   data = fm->GPP

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

   189   delete (data)

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

   192      data = fm->ER

   193   else

   194      data1 = fm->AR

   195      data2 = fm->HR

   196      data  = data1 + data2

   198      delete (data1)

   199      delete (data2)

   200   end if

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

   203   delete (data)

   205   delete (fm)

   207 ;************************************************************

   208 ; interpolate model data into observed station

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

   210 ;************************************************************

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

   213   ym(0) = -90.

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

   217 ; print (yy)

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

   221   delete (yy)

   223  end do

   224 ;************************************************************

   225 ; compute correlation coef and M score

   226 ;************************************************************

   228  score_max = 5.

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

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

   233  do n=0,nstation-1

   234  do m=0,nfield-1

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

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

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

   238  end do

   239  end do

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   257 ;*******************************************************************

   258 ; for station line plot

   259 ;*******************************************************************

   261 ; for x-axis in xyplot

   262   mon = ispan(1,12,1)

   263   mon@long_name = "month"

   265   res                   = True               ; plot mods desired

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

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

   268 ;-------------------------------------------------------------------------

   269 ; Add a boxed legend using the more simple method

   271   res@pmLegendDisplayMode    = "Always"

   272 ; res@pmLegendWidthF         = 0.1

   273   res@pmLegendWidthF         = 0.08

   274   res@pmLegendHeightF        = 0.06

   275 ; res@pmLegendOrthogonalPosF = -1.17

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

   277   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   279 ; res@pmLegendParallelPosF   =  0.18

   280   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   282 ; res@lgPerimOn             = False

   283   res@lgLabelFontHeightF     = 0.015

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

   285 ;-------------------------------------------------------------------

   286 ; for panel plot

   287   res@gsnFrame     = False                   ; Do not draw plot

   288   res@gsnDraw      = False                   ; Do not advance frame

   290   pres                            = True     ; panel plot mods desired

   291   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   292                                              ; indiv. plots in panel

   293   pres@gsnMaximize                = True     ; fill the page

   294 ;-------------------------------------------------------------------

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

   297   plot_data!0 = "case"

   298   plot_data!1 = "month"

   300   do n = 0,nstation-1

   301 ;----------------------------

   302 ; for observed

   304     plot_name = station(n)+"_ob"

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

   306     res@tiMainString = title

   308     wks = gsn_open_wks (plot_type,plot_name)

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

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

   312     plot_data@long_name = field(0)

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

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

   316     plot_data@long_name = field(1)

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

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

   320     plot_data@long_name = field(2)

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

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

   324     plot_data@long_name = field(3)

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

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

   328     plot_data@long_name = field(4)

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

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

   332     plot_data@long_name = field(5)

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

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

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

   338            "rm "+plot_name+"."+plot_type)

   340     clear (wks)

   341     delete (plot)

   342 ;----------------------------

   343 ; for model_vs_ob

   345     plot_name = station(n)+"_model_vs_ob"

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

   347     res@tiMainString = title

   349     wks = gsn_open_wks (plot_type,plot_name)

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

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

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

   354     plot_data@long_name = field(0)

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

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

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

   359     plot_data@long_name = field(1)

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

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

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

   364     plot_data@long_name = field(2)

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

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

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

   369     plot_data@long_name = field(3)

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

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

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

   375     plot_data@long_name = field(4)

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

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

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

   381     plot_data@long_name = field(5)

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

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

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

   387            "rm "+plot_name+"."+plot_type)

   389     clear (wks)

   390     delete (plot)

   391  end do

   393 ;*******************************************************************

   394 ; html table of site: observed

   395 ;*******************************************************************

   396   output_html = "line_ob.html"

   398   header = (/"<HTML>" \

   399             ,"<HEAD>" \

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

   401             ,"</HEAD>" \

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

   403             /)

   404   footer = "</HTML>"

   406   table_header = (/ \

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

   408        ,"<tr>" \

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

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

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

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

   413        ,"</tr>" \

   414        /)

   415   table_footer = "</table>"

   416   row_header = "<tr>"

   417   row_footer = "</tr>"

   419   lines = new(50000,string)

   420   nline = 0

   422   set_line(lines,nline,header)

   423   set_line(lines,nline,table_header)

   424 ;-----------------------------------------------

   425 ; row of table

   427   do n = 0,nstation-1

   428      set_line(lines,nline,row_header)

   430      txt0 = station(n)

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

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

   433      txt3 = year_ob(n)

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

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

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

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

   440      set_line(lines,nline,row_footer)

   441   end do

   442 ;-----------------------------------------------

   443   set_line(lines,nline,table_footer)

   444   set_line(lines,nline,footer)

   446 ; Now write to an HTML file.

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

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

   449     asciiwrite(output_html,lines(idx))

   450   else

   451    print ("error?")

   452   end if

   453   delete (idx)

   455 ;*******************************************************************

   456 ; score and line table : model vs observed

   457 ;*******************************************************************

   458   output_html = "score+line_vs_ob.html"

   460   header = (/"<HTML>" \

   461             ,"<HEAD>" \

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

   463             ,"</HEAD>" \

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

   465             /)

   466   footer = "</HTML>"

   468   delete (table_header)

   469   table_header = (/ \

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

   471        ,"<tr>" \

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

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

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

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

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

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

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

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

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

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

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

   483        ,"</tr>" \

   484        /)

   485   table_footer = "</table>"

   486   row_header = "<tr>"

   487   row_footer = "</tr>"

   489   lines = new(50000,string)

   490   nline = 0

   492   set_line(lines,nline,header)

   493   set_line(lines,nline,table_header)

   494 ;-----------------------------------------------

   495 ; row of table

   497   do n = 0,nstation-1

   498      set_line(lines,nline,row_header)

   500      txt0  = station(n)

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

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

   503      txt3  = year_ob(n)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   524      set_line(lines,nline,row_footer)

   525   end do

   527 ; last row, summary

   528   set_line(lines,nline,row_header)

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

   531   txt1  = "-"

   532   txt2  = "-"

   533   txt3  = "-"

   534   txt4  = M_energy_co2

   535   txt5  = M_energy_rad

   536   txt6  = M_energy_lh

   537   txt7  = M_energy_sh

   538   txt8  = M_energy_gpp

   539   txt9  = M_energy_er

   540   txt10 = M_energy_all

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

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

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

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

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

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

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

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

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

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

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

   554   set_line(lines,nline,row_footer)

   555 ;-----------------------------------------------

   556   set_line(lines,nline,table_footer)

   557   set_line(lines,nline,footer)

   559 ; Now write to an HTML file.

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

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

   562     asciiwrite(output_html,lines(idx))

   563   else

   564    print ("error?")

   565   end if

   566   delete (idx)

   568 ;***************************************************************************

   569 ; output plots

   570 ;***************************************************************************

   571   output_dir = model_name+"/ameriflux"

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

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

   575 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
permissions	-rw-r--r--