diff -r 000000000000 -r 0c6405ab2ff4 npp/15.scatter_bias.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/npp/15.scatter_bias.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,285 @@
+; ***********************************************
+; combine all scatter
+; ***********************************************
+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"
+;************************************************
+begin
+
+ plot_type     = "ps"
+ plot_type_new = "png"
+
+;************************************************
+; read in ob data
+;************************************************
+
+;(A) plot ob scatter_81
+
+ dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
+ fil81 = "data.81.nc"
+ f81   = addfile (dir81+fil81,"r")
+
+ x81   = f81->SITE_ID  
+ y81   = f81->TNPP_C
+ xo81  = f81->LONG_DD  
+ yo81  = f81->LAT_DD
+
+ x81@long_name = "SITE_ID"
+ y81@long_name = "NPP (gC/m2/year)"
+   
+ plot_name = "scatter_ob_81"
+ title     = "Observed 81 sites"
+
+ wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
+ res                   = True                  ; plot mods desired
+ res@tiMainString      = title                 ; add title
+ res@xyMarkLineModes   = "Markers"             ; choose which have markers
+ res@xyMarkers         =  16                   ; choose type of marker  
+ res@xyMarkerColor     = "red"                 ; Marker color
+ res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True                  ; nice tick mark labels
+
+ plot  = gsn_csm_xy (wks,x81,y81,res)              ; create plot
+ frame(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+ system("rm "+plot_name+"."+plot_type)
+ system("rm "+plot_name+"-1."+plot_type_new)
+ system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
+
+ clear (wks)
+ 
+;(B) plot ob scatter_933
+
+ dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
+ fil933 = "data.933.nc"
+ f933   = addfile (dir933+fil933,"r")
+
+ x933  = f933->SITE_ID  
+ y933  = f933->TNPP_C
+ xo933 = f933->LONG_DD  
+ yo933 = f933->LAT_DD 
+
+ x933@long_name = "SITE_ID"
+ y933@long_name = "NPP (gC/m2/year)"
+
+ plot_name = "scatter_ob_933"
+ title     = "Observed 933 sites"
+
+ wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
+ res                   = True                  ; plot mods desired
+ res@tiMainString      = title                 ; add title
+ res@xyMarkLineModes   = "Markers"             ; choose which have markers
+ res@xyMarkers         =  16                   ; choose type of marker  
+ res@xyMarkerColor     = "red"                 ; Marker color
+ res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True                  ; nice tick mark labels
+
+ plot  = gsn_csm_xy (wks,x933,y933,res)              ; create plot
+ frame(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+ system("rm "+plot_name+"."+plot_type)
+ system("rm "+plot_name+"-1."+plot_type_new)
+ system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
+
+ clear (wks)
+;************************************************
+; read in model data
+;************************************************
+
+;(C) plot model scatter_81
+
+ model_name = "b30.061n"
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
+ fm   = addfile (dirm+film,"r")
+  
+ ymod = fm->NPP
+ xm   = fm->lon  
+ ym   = fm->lat
+
+ nx = dimsizes(xo81)
+ do i= 0,nx-1
+    if (xo81(i) .lt. 0.) then
+        xo81(i) = xo81(i)+ 360.
+    end if
+ end do
+ print (xo81)
+
+ sec_to_year = 86400.*365.
+
+ ymod81 = linint2_points(xm,ym,ymod(0,:,:),True,xo81,yo81,0) * sec_to_year
+ xmod81 = x81
+
+ ymod81@long_name = "NPP (gC/m2/year)"
+ xmod81@long_name = "SITE_ID"
+print (ymod81)
+print (xmod81)
+  
+ plot_name = "scatter_model_81"
+ title     = "Model "+ model_name +" 81 sites"
+
+ wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
+ res                   = True                  ; plot mods desired
+ res@tiMainString      = title                 ; add title
+ res@xyMarkLineModes   = "Markers"             ; choose which have markers
+ res@xyMarkers         =  16                   ; choose type of marker  
+ res@xyMarkerColor     = "red"                 ; Marker color
+ res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True                  ; nice tick mark labels
+
+ plot  = gsn_csm_xy (wks,xmod81,ymod81,res)    ; create plot
+ frame(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+ system("rm "+plot_name+"."+plot_type)
+ system("rm "+plot_name+"-1."+plot_type_new)
+ system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
+
+ clear (wks)
+ 
+;(D) plot model scatter_933
+
+ model_name = "b30.061n"
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
+ fm   = addfile (dirm+film,"r")
+  
+ ymod = fm->NPP
+ xm   = fm->lon  
+ ym   = fm->lat
+
+ nx = dimsizes(xo933)
+ do i= 0,nx-1
+    if (xo933(i) .lt. 0.) then
+        xo933(i) = xo933(i)+ 360.
+    end if
+ end do
+ print (xo933)
+
+ sec_to_year = 86400.*365.
+
+ ymod933 = linint2_points(xm,ym,ymod(0,:,:),True,xo933,yo933,0) * sec_to_year
+ xmod933 = x933
+
+ ymod933@long_name = "NPP (gC/m2/year)"
+ xmod933@long_name = "SITE_ID"
+;print (ymod933)
+;print (xmod933)
+  
+ plot_name = "scatter_model_933"
+ title     = "Model "+ model_name +" 933 sites"
+
+ wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
+ res                   = True                  ; plot mods desired
+ res@tiMainString      = title                 ; add title
+ res@xyMarkLineModes   = "Markers"             ; choose which have markers
+ res@xyMarkers         =  16                   ; choose type of marker  
+ res@xyMarkerColor     = "red"                 ; Marker color
+ res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True                  ; nice tick mark labels
+
+ plot  = gsn_csm_xy (wks,xmod933,ymod933,res)    ; create plot
+ frame(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+ system("rm "+plot_name+"."+plot_type)
+ system("rm "+plot_name+"-1."+plot_type_new)
+ system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
+
+ clear (wks)
+
+;(E) scatter model vs ob 81
+
+ ymod81@long_name = "NPP (gC/m2/year)"
+ y81@long_name    = "NPP (gC/m2/year)"
+;print (ymod81)
+;print (y81)
+  
+ plot_name = "scatter_model_vs_ob_81"
+ title     = model_name +" vs ob 81 sites"
+
+ wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
+ res                   = True                  ; plot mods desired
+ res@tiMainString      = title                 ; add title
+ res@xyMarkLineModes   = "Markers"             ; choose which have markers
+ res@xyMarkers         =  16                   ; choose type of marker  
+ res@xyMarkerColor     = "red"                 ; Marker color
+ res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True                  ; nice tick mark labels
+;-------------------------------
+;compute correlation coef. and M
+ ccr81 = esccr(ymod81,y81,0)
+ print (ccr81)
+
+;new eq
+ bias = sum(abs(ymod81-y81)/(ymod81+y81))
+ M81  = (1. - (bias/dimsizes(y81)))*5.
+ print (M81)
+
+ tRes  = True
+ tRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
+;--------------------------------
+ plot  = gsn_csm_xy (wks,y81,ymod81,res)       ; create plot
+ frame(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+ system("rm "+plot_name+"."+plot_type)
+ system("rm "+plot_name+"-1."+plot_type_new)
+ system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
+
+ clear (wks)
+ 
+;(F) scatter model vs ob 933
+ 
+ ymod933@long_name = "NPP (gC/m2/year)"
+ y933@long_name    = "NPP (gC/m2/year)"
+;print (ymod33)
+;print (y933)
+  
+ plot_name = "scatter_model_vs_ob_933"
+ title     = model_name +" vs ob 933 sites"
+
+ wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
+ res                   = True                  ; plot mods desired
+ res@tiMainString      = title                 ; add title
+ res@xyMarkLineModes   = "Markers"             ; choose which have markers
+ res@xyMarkers         =  16                   ; choose type of marker  
+ res@xyMarkerColor     = "red"                 ; Marker color
+ res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True                  ; nice tick mark labels
+;-------------------------------
+;compute correlation coef. and M
+ ccr933 = esccr(ymod933,y933,0)
+ print (ccr933)
+
+;new eq
+ bias = sum(abs(ymod933-y933)/(ymod933+y933))
+ M933 = (1. - (bias/dimsizes(y933)))*5.
+ print (M933)
+
+ tRes  = True
+ tRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
+;--------------------------------
+ plot  = gsn_csm_xy (wks,y933,ymod933,res)    ; create plot
+ frame(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+ system("rm "+plot_name+"."+plot_type)
+ system("rm "+plot_name+"-1."+plot_type_new)
+ system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
+
+ clear (wks)
+end