Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;********************************************************
2 ; required command line input parameters:
3 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
5 ; histogram normalized by rain and compute correleration
6 ;**************************************************************
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
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 ;**************************************************************
32 ;************************************************
34 ;************************************************
40 model_name_i = "i01.07cn"
41 model_name_f = "i01.10cn"
43 ;model_name_i = "i01.07casa"
44 ;model_name_f = "i01.10casa"
46 model_name = model_name_f
48 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
49 film_i = model_name_i + "_1990-2004_ANN_climo.nc"
50 film_f = model_name_f + "_1990-2004_ANN_climo.nc"
52 fm_i = addfile (dirm+film_i,"r")
53 fm_f = addfile (dirm+film_f,"r")
61 ;Units for these variables are:
64 nsec_per_year = 60*60*24*365
66 npp_i = npp_i * nsec_per_year
67 npp_f = npp_f * nsec_per_year
69 ;===================================================
70 ; read data: observed -station
71 ;===================================================
73 station = (/"DukeFACE" \
79 lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
80 lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
81 lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
84 n_sta = dimsizes(station)
85 beta_4_ob = new((/n_sta/),float)
88 ;===================================================
89 ; get model data at station
90 ;===================================================
92 npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
94 npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
99 ;============================
101 ;============================
103 beta_4 = new((/n_sta/),float)
105 beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
107 beta_4_avg = avg(beta_4)
112 ;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
114 bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
115 M_beta = (1. - (bias/n_sta))*3.
119 ;=========================
120 ; for html table - station
121 ;=========================
123 output_html = "table_station.html"
125 ; column (not including header column)
127 col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
129 ncol = dimsizes(col_head)
131 ; row (not including header row)
132 row_head = (/"DukeFACE" \
138 nrow = dimsizes(row_head)
140 ; arrays to be passed to table.
141 text4 = new ((/nrow, ncol/),string )
144 text4(i,0) = sprintf("%.1f",lat_ob(i))
145 text4(i,1) = sprintf("%.1f",lon_ob(i))
146 text4(i,2) = sprintf("%.1f",co2_i)
147 text4(i,3) = sprintf("%.1f",co2_f)
148 text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
149 text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
150 text4(i,6) = sprintf("%.2f",beta_4(i))
153 text4(nrow-1,0) = "-"
154 text4(nrow-1,1) = "-"
155 text4(nrow-1,2) = "-"
156 text4(nrow-1,3) = "-"
157 text4(nrow-1,4) = "-"
158 text4(nrow-1,5) = "-"
159 text4(nrow-1,6) = sprintf("%.2f",beta_4_avg)
160 text4(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
166 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
168 header = (/"<HTML>" \
170 ,"<TITLE>CLAMP metrics</TITLE>" \
177 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
179 ," <th bgcolor=DDDDDD >Station</th>" \
180 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
181 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
182 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
183 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
184 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
185 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
186 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
187 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
190 table_footer = "</table>"
194 lines = new(50000,string)
197 set_line(lines,nline,header)
198 set_line(lines,nline,table_header)
199 ;-----------------------------------------------
203 set_line(lines,nline,row_header)
215 set_line(lines,nline,"<th>"+txt1+"</th>")
216 set_line(lines,nline,"<th>"+txt2+"</th>")
217 set_line(lines,nline,"<th>"+txt3+"</th>")
218 set_line(lines,nline,"<th>"+txt4+"</th>")
219 set_line(lines,nline,"<th>"+txt5+"</th>")
220 set_line(lines,nline,"<th>"+txt6+"</th>")
221 set_line(lines,nline,"<th>"+txt7+"</th>")
222 set_line(lines,nline,"<th>"+txt8+"</th>")
223 set_line(lines,nline,"<th>"+txt9+"</th>")
225 set_line(lines,nline,row_footer)
227 ;-----------------------------------------------
228 set_line(lines,nline,table_footer)
229 set_line(lines,nline,footer)
231 ; Now write to an HTML file.
232 idx = ind(.not.ismissing(lines))
233 if(.not.any(ismissing(idx))) then
234 asciiwrite(output_html,lines(idx))
242 delete (table_header)
245 ;************************************************
246 ; read data: observed-2
247 ;************************************************
249 ob_name = "MODIS MOD 15A2 2000-2005"
251 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
252 filo = "land_class_"+model_grid+".nc"
254 fo = addfile(diro+filo,"r")
256 classob = tofloat(fo->LAND_CLASS)
258 ; observed data has 20 land-type classes
261 ;*******************************************************************
262 ; Calculate "nice" bins for binning the data in equally spaced ranges
263 ;********************************************************************
266 range = fspan(0,nclassn-1,nclassn)
269 ; Use this range information to grab all the values in a
270 ; particular range, and then take an average.
274 xvalues = new((/2,nx/),float)
275 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
276 dx = xvalues(0,1) - xvalues(0,0) ; range width
277 dx4 = dx/4 ; 1/4 of the range
278 xvalues(1,:) = xvalues(0,:) - dx/5.
281 base_1D = ndtooned(classob)
282 data1_1D = ndtooned(npp_i)
283 data2_1D = ndtooned(npp_f)
287 yvalues = new((/2,nx/),float)
288 count = new((/2,nx/),float)
292 ; See if we are doing data1 (nd=0) or data2 (nd=1).
302 ; Loop through each range, using base.
305 if (i.ne.(nr-2)) then
307 ; print("In range ["+range(i)+","+range(i+1)+")")
308 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
311 ; print("In range ["+range(i)+",)")
312 idx = ind(base.ge.range(i))
317 if(.not.any(ismissing(idx))) then
318 yvalues(nd,i) = avg(data(idx))
319 count(nd,i) = dimsizes(idx)
321 yvalues(nd,i) = yvalues@_FillValue
325 ;#############################################################
326 ; set the following 4 classes to _FillValue:
327 ; Water Bodies(0), Urban and Build-Up(13),
328 ; Permenant Snow and Ice(15), Unclassified(17)
330 if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
331 yvalues(nd,i) = yvalues@_FillValue
334 ;#############################################################
336 ; print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
338 ; Clean up for next time in loop.
346 ;============================
348 ;============================
355 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
359 uu_count = u_count(good)
360 vv_count = v_count(good)
362 n_biome = dimsizes(uu)
363 beta_biome = new((/n_biome/),float)
365 beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
367 ;beta_biome_avg = avg(beta_biome)
368 beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
370 print (beta_biome_avg)
372 ;===========================
373 ; for html table - biome
374 ;===========================
376 output_html = "table_biome.html"
378 ; column (not including header column)
380 col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
382 ncol = dimsizes(col_head)
384 ; row (not including header row)
385 ; 3 classes removed: Water Bodies, Urban and Build-Up, Unclassified
386 ; function "good" removed the last 2 classes
387 ; text4(i,2) = sprintf("%.2f",uu(i+1)) remove the first class
389 row_head = (/"Evergreen Needleleaf Forests" \
390 ,"Evergreen Broadleaf Forests" \
391 ,"Deciduous Needleleaf Forest" \
392 ,"Deciduous Broadleaf Forests" \
394 ,"Closed Bushlands" \
396 ,"Woody Savannas (S. Hem.)" \
397 ,"Savannas (S. Hem.)" \
399 ,"Permanent Wetlands" \
401 ,"Cropland/Natural Vegetation Mosaic" \
402 ,"Barren or Sparsely Vegetated" \
403 ,"Woody Savannas (N. Hem.)" \
404 ,"Savannas (N. Hem.)" \
407 nrow = dimsizes(row_head)
409 ; arrays to be passed to table.
410 text4 = new ((/nrow, ncol/),string )
413 text4(i,0) = sprintf("%.1f",co2_i)
414 text4(i,1) = sprintf("%.1f",co2_f)
415 text4(i,2) = sprintf("%.1f",uu(i))
416 text4(i,3) = sprintf("%.1f",vv(i))
417 text4(i,4) = sprintf("%.2f",beta_biome(i))
419 text4(nrow-1,0) = "-"
420 text4(nrow-1,1) = "-"
421 text4(nrow-1,2) = "-"
422 text4(nrow-1,3) = "-"
423 text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
425 ;**************************************************
427 ;**************************************************
429 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
431 header = (/"<HTML>" \
433 ,"<TITLE>CLAMP metrics</TITLE>" \
440 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
442 ," <th bgcolor=DDDDDD >Biome Class</th>" \
443 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
444 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
445 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
446 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
447 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
450 table_footer = "</table>"
454 lines = new(50000,string)
457 set_line(lines,nline,header)
458 set_line(lines,nline,table_header)
459 ;-----------------------------------------------
463 set_line(lines,nline,row_header)
472 set_line(lines,nline,"<th>"+txt1+"</th>")
473 set_line(lines,nline,"<th>"+txt2+"</th>")
474 set_line(lines,nline,"<th>"+txt3+"</th>")
475 set_line(lines,nline,"<th>"+txt4+"</th>")
476 set_line(lines,nline,"<th>"+txt5+"</th>")
477 set_line(lines,nline,"<th>"+txt6+"</th>")
479 set_line(lines,nline,row_footer)
481 ;-----------------------------------------------
482 set_line(lines,nline,table_footer)
483 set_line(lines,nline,footer)
485 ; Now write to an HTML file.
486 idx = ind(.not.ismissing(lines))
487 if(.not.any(ismissing(idx))) then
488 asciiwrite(output_html,lines(idx))