Mercurial Mercurial > repos > public_software > clamp / file revision
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
ameriflux/21.plot.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB. 
     1 ;************************************************************

     2 ; add table

     3 ; add time series plot

     4 ; sort by latitude in decending order (N->S)

     5 ; find year_ob >= 4

     6 ; add year_ob_i and year_ob_f

     7 ; change long_name

     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

    16     

    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

    29 

    30   plot_type     = "ps"

    31   plot_type_new = "png"

    32 

    33 ;************************************************

    34 ; model data info

    35 ;************************************************

    36   model = "cn"

    37 ; model = "casa"

    38   

    39   model_name = "i01.10" + model

    40   

    41   ENERGY ="new"

    42 

    43 ;************************************************

    44 ; observed data info

    45 ;************************************************

    46 

    47  station = (/"ARM_Oklahoma" \

    48             ,"ARM_Oklahoma_burn" \

    49             ,"ARM_Oklahoma_control" \

    50             ,"Atqasuk" \

    51             ,"Audubon" \

    52             ,"AustinCary" \

    53             ,"Bartlett" \

    54             ,"Bondville" \

    55             ,"Brookings" \

    56             ,"Donaldson" \

    57             ,"Duke_Forest_Hardwoods" \

    58             ,"Duke_Forest_Open_Field" \

    59             ,"Duke_Forest_Pine" \

    60             ,"Fermi_Ag" \

    61             ,"Fermi_Prairie" \

    62             ,"Flagstaff_Managed" \

    63             ,"Flagstaff_Unmanaged" \

    64             ,"Flagstaff_Wildfire" \

    65             ,"FortPeck" \

    66             ,"FreemanRanch_mesquite" \

    67             ,"Goodwin_Creek" \

    68             ,"HarvardForest" \

    69             ,"HarvardForestHemlock" \

    70             ,"HowlandForestMain" \

    71             ,"HowlandForestWest" \

    72             ,"Ivotuk" \

    73             ,"KendallGrasslands" \

    74             ,"KennedySpaceCenterPine" \

    75             ,"KennedySpaceCenterScrub" \

    76             ,"LittleProspect" \

    77             ,"LostCreek" \

    78             ,"Mead-irrigated" \

    79             ,"Mead-irrigated-rotation" \

    80             ,"Mead-rainfed" \

    81             ,"Metolius_2nd_YoungPonderosaPine" \

    82             ,"MetoliusEyerly" \

    83             ,"MetoliusIntermediatePine" \

    84             ,"MetoliusOldPonderosaPine" \

    85             ,"MissouriOzark" \

    86             ,"Mize" \

    87             ,"MorganMonroe" \

    88             ,"NiwotRidge" \

    89             ,"NorthCarolina_cc" \

    90             ,"NorthCarolina_lp" \

    91             ,"ParkFalls" \

    92             ,"Rayonier" \

    93             ,"SantaRita" \

    94             ,"SkyOaks_Old" \

    95             ,"SkyOaks_PostFire" \

    96             ,"SkyOaks_Young" \

    97             ,"SylvaniaWilderness" \

    98             ,"Toledo" \

    99             ,"Tonzi" \

   100             ,"UCI_1850" \

   101             ,"UCI_1930" \

   102             ,"UCI_1964" \

   103             ,"UCI_1964wet" \

   104             ,"UCI_1981" \

   105             ,"UCI_1989" \

   106             ,"UCI_1998" \

   107             ,"UMBS" \

   108             ,"Vaira" \

   109             ,"WalkerBranch" \

   110             ,"WillowCreek" \

   111             ,"WindRiver" \

   112             ,"Wisconsin_ihw" \

   113             ,"Wisconsin_irp" \

   114             ,"Wisconsin_mrp" \

   115             ,"Wisconsin_myjp" \

   116             ,"Wisconsin_pb" \

   117             ,"Wisconsin_rpcc" \

   118             ,"Wisconsin_yhw" \

   119             ,"Wisconsin_yjp" \

   120             ,"Wisconsin_yrp" \

   121             /)

   122 

   123  year_ob = (/"2003-2006" \

   124             ,"2005-2006" \

   125             ,"2005-2006" \

   126             ,"1999-2006" \

   127             ,"2002-2006" \

   128             ,"2000-2005" \

   129             ,"2004-2005" \

   130             ,"1996-2006" \

   131             ,"2004-2006" \

   132             ,"1999-2004" \

   133             ,"2003-2005" \

   134             ,"2001-2005" \

   135             ,"2001-2005" \

   136             ,"2005-2006" \

   137             ,"2004-2006" \

   138             ,"2005-2006" \

   139             ,"2005-2006" \

   140             ,"2005-2006" \

   141             ,"2000-2006" \

   142             ,"2004-2006" \

   143             ,"2002-2006" \

   144             ,"1991-2004" \

   145             ,"2004-2004" \

   146             ,"1996-2004" \

   147             ,"1999-2004" \

   148             ,"2003-2006" \

   149             ,"2004-2006" \

   150             ,"2002-2002" \

   151             ,"2000-2006" \

   152             ,"2002-2005" \

   153             ,"2001-2005" \

   154             ,"2001-2005" \

   155             ,"2001-2005" \

   156             ,"2001-2005" \

   157             ,"2004-2005" \

   158             ,"2004-2005" \

   159             ,"2003-2005" \

   160             ,"1996-2000" \

   161             ,"2004-2006" \

   162             ,"1998-2004" \

   163             ,"1999-2005" \

   164             ,"1999-2003" \

   165             ,"2005-2006" \

   166             ,"2005-2006" \

   167             ,"1996-2003" \

   168             ,"1998-1998" \

   169             ,"2004-2006" \

   170             ,"1997-2006" \

   171             ,"2004-2006" \

   172             ,"1997-2006" \

   173             ,"2002-2006" \

   174             ,"2004-2005" \

   175             ,"2001-2006" \

   176             ,"2004-2005" \

   177             ,"2001-2005" \

   178             ,"2001-2005" \

   179             ,"2002-2004" \

   180             ,"2001-2005" \

   181             ,"2001-2005" \

   182             ,"2002-2005" \

   183             ,"1999-2003" \

   184             ,"2001-2006" \

   185             ,"1995-1999" \

   186             ,"1999-2005" \

   187             ,"1999-2004" \

   188             ,"2003-2003" \

   189             ,"2003-2003" \

   190             ,"2002-2005" \

   191             ,"2004-2004" \

   192             ,"2002-2002" \

   193             ,"2005-2005" \

   194             ,"2002-2002" \

   195             ,"2004-2005" \

   196             ,"2002-2002" \

   197             /)

   198 

   199  field   = (/"NEE Flux" \

   200             ,"Net Radiation" \

   201             ,"Latent Heat" \

   202             ,"Sensible Heat" \

   203             ,"GPP Flux" \

   204             ,"Respiration" \

   205             /)

   206 

   207  field_unit = (/"u mol/m2/s" \

   208                ,"W/m2" \

   209                ,"W/m2" \

   210                ,"W/m2" \

   211                ,"u mol/m2/s" \

   212                ,"u mol/m2/s" \

   213                /)

   214 

   215  nstation = dimsizes(station)

   216  nfield   = dimsizes(field)

   217 

   218 ;========================================================================

   219 ; get observed info: number of year, first/last year, lat, lon

   220 

   221  dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"

   222 

   223  year_station = new ((/nstation/),integer) ; number of year

   224  year_ob_i    = new ((/nstation/),integer) ; first year

   225  year_ob_f    = new ((/nstation/),integer) ; last year

   226  lat_ob       = new ((/nstation/),float)   ; latitude

   227  lon_ob       = new ((/nstation/),float)   ; longitude 

   228 

   229  do n = 0, nstation-1

   230 

   231     diro = dir_root + station(n)+"/"

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

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

   234  

   235     lat_ob(n) = fo->lat

   236     lon_ob(n) = fo->lon

   237 

   238     year      = fo->year

   239 

   240     year_station(n) = dimsizes(year)

   241     year_ob_i(n)    = year(0)

   242     year_ob_f(n)    = year(year_station(n)-1)

   243  

   244     delete (year) 

   245     

   246  end do

   247 

   248  delete (fo)

   249 ;print (year_ob_i)

   250 ;print (year_ob_f)

   251 

   252 ;--------------------------------------------------------------

   253 ; find (# of year observed) >=4 and year_ob_i <= 2001

   254 

   255  i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)

   256  

   257  station_long      = station(i_long_ob)

   258  lat_ob_long       = lat_ob(i_long_ob)

   259  lon_ob_long       = lat_ob(i_long_ob)

   260  year_ob_long      = year_ob(i_long_ob)

   261  year_ob_i_long    = year_ob_i(i_long_ob)

   262  year_ob_f_long    = year_ob_f(i_long_ob)

   263  year_station_long = year_station(i_long_ob)

   264 

   265  nstation_long     = dimsizes(station_long)

   266 

   267 ; print (i_long_ob)

   268 ; print (nstation_long)

   269 ; print (station_long)

   270 ; print (year_ob_i(i_long_ob))

   271 

   272 ;=========================================================

   273 ;   get model data at observed lat-lon

   274 

   275     dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"    

   276     film = model_name+"_ameriflux_1990-2004_monthly.nc"

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

   278 

   279     xm   = fm->lon

   280     ym   = fm->lat 

   281     date = fm->date

   282 

   283     date_dim = dimsizes(date)

   284     nyear  = date_dim(0)

   285 

   286     nmonth =12

   287 

   288     data_mod      = new ((/nfield,nyear,nmonth,nstation/),float)

   289     data_mod_ann  = new ((/nfield,nmonth,nstation/),float)

   290     data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)

   291 

   292 ;   change to unit of observed (u mol/m2/s)

   293 ;   Model_units [=] gC/m2/s

   294 ;   12. = molecular weight of C

   295 ;   u mol = 1e-6 mol

   296 

   297     factor = 1.e6 /12.

   298 

   299 ;------------------------------------------------------------

   300 ;   interpolate model data into observed station

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

   302 

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

   304     ym(0) = -90.

   305 ;------------------------------------------------------------  

   306 

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

   308 

   309   data = fm->NEE

   310   data_mod(0,:,:,:) = data(:,:,:) * factor

   311 

   312 ; data  = fm->LATENT

   313   data1 = fm->FCEV

   314   data2 = fm->FCTR

   315   data3 = fm->FGEV

   316   data_mod(2,:,:,:) = data1(:,:,:)+data2(:,:,:)+data3(:,:,:) 

   317   delete (data1)

   318   delete (data2)

   319   delete (data3)

   320 

   321 ; data = fm->SENSIBLE

   322   data  = fm->FSH

   323   data_mod(3,:,:,:) = data(:,:,:) 

   324 

   325 ; data  = fm->NETRAD

   326   data1 = fm->FSA

   327   data2 = fm->FIRA

   328   data_mod(1,:,:,:) = data1(:,:,:)-data2(:,:,:) 

   329   delete (data1)

   330   delete (data2)

   331 

   332 else

   333 

   334   data = fm->NEE

   335   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   336   data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor

   337   data_mod(0,:,:,:) = yy(:,:,:) * factor

   338 

   339 ; printVarSummary(yy)

   340 

   341   data = fm->NETRAD

   342   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   343   data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   344   data_mod(1,:,:,:) = yy(:,:,:) 

   345 

   346   data = fm->LATENT

   347   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   348   data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   349   data_mod(2,:,:,:) = yy(:,:,:) 

   350 

   351 ; data = fm->SENSIBLE

   352   data = fm->FSH

   353   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   354   data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 

   355   data_mod(3,:,:,:) = yy(:,:,:) 

   356 

   357 end if

   358 

   359   data = fm->GPP

   360   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   361   data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor

   362   data_mod(4,:,:,:) = yy(:,:,:) * factor 

   363 

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

   365      data = fm->ER

   366   else

   367      data1 = fm->AR

   368      data2 = fm->HR

   369      data  = data1 + data2

   370     

   371      delete (data1)

   372      delete (data2)

   373   end if

   374 

   375   yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)

   376   data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor

   377   data_mod(5,:,:,:) = yy(:,:,:) * factor

   378 

   379   data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)

   380 

   381   delete (data_mod)

   382   delete (fm) 

   383   delete (data)

   384   delete (yy)

   385 

   386 ;*******************************************************************

   387 ; for station line plot

   388 ;*******************************************************************

   389 

   390 ; for x-axis in xyplot

   391   mon = ispan(1,12,1)

   392   mon@long_name = "month"

   393 

   394   res                   = True               ; plot mods desired

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

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

   397 

   398   res@tmXBFormat  = "f"                      ; not to add trailing zeros

   399 

   400 ;-------------------------------------------------------------------------

   401 ; Add a boxed legend using the more simple method

   402 

   403   res@pmLegendDisplayMode    = "Always"

   404 ; res@pmLegendWidthF         = 0.1

   405   res@pmLegendWidthF         = 0.08

   406   res@pmLegendHeightF        = 0.06

   407 ; res@pmLegendOrthogonalPosF = -1.17

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

   409   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   410 

   411 ; res@pmLegendParallelPosF   =  0.18

   412   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   413 

   414 ; res@lgPerimOn             = False

   415   res@lgLabelFontHeightF     = 0.015

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

   417 ;-------------------------------------------------------------------

   418 ; for panel plot

   419   res@gsnFrame     = False                   ; Do not draw plot 

   420   res@gsnDraw      = False                   ; Do not advance frame

   421 

   422   pres                            = True     ; panel plot mods desired

   423   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   424                                              ; indiv. plots in panel

   425   pres@gsnMaximize                = True     ; fill the page

   426 ;-------------------------------------------------------------------

   427 

   428 ;==============================================================

   429 ; get ob data at each site with long observation

   430 

   431   dir_root = "/fis/cgd/cseg/people/jeff/clamp/ameriflux/data/"

   432 

   433   do n = 0,nstation_long-1

   434 ; do n = 0,0

   435 

   436     year_setback = 0

   437 

   438     nyear = year_station_long(n)

   439 

   440     if (year_ob_f_long(n).eq. 2006) then

   441        year_setback = 2006 -2004

   442     end if

   443     if (year_ob_f_long(n).eq. 2005) then

   444        year_setback = 2005 -2004

   445     end if

   446     

   447     ntime = (nyear - year_setback) * nmonth

   448     print (nyear)

   449 

   450     data_ob   = new ((/nfield, nyear, nmonth/),float)

   451 

   452     diro = dir_root + station_long(n)+"/"

   453     filo = year_ob_long(n)+"_L4_m.nc"

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

   455 

   456 ;   print (diro)

   457 ;   print (filo)

   458 

   459     data_ob(0,:,:) = fo->NEE_or_fMDS

   460     data_ob(1,:,:) = fo->Rg_f

   461     data_ob(2,:,:) = fo->LE_f

   462     data_ob(3,:,:) = fo->H_f

   463     data_ob(4,:,:) = fo->GPP_or_MDS

   464     data_ob(5,:,:) = fo->Reco_or

   465 

   466     delete (fo)

   467 

   468     timeI = new((/ntime/),integer)

   469     timeF = new((/ntime/),float)

   470     timeI = ispan(1,ntime,1)

   471     timeF = year_ob_i_long(n) + (timeI-1)/12.

   472     timeF@long_name = "year" 

   473 

   474     plot_data = new((/2,ntime/),float)

   475  

   476 ;----------------------------

   477 ; for model_vs_ob

   478 

   479     plot_name = station_long(n)+"_tseries_vs_ob"

   480     title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"    

   481     res@tiMainString = title

   482 

   483     wks = gsn_open_wks (plot_type,plot_name)

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

   485 

   486     i_year_mod_i = year_ob_i_long(n) - 1990 

   487     i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback

   488 

   489     i_year_ob_f =  nyear - year_setback - 1 

   490 

   491 ;   print (nyear)

   492 ;   print (i_year_ob_f)

   493 ;   print (i_year_mod_i)

   494 ;   print (i_year_mod_f)

   495 

   496     do i = 0,nfield-1                           

   497        plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))

   498        plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))

   499        plot_data@long_name = field(i)+" ("+field_unit(i)+")"   

   500        plot(i)=gsn_csm_xy(wks,timeF,plot_data,res)            ; create plot 

   501     end do

   502    

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

   504 

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

   506            "rm "+plot_name+"."+plot_type)

   507 

   508     clear (wks)  

   509     delete (plot)

   510 

   511     delete (data_ob)

   512     delete (timeI)    

   513     delete (timeF)

   514     delete (plot_data)

   515 

   516  end do

   517 

   518 ;###################################################################

   519 ; for the following tables,

   520 ; sort by latitude in decending order (N->S)

   521 

   522   isort = dim_pqsort(lat_ob_long,-1)

   523 

   524   station_sort = station_long(isort)

   525   year_ob_sort = year_ob_long(isort)

   526   lat_ob_sort  = lat_ob_long(isort)

   527   lon_ob_sort  = lon_ob_long(isort)

   528  

   529 ; print(isort)

   530 ; print(lat_ob_sort)

   531 

   532 ;*******************************************************************

   533 ; html table of site: observed

   534 ;*******************************************************************

   535   output_html = "tseries_vs_ob.html"

   536 

   537   header = (/"<HTML>" \

   538             ,"<HEAD>" \

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

   540             ,"</HEAD>" \

   541             ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \

   542             /) 

   543   footer = "</HTML>"

   544 

   545   table_header = (/ \

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

   547        ,"<tr>" \

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

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

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

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

   552        ,"</tr>" \

   553        /)

   554   table_footer = "</table>"

   555   row_header = "<tr>"

   556   row_footer = "</tr>"

   557 

   558   lines = new(50000,string)

   559   nline = 0

   560 

   561   set_line(lines,nline,header)

   562   set_line(lines,nline,table_header)

   563 ;-----------------------------------------------

   564 ; row of table

   565   

   566   do n = 0,nstation_long-1

   567 

   568      set_line(lines,nline,row_header)

   569 

   570      txt0 = station_sort(n)

   571      txt1 = sprintf("%5.2f", lat_ob_sort(n))

   572      txt2 = sprintf("%5.2f", lon_ob_sort(n))

   573      txt3 = year_ob_sort(n)

   574 

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

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

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

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

   579 

   580      set_line(lines,nline,row_footer)

   581   end do

   582 ;-----------------------------------------------

   583   set_line(lines,nline,table_footer)

   584   set_line(lines,nline,footer) 

   585 

   586 ; Now write to an HTML file.

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

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

   589     asciiwrite(output_html,lines(idx))

   590   else

   591    print ("error?")

   592   end if

   593   delete (idx)

   594 

   595 ;***************************************************************************

   596 ; output plots

   597 ;***************************************************************************

   598   output_dir = model_name+"/ameriflux"

   599 

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

   601 ;***************************************************************************

   602 end
clamp