Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;********************************************************
2 ; take into account landfrac
3 ; note: landfrac from lnd_T42.nc
4 ; <= lnd_diag_4.0 has correct landfrac
5 ; lnd_diag_3.1 has wrong landfrac
9 ; required command line input parameters:
10 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
12 ;**************************************************************
13 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
14 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
15 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
16 ;**************************************************************
17 procedure set_line(lines:string,nline:integer,newlines:string)
19 ; add line to ascci/html file
21 nnewlines = dimsizes(newlines)
22 if(nline+nnewlines-1.ge.dimsizes(lines))
23 print("set_line: bad index, not setting anything.")
26 lines(nline:nline+nnewlines-1) = newlines
27 ; print ("lines = " + lines(nline:nline+nnewlines-1))
28 nline = nline + nnewlines
31 ;**************************************************************
38 ;************************************************
40 ;************************************************
46 ;model_name_i = "i01.07cn"
47 ;model_name_f = "i01.10cn"
49 model_name_i = "i01.07casa"
50 model_name_f = "i01.10casa"
52 model_name = model_name_f
54 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
55 film_i = model_name_i + "_1990-2004_ANN_climo.nc"
56 film_f = model_name_f + "_1990-2004_ANN_climo.nc"
58 fm_i = addfile (dirm+film_i,"r")
59 fm_f = addfile (dirm+film_f,"r")
70 ;Units for these variables are:
73 nsec_per_year = 60*60*24*365
75 npp_i = npp_i * nsec_per_year
76 npp_f = npp_f * nsec_per_year
78 ;--------------------------------------------------
81 dirm= "/fis/cgd/cseg/people/jeff/surface_data/"
83 fm_l = addfile (dirm+film_l,"r")
85 landfrac = fm_l->landfrac
87 ;npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)
88 ;npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)
90 npp_i = npp_i * conform(npp_i, landfrac, (/1,2/))
91 npp_f = npp_f * conform(npp_f, landfrac, (/1,2/))
93 ;===================================================
94 ; read data: observed at stations
95 ;===================================================
97 station = (/"DukeFACE" \
103 lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
104 lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
105 lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
108 n_sta = dimsizes(station)
109 beta_4_ob = new((/n_sta/),float)
112 ;===================================================
113 ; get model data at station
114 ;===================================================
116 npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
118 npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
123 ;============================
125 ;============================
127 beta_4 = new((/n_sta/),float)
129 beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
131 beta_4_avg = avg(beta_4)
136 ;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
138 bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
139 M_beta = (1. - (bias/n_sta))*3.
143 ;=========================
144 ; for html table - station
145 ;=========================
147 output_html = "table_station.html"
149 ; column (not including header column)
151 col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
153 ncol = dimsizes(col_head)
155 ; row (not including header row)
156 row_head = (/"DukeFACE" \
162 nrow = dimsizes(row_head)
164 ; arrays to be passed to table.
165 text4 = new ((/nrow, ncol/),string )
168 text4(i,0) = sprintf("%.1f",lat_ob(i))
169 text4(i,1) = sprintf("%.1f",lon_ob(i))
170 text4(i,2) = sprintf("%.1f",co2_i)
171 text4(i,3) = sprintf("%.1f",co2_f)
172 text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
173 text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
174 text4(i,6) = sprintf("%.2f",beta_4(i))
177 text4(nrow-1,0) = "-"
178 text4(nrow-1,1) = "-"
179 text4(nrow-1,2) = "-"
180 text4(nrow-1,3) = "-"
181 text4(nrow-1,4) = "-"
182 text4(nrow-1,5) = "-"
183 text4(nrow-1,6) = sprintf("%.2f",beta_4_avg)
184 text4(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
190 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
192 header = (/"<HTML>" \
194 ,"<TITLE>CLAMP metrics</TITLE>" \
201 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
203 ," <th bgcolor=DDDDDD >Station</th>" \
204 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
205 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
206 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
207 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
208 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
209 ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
210 ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
211 ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
214 table_footer = "</table>"
218 lines = new(50000,string)
221 set_line(lines,nline,header)
222 set_line(lines,nline,table_header)
223 ;-----------------------------------------------
227 set_line(lines,nline,row_header)
239 set_line(lines,nline,"<th>"+txt1+"</th>")
240 set_line(lines,nline,"<th>"+txt2+"</th>")
241 set_line(lines,nline,"<th>"+txt3+"</th>")
242 set_line(lines,nline,"<th>"+txt4+"</th>")
243 set_line(lines,nline,"<th>"+txt5+"</th>")
244 set_line(lines,nline,"<th>"+txt6+"</th>")
245 set_line(lines,nline,"<th>"+txt7+"</th>")
246 set_line(lines,nline,"<th>"+txt8+"</th>")
247 set_line(lines,nline,"<th>"+txt9+"</th>")
249 set_line(lines,nline,row_footer)
251 ;-----------------------------------------------
252 set_line(lines,nline,table_footer)
253 set_line(lines,nline,footer)
255 ; Now write to an HTML file.
256 idx = ind(.not.ismissing(lines))
257 if(.not.any(ismissing(idx))) then
258 asciiwrite(output_html,lines(idx))
266 delete (table_header)
269 ;------------------------------------------------
270 ; read biome data: model
272 biome_name_mod = "Model PFT Class"
274 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
275 film = "class_pft_"+model_grid+".nc"
277 fm = addfile(dirm+film,"r")
279 classmod = fm->CLASS_PFT
283 ; model data has 17 land-type classes
287 ;------------------------------------------------
288 ; read biome data: observed
290 biome_name_ob = "MODIS LandCover"
292 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
293 filo = "land_class_"+model_grid+".nc"
295 fo = addfile(diro+filo,"r")
297 classob = tofloat(fo->LAND_CLASS)
301 ; input observed data has 20 land-type classes
305 ;********************************************************************
306 ; use land-type class to bin the data in equally spaced ranges
307 ;********************************************************************
309 ; using observed biome class
311 ; using model biome class
315 range = fspan(0,nclassn-1,nclassn)
318 ; Use this range information to grab all the values in a
319 ; particular range, and then take an average.
323 xvalues = new((/2,nx/),float)
324 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
325 dx = xvalues(0,1) - xvalues(0,0) ; range width
326 dx4 = dx/4 ; 1/4 of the range
327 xvalues(1,:) = xvalues(0,:) - dx/5.
329 ;==============================
331 ;==============================
333 ; using observed biome class
334 ; base_1D = ndtooned(classob)
335 ; using model biome class
336 base_1D = ndtooned(classmod)
338 data1_1D = ndtooned(npp_i)
339 data2_1D = ndtooned(npp_f)
343 yvalues = new((/2,nx/),float)
344 count = new((/2,nx/),float)
348 ; See if we are doing data1 (nd=0) or data2 (nd=1).
358 ; Loop through each range, using base.
361 if (i.ne.(nr-2)) then
363 ; print("In range ["+range(i)+","+range(i+1)+")")
364 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
367 ; print("In range ["+range(i)+",)")
368 idx = ind(base.ge.range(i))
373 if(.not.any(ismissing(idx))) then
374 yvalues(nd,i) = avg(data(idx))
375 count(nd,i) = dimsizes(idx)
377 yvalues(nd,i) = yvalues@_FillValue
381 ;#############################################################
382 ;using observed biome class:
384 ; set the following 4 classes to _FillValue:
385 ; Water Bodies(0), Urban and Build-Up(13),
386 ; Permenant Snow and Ice(15), Unclassified(17)
388 ; if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
389 ; yvalues(nd,i) = yvalues@_FillValue
392 ;#############################################################
394 ;#############################################################
395 ;using model biome class:
397 ; set the following 4 classes to _FillValue:
398 ; (3)Needleleaf Deciduous Boreal Tree,
399 ; (8)Broadleaf Deciduous Boreal Tree,
400 ; (9)Broadleaf Evergreen Shrub,
403 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
404 yvalues(nd,i) = yvalues@_FillValue
407 ;#############################################################
409 ; print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
411 ; Clean up for next time in loop.
419 ;============================
421 ;============================
428 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
432 uu_count = u_count(good)
433 vv_count = v_count(good)
435 n_biome = dimsizes(uu)
436 beta_biome = new((/n_biome/),float)
438 beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
440 ;beta_biome_avg = avg(beta_biome)
441 beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
443 ;print (beta_biome_avg)
445 ;===========================
446 ; for html table - biome
447 ;===========================
449 output_html = "table_biome.html"
451 ; column (not including header column)
453 col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
455 ncol = dimsizes(col_head)
457 ; row (not including header row)
459 ;----------------------------------------------------
460 ; using observed biome class:
462 ; row_head = (/"Evergreen Needleleaf Forests" \
463 ; ,"Evergreen Broadleaf Forests" \
464 ; ,"Deciduous Needleleaf Forest" \
465 ; ,"Deciduous Broadleaf Forests" \
467 ; ,"Closed Bushlands" \
468 ; ,"Open Bushlands" \
469 ; ,"Woody Savannas (S. Hem.)" \
470 ; ,"Savannas (S. Hem.)" \
472 ; ,"Permanent Wetlands" \
474 ; ,"Cropland/Natural Vegetation Mosaic" \
475 ; ,"Barren or Sparsely Vegetated" \
476 ; ,"Woody Savannas (N. Hem.)" \
477 ; ,"Savannas (N. Hem.)" \
481 ;----------------------------------------------------
482 ; using model biome class:
484 row_head = (/"Not Vegetated" \
485 ,"Needleleaf Evergreen Temperate Tree" \
486 ,"Needleleaf Evergreen Boreal Tree" \
487 ; ,"Needleleaf Deciduous Boreal Tree" \
488 ,"Broadleaf Evergreen Tropical Tree" \
489 ,"Broadleaf Evergreen Temperate Tree" \
490 ,"Broadleaf Deciduous Tropical Tree" \
491 ,"Broadleaf Deciduous Temperate Tree" \
492 ; ,"Broadleaf Deciduous Boreal Tree" \
493 ; ,"Broadleaf Evergreen Shrub" \
494 ,"Broadleaf Deciduous Temperate Shrub" \
495 ,"Broadleaf Deciduous Boreal Shrub" \
497 ,"C3 Non-Arctic Grass" \
504 nrow = dimsizes(row_head)
506 ; arrays to be passed to table.
507 text4 = new ((/nrow, ncol/),string )
510 text4(i,0) = sprintf("%.1f",co2_i)
511 text4(i,1) = sprintf("%.1f",co2_f)
512 text4(i,2) = sprintf("%.1f",uu(i))
513 text4(i,3) = sprintf("%.1f",vv(i))
514 text4(i,4) = sprintf("%.2f",beta_biome(i))
516 text4(nrow-1,0) = "-"
517 text4(nrow-1,1) = "-"
518 text4(nrow-1,2) = "-"
519 text4(nrow-1,3) = "-"
520 text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
522 ;**************************************************
524 ;**************************************************
526 header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
528 header = (/"<HTML>" \
530 ,"<TITLE>CLAMP metrics</TITLE>" \
537 "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
539 ," <th bgcolor=DDDDDD >Biome Class</th>" \
540 ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
541 ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
542 ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
543 ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
544 ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
547 table_footer = "</table>"
551 lines = new(50000,string)
554 set_line(lines,nline,header)
555 set_line(lines,nline,table_header)
556 ;-----------------------------------------------
560 set_line(lines,nline,row_header)
569 set_line(lines,nline,"<th>"+txt1+"</th>")
570 set_line(lines,nline,"<th>"+txt2+"</th>")
571 set_line(lines,nline,"<th>"+txt3+"</th>")
572 set_line(lines,nline,"<th>"+txt4+"</th>")
573 set_line(lines,nline,"<th>"+txt5+"</th>")
574 set_line(lines,nline,"<th>"+txt6+"</th>")
576 set_line(lines,nline,row_footer)
578 ;-----------------------------------------------
579 set_line(lines,nline,table_footer)
580 set_line(lines,nline,footer)
582 ; Now write to an HTML file.
583 idx = ind(.not.ismissing(lines))
584 if(.not.any(ismissing(idx))) then
585 asciiwrite(output_html,lines(idx))