forrest@0: ;*****************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;****************************************************************************
forrest@0: procedure set_line(lines:string,nline:integer,newlines:string)
forrest@0: begin
forrest@0: ; add line to ascci/html file
forrest@0:
forrest@0: nnewlines = dimsizes(newlines)
forrest@0: if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0: print("set_line: bad index, not setting anything.")
forrest@0: return
forrest@0: end if
forrest@0: lines(nline:nline+nnewlines-1) = newlines
forrest@0: ; print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0: nline = nline + nnewlines
forrest@0: return
forrest@0: end
forrest@0: ;****************************************************************************
forrest@0: procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
forrest@0: ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
forrest@0: ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
forrest@0: ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
forrest@0: ,dx4[1]:numeric \
forrest@0: )
forrest@0: begin
forrest@0: ; Calculaee "nice" bins for binning the data in equally spaced ranges.
forrest@0: ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
forrest@0: ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
forrest@0:
forrest@0: nbins = 15 ; Number of bins to use.
forrest@0:
forrest@0: nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
forrest@0: nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
forrest@0: range = fspan(nicevals(0),nicevals(1),nvals)
forrest@0:
forrest@0: ; Use this range information to grab all the values in a
forrest@0: ; particular range, and then take an average.
forrest@0:
forrest@0: nr = dimsizes(range)
forrest@0: nx = nr-1
forrest@0: ; print (nx)
forrest@0:
forrest@0: ; xvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0: xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
forrest@0: dx = xvalues(0,1) - xvalues(0,0) ; range width
forrest@0: dx4 = dx/4 ; 1/4 of the range
forrest@0: xvalues(1,:) = xvalues(0,:) - dx/5.
forrest@0: ; yvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0: ; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0: ; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
forrest@0:
forrest@0: do nd=0,1
forrest@0:
forrest@0: ; See if we are doing model or observational data.
forrest@0:
forrest@0: if(nd.eq.0) then
forrest@0: data = RAIN1_1D
forrest@0: npp_data = NPP1_1D
forrest@0: else
forrest@0: data = RAIN2_1D
forrest@0: npp_data = NPP2_1D
forrest@0: end if
forrest@0:
forrest@0: ; Loop through each range and check for values.
forrest@0:
forrest@0: do i=0,nr-2
forrest@0: if (i.ne.(nr-2)) then
forrest@0: ; print("")
forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0: idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
forrest@0: else
forrest@0: ; print("")
forrest@0: ; print("In range ["+range(i)+",)")
forrest@0: idx = ind(range(i).le.data)
forrest@0: end if
forrest@0:
forrest@0: ; Calculate average, and get min and max.
forrest@0:
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: yvalues(nd,i) = avg(npp_data(idx))
forrest@0: mn_yvalues(nd,i) = min(npp_data(idx))
forrest@0: mx_yvalues(nd,i) = max(npp_data(idx))
forrest@0: count = dimsizes(idx)
forrest@0: else
forrest@0: count = 0
forrest@0: yvalues(nd,i) = yvalues@_FillValue
forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0: end if
forrest@0:
forrest@0: ; Print out information.
forrest@0: ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0:
forrest@0: ; Clean up for next time in loop.
forrest@0:
forrest@0: delete(idx)
forrest@0: end do
forrest@0: delete(data)
forrest@0: delete(npp_data)
forrest@0: end do
forrest@0: end
forrest@0: ;****************************************************************************
forrest@0: ; Main code.
forrest@0: ;****************************************************************************
forrest@0: begin
forrest@0:
forrest@0: plot_type = "ps"
forrest@0: plot_type_new = "png"
forrest@0:
forrest@0: ;-----------------------------------------------------
forrest@0: ; edit table.html of current model for movel1_vs_model2
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: html_name2 = compare+"/table.html"
forrest@0: html_new2 = html_name2 +".new"
forrest@0: end if
forrest@0: ;------------------------------------------------------
forrest@0: ; edit table.html for current model
forrest@0:
forrest@0: html_name = model_name+"/table.html"
forrest@0: html_new = html_name +".new"
forrest@0:
forrest@0: ;------------------------------------------------------
forrest@0: ; read model data
forrest@0:
forrest@0: fm = addfile (dirm+film1,"r")
forrest@0:
forrest@0: nppmod0 = fm->NPP
forrest@0: rainmod0 = fm->RAIN
forrest@0: xm = fm->lon
forrest@0: ym = fm->lat
forrest@0:
forrest@0: delete (fm)
forrest@0:
forrest@0: ;----------------------------------------------------
forrest@0: ; read ob data
forrest@0:
forrest@0: ;(1) data at 81 sites
forrest@0:
forrest@0: dir81 = diro + "npp/"
forrest@0: fil81 = "data.81.nc"
forrest@0: f81 = addfile (dir81+fil81,"r")
forrest@0:
forrest@0: id81 = f81->SITE_ID
forrest@0: npp81 = f81->TNPP_C
forrest@0: rain81 = tofloat(f81->PREC_ANN)
forrest@0: x81 = f81->LONG_DD
forrest@0: y81 = f81->LAT_DD
forrest@0:
forrest@0: delete (f81)
forrest@0:
forrest@0: id81@long_name = "SITE_ID"
forrest@0:
forrest@0: ;change longitude from (-180,180) to (0,360)
forrest@0: ;for model data interpolation
forrest@0:
forrest@0: do i= 0,dimsizes(x81)-1
forrest@0: if (x81(i) .lt. 0.) then
forrest@0: x81(i) = x81(i)+ 360.
forrest@0: end if
forrest@0: end do
forrest@0:
forrest@0: ;-------------------------------------
forrest@0: ;(2) data at 933 sites
forrest@0:
forrest@0: dir933 = diro + "npp/"
forrest@0: fil933 = "data.933.nc"
forrest@0: f933 = addfile (dir933+fil933,"r")
forrest@0:
forrest@0: id933 = f933->SITE_ID
forrest@0: npp933 = f933->TNPP_C
forrest@0: rain933 = f933->PREC
forrest@0: x933 = f933->LONG_DD
forrest@0: y933 = f933->LAT_DD
forrest@0:
forrest@0: delete (f933)
forrest@0:
forrest@0: id933@long_name = "SITE_ID"
forrest@0:
forrest@0: ;change longitude from (-180,180) to (0,360)
forrest@0: ;for model data interpolation
forrest@0:
forrest@0: do i= 0,dimsizes(x933)-1
forrest@0: if (x933(i) .lt. 0.) then
forrest@0: x933(i) = x933(i)+ 360.
forrest@0: end if
forrest@0: end do
forrest@0:
forrest@0: ;----------------------------------------
forrest@0: ;(3) global data, interpolated into model grid
forrest@0:
forrest@0: dirg = diro + "npp/"
forrest@1: ;filg = "Npp_"+model_grid+"_mean.nc"
forrest@1: filg = "npp_"+model_grid+"_mean_2000-2004.nc"
forrest@0: fglobe = addfile (dirg+filg,"r")
forrest@0:
forrest@0: nppglobe0 = fglobe->NPP
forrest@0: nppglobe = nppglobe0
forrest@0:
forrest@0: delete (fglobe)
forrest@0:
forrest@0: ;***********************************************************************
forrest@0: ; interpolate model data into ob sites
forrest@0: ;***********************************************************************
forrest@0:
forrest@0: nppmod = nppmod0(0,:,:)
forrest@0: rainmod = rainmod0(0,:,:)
forrest@0: delete (nppmod0)
forrest@0: delete (rainmod0)
forrest@0:
forrest@0: nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
forrest@0:
forrest@0: nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
forrest@0:
forrest@0: rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
forrest@0:
forrest@0: rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
forrest@0:
forrest@0: ;**********************************************************
forrest@0: ; unified units
forrest@0: ;**********************************************************
forrest@0: ; Units for these variables are:
forrest@0: ;
forrest@0: ; rain81 : mm/year
forrest@0: ; rainmod : mm/s
forrest@1: ; npp81 : gC/m^2/year
forrest@1: ; nppmod81: gC/m^2/s
forrest@1: ; nppglobe: gC/m^2/year
forrest@0: ;
forrest@1: ; We want to convert these to "m/year" and "gC/m^2/year".
forrest@0:
forrest@0: nsec_per_year = 60*60*24*365
forrest@0:
forrest@0: rain81 = rain81 / 1000.
forrest@0: rainmod81 = (rainmod81/ 1000.) * nsec_per_year
forrest@0: nppmod81 = nppmod81 * nsec_per_year
forrest@0:
forrest@0: rain933 = rain933 / 1000.
forrest@0: rainmod933 = (rainmod933/ 1000.) * nsec_per_year
forrest@0: nppmod933 = nppmod933 * nsec_per_year
forrest@0:
forrest@0: nppmod = nppmod * nsec_per_year
forrest@0:
forrest@1: npp81@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: nppmod81@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: npp933@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: nppmod933@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: nppmod@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: nppglobe@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: rain81@units = "m y~S~-1~N~"
forrest@1: rainmod81@units = "m y~S~-1~N~"
forrest@1: rain933@units = "m y~S~-1~N~"
forrest@1: rainmod933@units = "m y~S~-1~N~"
forrest@0:
forrest@1: npp81@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: npp933@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: nppmod81@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: nppmod933@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: nppglobe@long_name = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: rain81@long_name = "PREC (m y~S~-1~N~)"
forrest@1: rain933@long_name = "PREC (m y~S~-1~N~)"
forrest@1: rainmod81@long_name = "PREC (m y~S~-1~N~)"
forrest@1: rainmod933@long_name = "PREC (m y~S~-1~N~)"
forrest@0:
forrest@0: ; change longitude from 0-360 to -180-180
forrest@0: x81 = where(x81 .gt. 180., x81 -360., x81 )
forrest@0: x933 = where(x933 .gt. 180., x933-360., x933)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ;(A)-1 html table of site81 -- observed
forrest@0: ;*******************************************************************
forrest@0:
forrest@0: output_html = "table_site81_ob.html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"EMDI Observations Class A (81 sites)
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," | Site ID | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@1: ," NPP (gC m-2 y-1) | " \
forrest@1: ," PPT (m y-1) | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: nrow = dimsizes(id81)
forrest@0: do n = 0,nrow-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt1 = sprintf("%5.0f", id81(n))
forrest@0: txt2 = sprintf("%5.2f", y81(n))
forrest@0: txt3 = sprintf("%5.2f", x81(n))
forrest@0: txt4 = sprintf("%5.2f", npp81(n))
forrest@0: txt5 = sprintf("%5.2f", rain81(n))
forrest@0:
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ;(A)-2 html table of site933 -- observed
forrest@0: ;*******************************************************************
forrest@0: output_html = "table_site933_ob.html"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,"EMDI Observations Class B (933 sites)
" \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," | Site ID | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@1: ," NPP (gC m-2 y-1) | " \
forrest@1: ," PPT (m y-1) | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: delete (lines)
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: nrow = dimsizes(id933)
forrest@0: do n = 0,nrow-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt1 = sprintf("%5.0f", id933(n))
forrest@0: txt2 = sprintf("%5.2f", y933(n))
forrest@0: txt3 = sprintf("%5.2f", x933(n))
forrest@0: txt4 = sprintf("%5.2f", npp933(n))
forrest@0: txt5 = sprintf("%5.2f", rain933(n))
forrest@0:
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: delete (idx)
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ;(A)-3 html table of site81 -- model vs ob
forrest@0: ;*******************************************************************
forrest@0: output_html = "table_site81_model_vs_ob.html"
forrest@0:
forrest@0: header_text = "Model "+model_name+" vs Class A Observations (81 sites)
"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,header_text \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: delete (table_header)
forrest@0: table_header_text = " "+model_name+" | "
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," | Site ID | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@1: ," NPP (gC m-2 y-1) | " \
forrest@1: ," RAIN (m y-1) | " \
forrest@0: ,"
" \
forrest@0: ,"" \
forrest@0: ," | observed | " \
forrest@0: , table_header_text \
forrest@0: ," observed | " \
forrest@0: , table_header_text \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: delete (lines)
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: nrow = dimsizes(id81)
forrest@0: do n = 0,nrow-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt1 = sprintf("%5.0f", id81(n))
forrest@0: txt2 = sprintf("%5.2f", y81(n))
forrest@0: txt3 = sprintf("%5.2f", x81(n))
forrest@0: txt4 = sprintf("%5.2f", npp81(n))
forrest@0: txt5 = sprintf("%5.2f", nppmod81(n))
forrest@0: txt6 = sprintf("%5.2f", rain81(n))
forrest@0: txt7 = sprintf("%5.2f", rainmod81(n))
forrest@0:
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0: set_line(lines,nline,""+txt7+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: delete (idx)
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ;(A)-4 html table of site933 -- model vs ob
forrest@0: ;*******************************************************************
forrest@0: output_html = "table_site933_model_vs_ob.html"
forrest@0:
forrest@0: header_text = "Model "+model_name+" vs Class B Observations (933 sites)
"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,header_text \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: delete (table_header)
forrest@0: table_header_text = " "+model_name+" | "
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," | Site ID | " \
forrest@0: ," Latitude | " \
forrest@0: ," Longitude | " \
forrest@1: ," NPP (gC m-2 y-1) | " \
forrest@1: ," RAIN (m y-1) | " \
forrest@0: ,"
" \
forrest@0: ,"" \
forrest@0: ," | observed | " \
forrest@0: , table_header_text \
forrest@0: ," observed | " \
forrest@0: , table_header_text \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: delete (lines)
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: nrow = dimsizes(id933)
forrest@0: do n = 0,nrow-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt1 = sprintf("%5.0f", id933(n))
forrest@0: txt2 = sprintf("%5.2f", y933(n))
forrest@0: txt3 = sprintf("%5.2f", x933(n))
forrest@0: txt4 = sprintf("%5.2f", npp933(n))
forrest@0: txt5 = sprintf("%5.2f", nppmod933(n))
forrest@0: txt6 = sprintf("%5.2f", rain933(n))
forrest@0: txt7 = sprintf("%5.2f", rainmod933(n))
forrest@0:
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0: set_line(lines,nline,""+txt7+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: delete (idx)
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ;(A)-5 scatter plot, model vs ob 81
forrest@0: ;***************************************************************************
forrest@0:
forrest@0: plot_name = "scatter_model_vs_ob_81"
forrest@1: ;title = model_name +" vs Class A observations (81 sites)"
forrest@1: title = model_name +" (1975-2000) vs Class A observations (81 sites)"
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: res = True ; plot mods desired
forrest@0: res@tiMainString = title ; add title
forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers
forrest@0: res@xyMarkers = 16 ; choose type of marker
forrest@0: res@xyMarkerColor = "red" ; Marker color
forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels
forrest@0:
forrest@0: res@gsnDraw = False
forrest@0: res@gsnFrame = False ; don't advance frame yet
forrest@0: ;-------------------------------
forrest@0: ;compute correlation coef. and M
forrest@0: ccr81 = esccr(nppmod81,npp81,0)
forrest@0: ;print (ccr81)
forrest@0:
forrest@1: score_max = 1.0
forrest@0:
forrest@0: bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
forrest@0: M81s = (1. - (bias/dimsizes(y81)))*score_max
forrest@0: M_npp_S81 = sprintf("%.2f", M81s)
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_npp_S81/s/M_npp_S81/"+M_npp_S81+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: tRes = True
forrest@0: tRes@txFontHeightF = 0.025
forrest@0:
forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
forrest@0:
forrest@0: gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
forrest@0: ;--------------------------------
forrest@0: plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
forrest@0: ;-------------------------------
forrest@0: ; add polyline
forrest@0:
forrest@0: dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
forrest@0: ;-------------------------------
forrest@0: draw (plot)
forrest@0: delete (wks)
forrest@0: delete (dum)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ;(A)-6 scatter plot, model vs ob 933
forrest@0: ;***************************************************************************
forrest@0:
forrest@0: plot_name = "scatter_model_vs_ob_933"
forrest@1: title = model_name +" (1975-2000) vs Class B Observations (933 sites)"
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: res = True ; plot mods desired
forrest@0: res@tiMainString = title ; add title
forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers
forrest@0: res@xyMarkers = 16 ; choose type of marker
forrest@0: res@xyMarkerColor = "red" ; Marker color
forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels
forrest@0:
forrest@0: res@gsnDraw = False
forrest@0: res@gsnFrame = False ; don't advance frame yet
forrest@0: ;-------------------------------
forrest@0: ;compute correlation coef. and M
forrest@0:
forrest@0: ccr933 = esccr(nppmod933,npp933,0)
forrest@0:
forrest@1: score_max = 1.0
forrest@0:
forrest@0: bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
forrest@0: M933s = (1. - (bias/dimsizes(y933)))*score_max
forrest@0: M_npp_S933 = sprintf("%.2f", M933s)
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_npp_S933/s/M_npp_S933/"+M_npp_S933+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: tRes = True
forrest@0: tRes@txFontHeightF = 0.025
forrest@0:
forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
forrest@0:
forrest@0: gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
forrest@0: ;--------------------------------
forrest@0: plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
forrest@0: ;-------------------------------
forrest@0: ; add polyline
forrest@0:
forrest@0: dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
forrest@0: ;-------------------------------
forrest@0: draw (plot)
forrest@0: delete (wks)
forrest@0: delete (dum)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;**************************************************************************
forrest@0: ;(B) histogram 81
forrest@0: ;**************************************************************************
forrest@0: ;get data
forrest@0: RAIN1_1D = ndtooned(rain81)
forrest@0: RAIN2_1D = ndtooned(rainmod81)
forrest@0: NPP1_1D = ndtooned(npp81)
forrest@0: NPP2_1D = ndtooned(nppmod81)
forrest@0:
forrest@0: ; number of bin
forrest@0: nx = 8
forrest@0:
forrest@0: xvalues = new((/2,nx/),float)
forrest@0: yvalues = new((/2,nx/),float)
forrest@0: mn_yvalues = new((/2,nx/),float)
forrest@0: mx_yvalues = new((/2,nx/),float)
forrest@0: dx4 = new((/1/),float)
forrest@0:
forrest@0: get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
forrest@0: ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
forrest@0:
forrest@0: ;----------------------------------------
forrest@0: ;compute correlation coeff and M score
forrest@0:
forrest@0: u = yvalues(0,:)
forrest@0: v = yvalues(1,:)
forrest@0:
forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0: uu = u(good)
forrest@0: vv = v(good)
forrest@0:
forrest@0: ccr81h = esccr(uu,vv,0)
forrest@0:
forrest@1: score_max = 2.0
forrest@0:
forrest@0: bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
forrest@0: M81h = (1.- (bias/dimsizes(uu)))*score_max
forrest@0: M_npp_H81 = sprintf("%.2f", M81h)
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_npp_H81/s/M_npp_H81/"+M_npp_H81+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: delete (u)
forrest@0: delete (v)
forrest@0: delete (uu)
forrest@0: delete (vv)
forrest@0: ;----------------------------------------------------------------------
forrest@0: ; histogram res
forrest@0: resh = True
forrest@0: resh@gsnMaximize = True
forrest@0: resh@gsnDraw = False
forrest@0: resh@gsnFrame = False
forrest@0: resh@xyMarkLineMode = "Markers"
forrest@0: resh@xyMarkerSizeF = 0.014
forrest@0: resh@xyMarker = 16
forrest@0: resh@xyMarkerColors = (/"Brown","Blue"/)
forrest@0: resh@trYMinF = min(mn_yvalues) - 10.
forrest@0: resh@trYMaxF = max(mx_yvalues) + 10.
forrest@0:
forrest@1: resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: resh@tiXAxisString = "Precipitation (m y~S~-1~N~)"
forrest@0:
forrest@0: max_bar = new((/2,nx/),graphic)
forrest@0: min_bar = new((/2,nx/),graphic)
forrest@0: max_cap = new((/2,nx/),graphic)
forrest@0: min_cap = new((/2,nx/),graphic)
forrest@0:
forrest@0: lnres = True
forrest@0: line_colors = (/"brown","blue"/)
forrest@0:
forrest@0: ;**************************************************************************
forrest@0: ;(B)-1 histogram plot, ob 81 site
forrest@0: ;**************************************************************************
forrest@0:
forrest@0: plot_name = "histogram_ob_81"
forrest@0: title = "Class A Observations (81 sites)"
forrest@0: resh@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name)
forrest@0:
forrest@0: xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
forrest@0:
forrest@1: ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
forrest@1: print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))
forrest@1:
forrest@0: ;-------------------------------
forrest@0: ;Attach the vertical bar and the horizontal cap line
forrest@0:
forrest@0: do nd=0,0
forrest@0: lnres@gsLineColor = line_colors(nd)
forrest@0: do i=0,nx-1
forrest@0:
forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0: .not.ismissing(mx_yvalues(nd,i))) then
forrest@0: ;
forrest@0: ; Attach the vertical bar, both above and below the marker.
forrest@0: ;
forrest@0: x1 = xvalues(nd,i)
forrest@0: y1 = yvalues(nd,i)
forrest@0: y2 = mn_yvalues(nd,i)
forrest@0: plx = (/x1,x1/)
forrest@0: ply = (/y1,y2/)
forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0:
forrest@0: y2 = mx_yvalues(nd,i)
forrest@0: plx = (/x1,x1/)
forrest@0: ply = (/y1,y2/)
forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0: ;
forrest@0: ; Attach the horizontal cap line, both above and below the marker.
forrest@0: ;
forrest@0: x1 = xvalues(nd,i) - dx4
forrest@0: x2 = xvalues(nd,i) + dx4
forrest@0: y1 = mn_yvalues(nd,i)
forrest@0: plx = (/x1,x2/)
forrest@0: ply = (/y1,y1/)
forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0:
forrest@0: y1 = mx_yvalues(nd,i)
forrest@0: plx = (/x1,x2/)
forrest@0: ply = (/y1,y1/)
forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
forrest@0: end if
forrest@0: end do
forrest@0: end do
forrest@0:
forrest@0: draw(xy)
forrest@0: delete (wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;****************************************************************************
forrest@0: ;(B)-2 histogram plot, model vs ob 81 site
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: plot_name = "histogram_model_vs_ob_81"
forrest@1: title = model_name+ " (1975-2000) vs Class A Observations (81 sites)"
forrest@0: resh@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0:
forrest@0: ;-----------------------------
forrest@0: ; Add a boxed legend using the more simple method, which won't have
forrest@0: ; vertical lines going through the markers.
forrest@0:
forrest@0: resh@pmLegendDisplayMode = "Always"
forrest@0: ; resh@pmLegendWidthF = 0.1
forrest@0: resh@pmLegendWidthF = 0.08
forrest@0: resh@pmLegendHeightF = 0.05
forrest@0: resh@pmLegendOrthogonalPosF = -1.17
forrest@0: ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0: ; resh@pmLegendParallelPosF = 0.18
forrest@0: resh@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0:
forrest@0: ; resh@lgPerimOn = False
forrest@0: resh@lgLabelFontHeightF = 0.015
forrest@0: resh@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0: ;-----------------------------
forrest@0: tRes = True
forrest@0: tRes@txFontHeightF = 0.025
forrest@0:
forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
forrest@0:
forrest@0: gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0:
forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
forrest@1:
forrest@1: ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
forrest@1: print(model_name + ": " + xvalues(1,:) + " " + yvalues(1,:))
forrest@1: ;print("All: " + xvalues(:,:) + " " + yvalues(:,:))
forrest@0: ;-------------------------------
forrest@0: ;Attach the vertical bar and the horizontal cap line
forrest@0:
forrest@0: do nd=0,1
forrest@0: lnres@gsLineColor = line_colors(nd)
forrest@0: do i=0,nx-1
forrest@0:
forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0: .not.ismissing(mx_yvalues(nd,i))) then
forrest@0: ;
forrest@0: ; Attach the vertical bar, both above and below the marker.
forrest@0: ;
forrest@0: x1 = xvalues(nd,i)
forrest@0: y1 = yvalues(nd,i)
forrest@0: y2 = mn_yvalues(nd,i)
forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0:
forrest@0: y2 = mx_yvalues(nd,i)
forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0: ;
forrest@0: ; Attach the horizontal cap line, both above and below the marker.
forrest@0: ;
forrest@0: x1 = xvalues(nd,i) - dx4
forrest@0: x2 = xvalues(nd,i) + dx4
forrest@0: y1 = mn_yvalues(nd,i)
forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0:
forrest@0: y1 = mx_yvalues(nd,i)
forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0: end if
forrest@0: end do
forrest@0: end do
forrest@0:
forrest@0: draw(xy)
forrest@0: delete(wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: delete (RAIN1_1D)
forrest@0: delete (RAIN2_1D)
forrest@0: delete (NPP1_1D)
forrest@0: delete (NPP2_1D)
forrest@0: ;delete (range)
forrest@0: delete (xvalues)
forrest@0: delete (yvalues)
forrest@0: delete (mn_yvalues)
forrest@0: delete (mx_yvalues)
forrest@0: delete (good)
forrest@0: delete (max_bar)
forrest@0: delete (min_bar)
forrest@0: delete (max_cap)
forrest@0: delete (min_cap)
forrest@0:
forrest@0: ;**************************************************************************
forrest@0: ;(B) histogram 933
forrest@0: ;**************************************************************************
forrest@0:
forrest@0: ; get data
forrest@0: RAIN1_1D = ndtooned(rain933)
forrest@0: RAIN2_1D = ndtooned(rainmod933)
forrest@0: NPP1_1D = ndtooned(npp933)
forrest@0: NPP2_1D = ndtooned(nppmod933)
forrest@0:
forrest@0: ; number of bin
forrest@0: nx = 10
forrest@0:
forrest@0: xvalues = new((/2,nx/),float)
forrest@0: yvalues = new((/2,nx/),float)
forrest@0: mn_yvalues = new((/2,nx/),float)
forrest@0: mx_yvalues = new((/2,nx/),float)
forrest@0: dx4 = new((/1/),float)
forrest@0:
forrest@0: get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
forrest@0: ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
forrest@0:
forrest@0: ;----------------------------------------
forrest@0: ;compute correlation coeff and M score
forrest@0:
forrest@0: u = yvalues(0,:)
forrest@0: v = yvalues(1,:)
forrest@0:
forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0: uu = u(good)
forrest@0: vv = v(good)
forrest@0:
forrest@0: ccr933h = esccr(uu,vv,0)
forrest@0:
forrest@1: score_max = 2.0
forrest@0:
forrest@0: bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
forrest@0: M933h = (1.- (bias/dimsizes(uu)))*score_max
forrest@0: M_npp_H933 = sprintf("%.2f", M933h)
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_npp_H933/s/M_npp_H933/"+M_npp_H933+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: delete (u)
forrest@0: delete (v)
forrest@0: delete (uu)
forrest@0: delete (vv)
forrest@0: ;----------------------------------------------------------------------
forrest@0: ; histogram res
forrest@0: delete (resh)
forrest@0: resh = True
forrest@0: resh@gsnMaximize = True
forrest@0: resh@gsnDraw = False
forrest@0: resh@gsnFrame = False
forrest@0: resh@xyMarkLineMode = "Markers"
forrest@0: resh@xyMarkerSizeF = 0.014
forrest@0: resh@xyMarker = 16
forrest@0: resh@xyMarkerColors = (/"Brown","Blue"/)
forrest@0: resh@trYMinF = min(mn_yvalues) - 10.
forrest@0: resh@trYMaxF = max(mx_yvalues) + 10.
forrest@0:
forrest@1: resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@1: resh@tiXAxisString = "Precipitation (m y~S~-1~N~)"
forrest@0:
forrest@0: max_bar = new((/2,nx/),graphic)
forrest@0: min_bar = new((/2,nx/),graphic)
forrest@0: max_cap = new((/2,nx/),graphic)
forrest@0: min_cap = new((/2,nx/),graphic)
forrest@0:
forrest@0: lnres = True
forrest@0: line_colors = (/"brown","blue"/)
forrest@0:
forrest@0: ;**************************************************************************
forrest@0: ;(B)-3 histogram plot, ob 933 site
forrest@0: ;**************************************************************************
forrest@0:
forrest@0: plot_name = "histogram_ob_933"
forrest@0: title = "Class B Observations (933 sites)"
forrest@0: resh@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name)
forrest@0:
forrest@0: xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
forrest@0:
forrest@1: ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
forrest@1: print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))
forrest@1:
forrest@0: ;-------------------------------
forrest@0: ;Attach the vertical bar and the horizontal cap line
forrest@0:
forrest@0: do nd=0,0
forrest@0: lnres@gsLineColor = line_colors(nd)
forrest@0: do i=0,nx-1
forrest@0:
forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0: .not.ismissing(mx_yvalues(nd,i))) then
forrest@0:
forrest@0: ; Attach the vertical bar, both above and below the marker.
forrest@0:
forrest@0: x1 = xvalues(nd,i)
forrest@0: y1 = yvalues(nd,i)
forrest@0: y2 = mn_yvalues(nd,i)
forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0:
forrest@0: y2 = mx_yvalues(nd,i)
forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0:
forrest@0: ; Attach the horizontal cap line, both above and below the marker.
forrest@0:
forrest@0: x1 = xvalues(nd,i) - dx4
forrest@0: x2 = xvalues(nd,i) + dx4
forrest@0: y1 = mn_yvalues(nd,i)
forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0:
forrest@0: y1 = mx_yvalues(nd,i)
forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0: end if
forrest@0: end do
forrest@0: end do
forrest@0:
forrest@0: draw(xy)
forrest@0: delete (xy)
forrest@0: delete (wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;****************************************************************************
forrest@0: ;(B)-4 histogram plot, model vs ob 933 site
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: plot_name = "histogram_model_vs_ob_933"
forrest@1: title = model_name+ " (1975-2000) vs Class B Observations (933 sites)"
forrest@0: resh@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0:
forrest@0: ;-----------------------------
forrest@0: ; Add a boxed legend using the more simple method, which won't have
forrest@0: ; vertical lines going through the markers.
forrest@0:
forrest@0: resh@pmLegendDisplayMode = "Always"
forrest@0: ; resh@pmLegendWidthF = 0.1
forrest@0: resh@pmLegendWidthF = 0.08
forrest@0: resh@pmLegendHeightF = 0.05
forrest@0: resh@pmLegendOrthogonalPosF = -1.17
forrest@0: ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
forrest@0: ; resh@pmLegendParallelPosF = 0.18
forrest@0: resh@pmLegendParallelPosF = 0.88 ;(rightward)
forrest@0:
forrest@0: ; resh@lgPerimOn = False
forrest@0: resh@lgLabelFontHeightF = 0.015
forrest@0: resh@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0: ;-----------------------------
forrest@0: tRes = True
forrest@0: tRes@txFontHeightF = 0.025
forrest@0:
forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
forrest@0:
forrest@0: gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
forrest@0:
forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
forrest@1:
forrest@1: ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
forrest@1: print("Observations: " + xvalues(1,:) + " " + yvalues(1,:))
forrest@0: ;-------------------------------
forrest@0: ;Attach the vertical bar and the horizontal cap line
forrest@0:
forrest@0: do nd=0,1
forrest@0: lnres@gsLineColor = line_colors(nd)
forrest@0: do i=0,nx-1
forrest@0:
forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \
forrest@0: .not.ismissing(mx_yvalues(nd,i))) then
forrest@0: ;
forrest@0: ; Attach the vertical bar, both above and below the marker.
forrest@0: ;
forrest@0: x1 = xvalues(nd,i)
forrest@0: y1 = yvalues(nd,i)
forrest@0: y2 = mn_yvalues(nd,i)
forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0:
forrest@0: y2 = mx_yvalues(nd,i)
forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
forrest@0: ;
forrest@0: ; Attach the horizontal cap line, both above and below the marker.
forrest@0: ;
forrest@0: x1 = xvalues(nd,i) - dx4
forrest@0: x2 = xvalues(nd,i) + dx4
forrest@0: y1 = mn_yvalues(nd,i)
forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0:
forrest@0: y1 = mx_yvalues(nd,i)
forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
forrest@0: end if
forrest@0: end do
forrest@0: end do
forrest@0:
forrest@0: draw(xy)
forrest@0: delete(xy)
forrest@0: delete(wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@1:
forrest@1: ;***************************************************************************
forrest@1: ; Read 2000-2004 dataset for MODIS comparison
forrest@1: ;***************************************************************************
forrest@1:
forrest@1: ;------------------------------------------------------
forrest@1: ; read model data
forrest@1:
forrest@1: fm = addfile (dirm+film9,"r")
forrest@1:
forrest@1: nppmod0 = fm->NPP
forrest@1: xm = fm->lon
forrest@1: ym = fm->lat
forrest@1:
forrest@1: delete (fm)
forrest@1:
forrest@1: nppmod = nppmod0(0,:,:)
forrest@1: delete (nppmod0)
forrest@1: nppmod = nppmod * nsec_per_year
forrest@1: nppmod@units = "gC m~S~-2~N~ y~S~-1~N~"
forrest@1: nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ;(C) global contour
forrest@0: ;***************************************************************************
forrest@0:
forrest@0: ;res
forrest@0: resg = True ; Use plot options
forrest@0: resg@cnFillOn = True ; Turn on color fill
forrest@0: resg@gsnSpreadColors = True ; use full colormap
forrest@0: ; resg@cnFillMode = "RasterFill" ; Turn on raster color
forrest@0: ; resg@lbLabelAutoStride = True
forrest@0: resg@cnLinesOn = False ; Turn off contourn lines
forrest@0: resg@mpFillOn = False ; Turn off map fill
forrest@0:
forrest@0: resg@gsnSpreadColors = True ; use full colormap
forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
forrest@0: resg@cnMinLevelValF = 0. ; Min level
forrest@0: resg@cnMaxLevelValF = 2200. ; Max level
forrest@0: resg@cnLevelSpacingF = 200. ; interval
forrest@0:
forrest@0: ;****************************************************************************
forrest@0: ;(C)-1 global contour plot, ob
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: delta = 0.00000000001
forrest@0: nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
forrest@0:
forrest@0: plot_name = "global_ob"
forrest@1: title = "MODIS MOD17A3 (2000-2004)"
forrest@0: resg@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0:
forrest@0: plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
forrest@0:
forrest@0: delete (wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0: ;****************************************************************************
forrest@0: ;(C)-2 global contour plot, model
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: plot_name = "global_model"
forrest@1: title = "Model "+ model_name + " (2000-2004)"
forrest@0: resg@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0:
forrest@0: plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
forrest@0:
forrest@0: delete (wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;****************************************************************************
forrest@0: ;(C)-3 global contour plot, model vs ob
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: plot_name = "global_model_vs_ob"
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap
forrest@0:
forrest@0: delete (plot)
forrest@0: plot=new(3,graphic) ; create graphic array
forrest@0:
forrest@0: resg@gsnFrame = False ; Do not draw plot
forrest@0: resg@gsnDraw = False ; Do not advance frame
forrest@0:
forrest@0: ; compute correlation coef and M score
forrest@0:
forrest@0: uu1 = ndtooned(nppmod)
forrest@0: vv1 = ndtooned(nppglobe)
forrest@0:
forrest@0: delete (good)
forrest@0: good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
forrest@0:
forrest@0: ug = uu1(good)
forrest@0: vg = vv1(good)
forrest@0:
forrest@0: ccrG = esccr(ug,vg,0)
forrest@0:
forrest@1: score_max = 2.0
forrest@0:
forrest@0: MG = (ccrG*ccrG)* score_max
forrest@0: M_npp_G = sprintf("%.2f", MG)
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_npp_G/s/M_npp_G/"+M_npp_G+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: ; plot correlation coef
forrest@0:
forrest@0: gRes = True
forrest@0: gRes@txFontHeightF = 0.02
forrest@0: gRes@txAngleF = 90
forrest@0:
forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
forrest@0:
forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0: ;--------------------------------------------------------------------
forrest@0: ;(a) ob
forrest@0:
forrest@1: title = "MODIS MOD17A3 (2000-2004)"
forrest@0: resg@tiMainString = title
forrest@0:
forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
forrest@0:
forrest@0: ;(b) model
forrest@0:
forrest@1: title = "Model "+ model_name + " (2000-2004)"
forrest@0: resg@tiMainString = title
forrest@0:
forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
forrest@0:
forrest@0: ;(c) model-ob
forrest@0:
forrest@0: zz = nppmod
forrest@0: zz = nppmod - nppglobe
forrest@1: title = "Model "+model_name+" - MODIS MOD17A3"
forrest@0:
forrest@0: resg@cnMinLevelValF = -500 ; Min level
forrest@0: resg@cnMaxLevelValF = 500. ; Max level
forrest@0: resg@cnLevelSpacingF = 50. ; interval
forrest@0: resg@tiMainString = title
forrest@0:
forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
forrest@0:
forrest@0: pres = True ; panel plot mods desired
forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
forrest@0: ; indiv. plots in panel
forrest@0: pres@gsnMaximize = True ; fill the page
forrest@0:
forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
forrest@0:
forrest@0: delete (wks)
forrest@0: delete (plot)
forrest@0:
forrest@1: system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ;(D)-1 zonal line plot, ob
forrest@0: ;***************************************************************************
forrest@0:
forrest@0: vv = zonalAve(nppglobe)
forrest@0: vv@long_name = nppglobe@long_name
forrest@0:
forrest@0: plot_name = "zonal_ob"
forrest@1: title = "MODIS MOD17A3 (2000-2004)"
forrest@0:
forrest@0: resz = True
forrest@0: resz@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation
forrest@0:
forrest@0: plot = gsn_csm_xy (wks,ym,vv,resz)
forrest@0:
forrest@0: delete (wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;****************************************************************************
forrest@0: ;(D)-2 zonal line plot, model vs ob
forrest@0: ;****************************************************************************
forrest@0:
forrest@0: s = new ((/2,dimsizes(ym)/), float)
forrest@0:
forrest@0: s(0,:) = zonalAve(nppglobe)
forrest@0: s(1,:) = zonalAve(nppmod)
forrest@0:
forrest@0: s@long_name = nppglobe@long_name
forrest@0: ;-------------------------------------------
forrest@0: ; compute correlation coef and M score
forrest@0:
forrest@1: score_max = 2.0
forrest@0:
forrest@0: ccrZ = esccr(s(1,:), s(0,:),0)
forrest@0: MZ = (ccrZ*ccrZ)* score_max
forrest@0: M_npp_Z = sprintf("%.2f", MZ)
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_npp_Z/s/M_npp_Z/"+M_npp_Z+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0: ;-------------------------------------------
forrest@0: plot_name = "zonal_model_vs_ob"
forrest@1: title = "Zonal Average (2000-2004)"
forrest@0: resz@tiMainString = title
forrest@0:
forrest@0: wks = gsn_open_wks (plot_type,plot_name)
forrest@0:
forrest@0: ; resz@vpHeightF = 0.4 ; change aspect ratio of plot
forrest@0: ; resz@vpWidthF = 0.7
forrest@0:
forrest@0: resz@xyMonoLineColor = "False" ; want colored lines
forrest@0: resz@xyLineColors = (/"black","red"/) ; colors chosen
forrest@0: ; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
forrest@0: resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
forrest@0: resz@xyDashPatterns = (/0.,0./) ; make all lines solid
forrest@0:
forrest@0: resz@tiYAxisString = s@long_name ; add a axis title
forrest@0: resz@txFontHeightF = 0.0195 ; change title font heights
forrest@0:
forrest@0: ; Legent
forrest@0: resz@pmLegendDisplayMode = "Always" ; turn on legend
forrest@0: resz@pmLegendSide = "Top" ; Change location of
forrest@0: ; resz@pmLegendParallelPosF = .45 ; move units right
forrest@1: resz@pmLegendParallelPosF = .75 ; move units right
forrest@0: resz@pmLegendOrthogonalPosF = -0.4 ; move units down
forrest@0:
forrest@0: resz@pmLegendWidthF = 0.10 ; Change width and
forrest@0: resz@pmLegendHeightF = 0.10 ; height of legend.
forrest@1: resz@lgLabelFontHeightF = .015 ; change font height
forrest@0: ; resz@lgTitleOn = True ; turn on legend title
forrest@0: ; resz@lgTitleString = "Example" ; create legend title
forrest@0: ; resz@lgTitleFontHeightF = .025 ; font of legend title
forrest@1: resz@xyExplicitLegendLabels = (/"MODIS MOD17A3",model_name/) ; explicit labels
forrest@0: ;--------------------------------------------------------------------
forrest@0: zRes = True
forrest@0: zRes@txFontHeightF = 0.025
forrest@0:
forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
forrest@0:
forrest@0: gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
forrest@0: ;--------------------------------------------------------------------
forrest@0:
forrest@0: plot = gsn_csm_xy (wks,ym,s,resz)
forrest@0:
forrest@0: delete (wks)
forrest@0:
forrest@1: system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0: "rm "+plot_name+"."+plot_type)
forrest@1: ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ; add total score and write to file
forrest@0: ;***************************************************************************
forrest@0: M_total = M81s + M81h + M933s + M933h + MG + MZ
forrest@0:
forrest@0: asciiwrite("M_save.npp", M_total)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ; output plots
forrest@0: ;***************************************************************************
forrest@0: output_dir = model_name+"/npp"
forrest@0:
forrest@0: system("mv table_*.html " + output_dir +";"+ \
forrest@0: "mv *.png " + output_dir)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: end
forrest@0: