diff -r 000000000000 -r 0c6405ab2ff4 beta/06.biome_model.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/beta/06.biome_model.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,571 @@ +;******************************************************** +; 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 + +;=================================================== +; 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 = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
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)+"
" + 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 = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Biome Class"+col_head(0)+""+col_head(1)+""+col_head(2)+""+col_head(3)+""+col_head(4)+"
" + 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 +