diff -r 000000000000 -r 0c6405ab2ff4 lai/99.all.ncl.new --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lai/99.all.ncl.new Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,735 @@ +;************************************************************* +; remove histogram plots. +; +; required command line input parameters: +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl +; +; histogram normalized by rain and compute correleration +;************************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" +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 +;************************************************ + +;film = "b30.061n_1995-2004_MONS_climo_lnd.nc" +;model_name = "b30.061n" +;model_grid = "T31" + +;film = "newcn05_ncep_1i_MONS_climo_lnd.nc" +;model_name = "newcn" +;model_grid = "1.9" + +;film = "i01.06cn_1798-2004_MONS_climo.nc" +;model_name = "06cn" +;model_grid = "T42" + +;film = "i01.06casa_1798-2004_MONS_climo.nc" +;model_name = "06casa" +;model_grid = "T42" + + film = "i01.10cn_1948-2004_MONS_climo.nc" + model_name = "10cn" + model_grid = "T42" + +;film = "i01.10casa_1948-2004_MONS_climo.nc" +;model_name = "10casa" +;model_grid = "T42" + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + fm = addfile(dirm+film,"r") + + laimod = fm->TLAI + +;************************************************ +; read data: observed +;************************************************ + + ob_name = "MODIS MOD 15A2 2000-2005" + + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" + filo1 = "land_class_"+model_grid+".nc" + filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc" + + fo1 = addfile(diro+filo1,"r") + fo2 = addfile(diro+filo2,"r") + + classob = tofloat(fo1->LAND_CLASS) + laiob = fo2->LAI + +; input observed data has 20 land-type classes + nclass = 20 + +;************************************************ +; global res +;************************************************ + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + +;************************************************ +; plot global land class: observed +;************************************************ + + resg@cnMinLevelValF = 1. ; Min level + resg@cnMaxLevelValF = 19. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +; global contour ob + classob@_FillValue = 1.e+36 + classob = where(classob.eq.0,classob@_FillValue,classob) + + plot_name = "global_class_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,classob,resg) + frame(wks) + + clear (wks) + +;******************************************************************* +; for html table : all 4 components (Mean, Max, Phase, Growth) +;******************************************************************* + +; column (not including header column) + + component = (/"Mean","Max","Phase","Growth"/) + col_head2 = (/"observed",model_name,"M_score" \ + ,"observed",model_name,"M_score" \ + ,"observed",model_name,"M_score" \ + ,"observed",model_name,"M_score" \ + /) + + n_comp = dimsizes(component) + ncol = dimsizes(col_head2) + +; row (not including header row) + 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" \ + /) + nrow = dimsizes(row_head) + +; arrays to be passed to table. + text4 = new ((/nrow, ncol/),string ) + +; M_comp + M_comp = (/"M_lai_mean","M_lai_max","M_lai_phase","M_lai_grow"/) + +; total M_score + M_total = 0. + +;******************************************************************** +; use land-type class to bin the data in equally spaced ranges +;******************************************************************** + + 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. + +;************************************************************************ +; go through all components +;************************************************************************ + + do n = 0,n_comp-1 + +;=================== +; get data: +;=================== +; (A) Mean + + if (n .eq. 0) then + data_ob = dim_avg_Wrap(laiob (lat|:,lon|:,time|:)) + data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:)) + end if + +; (B) Max + + if (n .eq. 1) then + +; observed + data_ob = laiob(0,:,:) + s = laiob(:,0,0) + data_ob@long_name = "Leaf Area Index Max" + + dsizes_z = dimsizes(laiob) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laiob(:,j,i) + data_ob(j,i) = max(s) + end do + end do + + delete (s) + delete (dsizes_z) + +; model + data_mod = laimod(0,:,:) + s = laimod(:,0,0) + data_mod@long_name = "Leaf Area Index Max" + + dsizes_z = dimsizes(laimod) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laimod(:,j,i) + data_mod(j,i) = max(s) + end do + end do + + delete (s) + delete (dsizes_z) + end if + +; (C) phase + + if (n .eq. 2) then + +; observed + data_ob = laiob(0,:,:) + s = laiob(:,0,0) + data_ob@long_name = "Leaf Area Index Max Month" + + dsizes_z = dimsizes(laiob) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laiob(:,j,i) + data_ob(j,i) = maxind(s) + 1 + end do + end do + + delete (s) + delete (dsizes_z) + +; model + data_mod = laimod(0,:,:) + s = laimod(:,0,0) + data_mod@long_name = "Leaf Area Index Max Month" + + dsizes_z = dimsizes(laimod) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laimod(:,j,i) + data_mod(j,i) = maxind(s) + 1 + end do + end do + + delete (s) + delete (dsizes_z) + end if + +; (D) grow day + + if (n .eq. 3) then + + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) + +; observed + data_ob = laiob(0,:,:) + data_ob@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(laiob) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + + data_ob(j,i) = nday + end do + end do + + delete (dsizes_z) + +; model + data_mod = laimod(0,:,:) + data_mod@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(laimod) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + + data_mod(j,i) = nday + end do + end do + + delete (dsizes_z) + end if + +;============================== +; put data into bins +;============================== + + base_1D = ndtooned(classob) + data1_1D = ndtooned(data_ob) + data2_1D = ndtooned(data_mod) + +; output for data in bins + + yvalues = new((/2,nx/),float) + count = new((/2,nx/),float) + +; put data into bins + + 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 + +;############################################################# +; 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 +;############################################################# + +; 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 + + delete (base) + delete (base_1D) + delete (data1_1D) + delete (data2_1D) + +;===================================== +; compute correlation coef and M score +;===================================== + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + +; compute correlation coef + cc = esccr(uu,vv,0) + + if (n .eq. 2) then + bias = avg(abs(vv-uu)) + bias = where((bias.gt. 6.),12.-bias,bias) + Mscore = ((6. - bias)/6.)*5. + M_score = sprintf("%.2f", Mscore) + else + bias = sum(abs(vv-uu)/abs(vv+uu)) + Mscore = (1.- (bias/dimsizes(uu)))*5. + M_score = sprintf("%.2f", Mscore) + end if + +; compute M_total + + M_total = M_total + Mscore + +;================== +; output M_score +;================== + + print (Mscore) +;======================= +; output to html table +;======================= + + nn = n*3 + + do i=0,nrow-2 + text4(i,nn) = sprintf("%.2f",u(i)) + text4(i,nn+1) = sprintf("%.2f",v(i)) + text4(i,nn+2) = "-" + end do + text4(nrow-1,nn) = sprintf("%.2f",avg(u)) + text4(nrow-1,nn+1) = sprintf("%.2f",avg(v)) + text4(nrow-1,nn+2) = M_score + + delete (u) + delete (v) + delete (uu) + delete (vv) + delete (yvalues) + delete (good) + +;======================================== +; global res changes for each component +;======================================== + delta = 0.00001 + + if (n .eq. 0) then + resg@cnMinLevelValF = 0. + resg@cnMaxLevelValF = 10. + resg@cnLevelSpacingF = 1. + + data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) + end if + + if (n .eq. 1) then + resg@cnMinLevelValF = 0. + resg@cnMaxLevelValF = 10. + resg@cnLevelSpacingF = 1. + + data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) + end if + + if (n .eq. 2) then + resg@cnMinLevelValF = 1. + resg@cnMaxLevelValF = 12. + resg@cnLevelSpacingF = 1. + + data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) + end if + + if (n .eq. 3) then + resg@cnMinLevelValF = 60. + resg@cnMaxLevelValF = 360. + resg@cnLevelSpacingF = 20. + + data_ob@_FillValue = 1.e+36 + data_ob = where(data_ob .lt. 10.,data_ob@_FillValue,data_ob) + + data_mod@_FillValue = 1.e+36 + data_mod = where(data_mod .lt. 10.,data_mod@_FillValue,data_mod) + end if + +;========================= +; global contour : ob +;========================= + + plot_name = "global_"+component(n)+"_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,data_ob,resg) + frame(wks) + + clear (wks) + delete (plot) + +;============================ +; global contour : model +;============================ + + plot_name = "global_"+component(n)+"_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") + + plot = gsn_csm_contour_map_ce(wks,data_mod,resg) + frame(wks) + + clear (wks) + delete (plot) + +;================================ +; global contour: model vs ob +;================================ + + plot_name = "global_"+component(n)+"_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") + + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) + +; plot ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg) + +; plot model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg) + +; plot model-ob + + if (n .eq. 0) then + resg@cnMinLevelValF = -2. + resg@cnMaxLevelValF = 2. + resg@cnLevelSpacingF = 0.4 + end if + + if (n .eq. 1) then + resg@cnMinLevelValF = -6. + resg@cnMaxLevelValF = 6. + resg@cnLevelSpacingF = 1. + end if + + if (n .eq. 2) then + resg@cnMinLevelValF = -6. + resg@cnMaxLevelValF = 6. + resg@cnLevelSpacingF = 1. + end if + + if (n .eq. 3) then + resg@cnMinLevelValF = -100. + resg@cnMaxLevelValF = 100. + resg@cnLevelSpacingF = 20. + end if + + zz = data_mod + zz = data_mod - data_ob + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + +; plot panel + + pres = True ; panel plot mods desired + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + clear (wks) + delete (plot) + + end do +;************************************************** +; html table +;************************************************** + output_html = "table_model_vs_ob.html" + + header_text = "

LAI: Model "+model_name+" vs Observed

" + + header = (/"" \ + ,"" \ + ,"CLAMP metrics" \ + ,"" \ + ,header_text \ + /) + footer = "" + + table_header = (/ \ + "" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + ,"" \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ," " \ + ,"" \ + /) + table_footer = "
Biome Class"+component(0)+""+component(1)+""+component(2)+""+component(3)+"
observed"+model_name+"M_scoreobserved"+model_name+"M_scoreobserved"+model_name+"M_scoreobserved"+model_name+"M_score
" + 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) + txt10 = text4(n,8) + txt11 = text4(n,9) + txt12 = text4(n,10) + txt13 = text4(n,11) + + 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,""+txt10+"") + set_line(lines,nline,""+txt11+"") + set_line(lines,nline,""+txt12+"") + set_line(lines,nline,""+txt13+"") + + 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 + +;*************************************************************************** +; write total score to file +;*************************************************************************** + + asciiwrite("M_save.lai", M_total) + +;*************************************************************************** +end +