diff -r 000000000000 -r 0c6405ab2ff4 beta/07.landfrac.ncl
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/beta/07.landfrac.ncl Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,591 @@
+;********************************************************
+; take into account landfrac
+; note: landfrac from lnd_T42.nc
+; <= lnd_diag_4.0 has correct landfrac
+; lnd_diag_3.1 has wrong landfrac
+;
+; using model biome
+;
+; required command line input parameters:
+; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
+;
+;**************************************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+;**************************************************************
+procedure set_line(lines:string,nline:integer,newlines:string)
+begin
+; add line to ascci/html file
+
+ nnewlines = dimsizes(newlines)
+ if(nline+nnewlines-1.ge.dimsizes(lines))
+ print("set_line: bad index, not setting anything.")
+ return
+ end if
+ lines(nline:nline+nnewlines-1) = newlines
+; print ("lines = " + lines(nline:nline+nnewlines-1))
+ nline = nline + nnewlines
+ return
+end
+;**************************************************************
+; Main code.
+begin
+
+ plot_type = "ps"
+ plot_type_new = "png"
+
+;************************************************
+; read data: model
+;************************************************
+ co2_i = 283.1878
+ co2_f = 364.1252
+
+ model_grid = "T42"
+
+;model_name_i = "i01.07cn"
+;model_name_f = "i01.10cn"
+
+ model_name_i = "i01.07casa"
+ model_name_f = "i01.10casa"
+
+ model_name = model_name_f
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film_i = model_name_i + "_1990-2004_ANN_climo.nc"
+ film_f = model_name_f + "_1990-2004_ANN_climo.nc"
+
+ fm_i = addfile (dirm+film_i,"r")
+ fm_f = addfile (dirm+film_f,"r")
+
+ xm = fm_f->lon
+ ym = fm_f->lat
+
+ npp_i = fm_i->NPP
+ npp_f = fm_f->NPP
+
+ delete (fm_i)
+ delete (fm_f)
+
+;Units for these variables are:
+;npp_i: g C/m^2/s
+
+ nsec_per_year = 60*60*24*365
+
+ npp_i = npp_i * nsec_per_year
+ npp_f = npp_f * nsec_per_year
+
+;--------------------------------------------------
+;get landfrac data
+
+ dirm= "/fis/cgd/cseg/people/jeff/surface_data/"
+ film_l = "lnd_T42.nc"
+ fm_l = addfile (dirm+film_l,"r")
+
+ landfrac = fm_l->landfrac
+
+;npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)
+;npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)
+
+ npp_i = npp_i * conform(npp_i, landfrac, (/1,2/))
+ npp_f = npp_f * conform(npp_f, landfrac, (/1,2/))
+
+;===================================================
+; read data: observed at stations
+;===================================================
+
+ station = (/"DukeFACE" \
+ ,"AspenFACE" \
+ ,"ORNL-FACE" \
+ ,"POP-EUROFACE" \
+ /)
+
+ lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
+ lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
+ lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
+;print (lon_ob)
+
+ n_sta = dimsizes(station)
+ beta_4_ob = new((/n_sta/),float)
+ beta_4_ob = 0.60
+
+;===================================================
+; get model data at station
+;===================================================
+
+ npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
+
+ npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
+
+;print (npp_i_4)
+;print (npp_f_4)
+
+;============================
+;compute beta_4
+;============================
+
+ beta_4 = new((/n_sta/),float)
+
+ beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
+
+ beta_4_avg = avg(beta_4)
+
+;print (beta_4)
+;print (beta_4_avg)
+
+;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
+
+ bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
+ M_beta = (1. - (bias/n_sta))*3.
+
+ print (M_beta)
+
+;=========================
+; for html table - station
+;=========================
+
+ output_html = "table_station.html"
+
+; column (not including header column)
+
+ col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
+
+ ncol = dimsizes(col_head)
+
+; row (not including header row)
+ row_head = (/"DukeFACE" \
+ ,"AspenFACE" \
+ ,"ORNL-FACE" \
+ ,"POP-EUROFACE" \
+ ,"All Station" \
+ /)
+ nrow = dimsizes(row_head)
+
+; arrays to be passed to table.
+ text4 = new ((/nrow, ncol/),string )
+
+ do i=0,nrow-2
+ text4(i,0) = sprintf("%.1f",lat_ob(i))
+ text4(i,1) = sprintf("%.1f",lon_ob(i))
+ text4(i,2) = sprintf("%.1f",co2_i)
+ text4(i,3) = sprintf("%.1f",co2_f)
+ text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
+ text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
+ text4(i,6) = sprintf("%.2f",beta_4(i))
+ text4(i,7) = "-"
+ end do
+ text4(nrow-1,0) = "-"
+ text4(nrow-1,1) = "-"
+ text4(nrow-1,2) = "-"
+ text4(nrow-1,3) = "-"
+ text4(nrow-1,4) = "-"
+ text4(nrow-1,5) = "-"
+ text4(nrow-1,6) = sprintf("%.2f",beta_4_avg)
+ text4(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
+
+;-----------
+; html table
+;-----------
+
+ header_text = "
Beta Factor: Model "+model_name+"
"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,header_text \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Station | " \
+ ," "+col_head(0)+" | " \
+ ," "+col_head(1)+" | " \
+ ," "+col_head(2)+" | " \
+ ," "+col_head(3)+" | " \
+ ," "+col_head(4)+" | " \
+ ," "+col_head(5)+" | " \
+ ," "+col_head(6)+" | " \
+ ," "+col_head(7)+" | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+ do n = 0,nrow-1
+ set_line(lines,nline,row_header)
+
+ txt1 = row_head(n)
+ txt2 = text4(n,0)
+ txt3 = text4(n,1)
+ txt4 = text4(n,2)
+ txt5 = text4(n,3)
+ txt6 = text4(n,4)
+ txt7 = text4(n,5)
+ txt8 = text4(n,6)
+ txt9 = text4(n,7)
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+ set_line(lines,nline,""+txt7+" | ")
+ set_line(lines,nline,""+txt8+" | ")
+ set_line(lines,nline,""+txt9+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+
+ delete (col_head)
+ delete (row_head)
+ delete (text4)
+ delete (table_header)
+ delete (idx)
+
+;------------------------------------------------
+; read biome data: model
+
+ biome_name_mod = "Model PFT Class"
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = "class_pft_"+model_grid+".nc"
+
+ fm = addfile(dirm+film,"r")
+
+ classmod = fm->CLASS_PFT
+
+ delete (fm)
+
+; model data has 17 land-type classes
+
+ nclass_mod = 17
+
+;------------------------------------------------
+; read biome data: observed
+
+ biome_name_ob = "MODIS LandCover"
+
+ diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
+ filo = "land_class_"+model_grid+".nc"
+
+ fo = addfile(diro+filo,"r")
+
+ classob = tofloat(fo->LAND_CLASS)
+
+ delete (fo)
+
+; input observed data has 20 land-type classes
+
+ nclass_ob = 20
+
+;********************************************************************
+; use land-type class to bin the data in equally spaced ranges
+;********************************************************************
+
+; using observed biome class
+; nclass = nclass_ob
+; using model biome class
+ nclass = nclass_mod
+
+ nclassn = nclass + 1
+ range = fspan(0,nclassn-1,nclassn)
+; print (range)
+
+; Use this range information to grab all the values in a
+; particular range, and then take an average.
+
+ nr = dimsizes(range)
+ nx = nr-1
+ xvalues = new((/2,nx/),float)
+ xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
+ dx = xvalues(0,1) - xvalues(0,0) ; range width
+ dx4 = dx/4 ; 1/4 of the range
+ xvalues(1,:) = xvalues(0,:) - dx/5.
+
+;==============================
+; put data into bins
+;==============================
+
+; using observed biome class
+; base_1D = ndtooned(classob)
+; using model biome class
+ base_1D = ndtooned(classmod)
+
+ data1_1D = ndtooned(npp_i)
+ data2_1D = ndtooned(npp_f)
+
+; output
+
+ yvalues = new((/2,nx/),float)
+ count = new((/2,nx/),float)
+
+ do nd=0,1
+
+; See if we are doing data1 (nd=0) or data2 (nd=1).
+
+ base = base_1D
+
+ if(nd.eq.0) then
+ data = data1_1D
+ else
+ data = data2_1D
+ end if
+
+; Loop through each range, using base.
+
+ do i=0,nr-2
+ if (i.ne.(nr-2)) then
+; print("")
+; print("In range ["+range(i)+","+range(i+1)+")")
+ idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
+ else
+; print("")
+; print("In range ["+range(i)+",)")
+ idx = ind(base.ge.range(i))
+ end if
+
+; Calculate average
+
+ if(.not.any(ismissing(idx))) then
+ yvalues(nd,i) = avg(data(idx))
+ count(nd,i) = dimsizes(idx)
+ else
+ yvalues(nd,i) = yvalues@_FillValue
+ count(nd,i) = 0
+ end if
+
+;#############################################################
+;using observed biome class:
+;
+; set the following 4 classes to _FillValue:
+; Water Bodies(0), Urban and Build-Up(13),
+; Permenant Snow and Ice(15), Unclassified(17)
+
+; if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
+; yvalues(nd,i) = yvalues@_FillValue
+; count(nd,i) = 0
+; end if
+;#############################################################
+
+;#############################################################
+;using model biome class:
+;
+; set the following 4 classes to _FillValue:
+; (3)Needleleaf Deciduous Boreal Tree,
+; (8)Broadleaf Deciduous Boreal Tree,
+; (9)Broadleaf Evergreen Shrub,
+; (16)Wheat
+
+ if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
+ yvalues(nd,i) = yvalues@_FillValue
+ count(nd,i) = 0
+ end if
+;#############################################################
+
+; print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
+
+; Clean up for next time in loop.
+
+ delete(idx)
+ end do
+
+ delete(data)
+ end do
+
+;============================
+;compute beta
+;============================
+
+ u = yvalues(0,:)
+ v = yvalues(1,:)
+ u_count = count(0,:)
+ v_count = count(1,:)
+
+ good = ind(.not.ismissing(u) .and. .not.ismissing(v))
+
+ uu = u(good)
+ vv = v(good)
+ uu_count = u_count(good)
+ vv_count = v_count(good)
+
+ n_biome = dimsizes(uu)
+ beta_biome = new((/n_biome/),float)
+
+ beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
+
+;beta_biome_avg = avg(beta_biome)
+ beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
+
+;print (beta_biome_avg)
+
+;===========================
+; for html table - biome
+;===========================
+
+ output_html = "table_biome.html"
+
+; column (not including header column)
+
+ col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
+
+ ncol = dimsizes(col_head)
+
+; row (not including header row)
+
+;----------------------------------------------------
+; using observed biome class:
+;
+; row_head = (/"Evergreen Needleleaf Forests" \
+; ,"Evergreen Broadleaf Forests" \
+; ,"Deciduous Needleleaf Forest" \
+; ,"Deciduous Broadleaf Forests" \
+; ,"Mixed Forests" \
+; ,"Closed Bushlands" \
+; ,"Open Bushlands" \
+; ,"Woody Savannas (S. Hem.)" \
+; ,"Savannas (S. Hem.)" \
+; ,"Grasslands" \
+; ,"Permanent Wetlands" \
+; ,"Croplands" \
+; ,"Cropland/Natural Vegetation Mosaic" \
+; ,"Barren or Sparsely Vegetated" \
+; ,"Woody Savannas (N. Hem.)" \
+; ,"Savannas (N. Hem.)" \
+; ,"All Biome" \
+; /)
+
+;----------------------------------------------------
+; using model biome class:
+;
+ row_head = (/"Not Vegetated" \
+ ,"Needleleaf Evergreen Temperate Tree" \
+ ,"Needleleaf Evergreen Boreal Tree" \
+; ,"Needleleaf Deciduous Boreal Tree" \
+ ,"Broadleaf Evergreen Tropical Tree" \
+ ,"Broadleaf Evergreen Temperate Tree" \
+ ,"Broadleaf Deciduous Tropical Tree" \
+ ,"Broadleaf Deciduous Temperate Tree" \
+; ,"Broadleaf Deciduous Boreal Tree" \
+; ,"Broadleaf Evergreen Shrub" \
+ ,"Broadleaf Deciduous Temperate Shrub" \
+ ,"Broadleaf Deciduous Boreal Shrub" \
+ ,"C3 Arctic Grass" \
+ ,"C3 Non-Arctic Grass" \
+ ,"C4 Grass" \
+ ,"Corn" \
+; ,"Wheat" \
+ ,"All Biome" \
+ /)
+
+ nrow = dimsizes(row_head)
+
+; arrays to be passed to table.
+ text4 = new ((/nrow, ncol/),string )
+
+ do i=0,nrow-2
+ text4(i,0) = sprintf("%.1f",co2_i)
+ text4(i,1) = sprintf("%.1f",co2_f)
+ text4(i,2) = sprintf("%.1f",uu(i))
+ text4(i,3) = sprintf("%.1f",vv(i))
+ text4(i,4) = sprintf("%.2f",beta_biome(i))
+ end do
+ text4(nrow-1,0) = "-"
+ text4(nrow-1,1) = "-"
+ text4(nrow-1,2) = "-"
+ text4(nrow-1,3) = "-"
+ text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
+
+;**************************************************
+; html table
+;**************************************************
+
+ header_text = "Beta Factor: Model "+model_name+"
"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,header_text \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Biome Class | " \
+ ," "+col_head(0)+" | " \
+ ," "+col_head(1)+" | " \
+ ," "+col_head(2)+" | " \
+ ," "+col_head(3)+" | " \
+ ," "+col_head(4)+" | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+ do n = 0,nrow-1
+ set_line(lines,nline,row_header)
+
+ txt1 = row_head(n)
+ txt2 = text4(n,0)
+ txt3 = text4(n,1)
+ txt4 = text4(n,2)
+ txt5 = text4(n,3)
+ txt6 = text4(n,4)
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+
+end
+