1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/15.scatter_bias.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,285 @@
1.4 +; ***********************************************
1.5 +; combine all scatter
1.6 +; ***********************************************
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.10 +;************************************************
1.11 +begin
1.12 +
1.13 + plot_type = "ps"
1.14 + plot_type_new = "png"
1.15 +
1.16 +;************************************************
1.17 +; read in ob data
1.18 +;************************************************
1.19 +
1.20 +;(A) plot ob scatter_81
1.21 +
1.22 + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.23 + fil81 = "data.81.nc"
1.24 + f81 = addfile (dir81+fil81,"r")
1.25 +
1.26 + x81 = f81->SITE_ID
1.27 + y81 = f81->TNPP_C
1.28 + xo81 = f81->LONG_DD
1.29 + yo81 = f81->LAT_DD
1.30 +
1.31 + x81@long_name = "SITE_ID"
1.32 + y81@long_name = "NPP (gC/m2/year)"
1.33 +
1.34 + plot_name = "scatter_ob_81"
1.35 + title = "Observed 81 sites"
1.36 +
1.37 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.38 + res = True ; plot mods desired
1.39 + res@tiMainString = title ; add title
1.40 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.41 + res@xyMarkers = 16 ; choose type of marker
1.42 + res@xyMarkerColor = "red" ; Marker color
1.43 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.44 + res@tmLabelAutoStride = True ; nice tick mark labels
1.45 +
1.46 + plot = gsn_csm_xy (wks,x81,y81,res) ; create plot
1.47 + frame(wks)
1.48 +
1.49 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.50 + system("rm "+plot_name+"."+plot_type)
1.51 + system("rm "+plot_name+"-1."+plot_type_new)
1.52 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.53 +
1.54 + clear (wks)
1.55 +
1.56 +;(B) plot ob scatter_933
1.57 +
1.58 + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.59 + fil933 = "data.933.nc"
1.60 + f933 = addfile (dir933+fil933,"r")
1.61 +
1.62 + x933 = f933->SITE_ID
1.63 + y933 = f933->TNPP_C
1.64 + xo933 = f933->LONG_DD
1.65 + yo933 = f933->LAT_DD
1.66 +
1.67 + x933@long_name = "SITE_ID"
1.68 + y933@long_name = "NPP (gC/m2/year)"
1.69 +
1.70 + plot_name = "scatter_ob_933"
1.71 + title = "Observed 933 sites"
1.72 +
1.73 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.74 + res = True ; plot mods desired
1.75 + res@tiMainString = title ; add title
1.76 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.77 + res@xyMarkers = 16 ; choose type of marker
1.78 + res@xyMarkerColor = "red" ; Marker color
1.79 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.80 + res@tmLabelAutoStride = True ; nice tick mark labels
1.81 +
1.82 + plot = gsn_csm_xy (wks,x933,y933,res) ; create plot
1.83 + frame(wks)
1.84 +
1.85 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.86 + system("rm "+plot_name+"."+plot_type)
1.87 + system("rm "+plot_name+"-1."+plot_type_new)
1.88 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.89 +
1.90 + clear (wks)
1.91 +;************************************************
1.92 +; read in model data
1.93 +;************************************************
1.94 +
1.95 +;(C) plot model scatter_81
1.96 +
1.97 + model_name = "b30.061n"
1.98 +
1.99 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.100 + film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.101 + fm = addfile (dirm+film,"r")
1.102 +
1.103 + ymod = fm->NPP
1.104 + xm = fm->lon
1.105 + ym = fm->lat
1.106 +
1.107 + nx = dimsizes(xo81)
1.108 + do i= 0,nx-1
1.109 + if (xo81(i) .lt. 0.) then
1.110 + xo81(i) = xo81(i)+ 360.
1.111 + end if
1.112 + end do
1.113 + print (xo81)
1.114 +
1.115 + sec_to_year = 86400.*365.
1.116 +
1.117 + ymod81 = linint2_points(xm,ym,ymod(0,:,:),True,xo81,yo81,0) * sec_to_year
1.118 + xmod81 = x81
1.119 +
1.120 + ymod81@long_name = "NPP (gC/m2/year)"
1.121 + xmod81@long_name = "SITE_ID"
1.122 +print (ymod81)
1.123 +print (xmod81)
1.124 +
1.125 + plot_name = "scatter_model_81"
1.126 + title = "Model "+ model_name +" 81 sites"
1.127 +
1.128 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.129 + res = True ; plot mods desired
1.130 + res@tiMainString = title ; add title
1.131 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.132 + res@xyMarkers = 16 ; choose type of marker
1.133 + res@xyMarkerColor = "red" ; Marker color
1.134 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.135 + res@tmLabelAutoStride = True ; nice tick mark labels
1.136 +
1.137 + plot = gsn_csm_xy (wks,xmod81,ymod81,res) ; create plot
1.138 + frame(wks)
1.139 +
1.140 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.141 + system("rm "+plot_name+"."+plot_type)
1.142 + system("rm "+plot_name+"-1."+plot_type_new)
1.143 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.144 +
1.145 + clear (wks)
1.146 +
1.147 +;(D) plot model scatter_933
1.148 +
1.149 + model_name = "b30.061n"
1.150 +
1.151 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.152 + film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.153 + fm = addfile (dirm+film,"r")
1.154 +
1.155 + ymod = fm->NPP
1.156 + xm = fm->lon
1.157 + ym = fm->lat
1.158 +
1.159 + nx = dimsizes(xo933)
1.160 + do i= 0,nx-1
1.161 + if (xo933(i) .lt. 0.) then
1.162 + xo933(i) = xo933(i)+ 360.
1.163 + end if
1.164 + end do
1.165 + print (xo933)
1.166 +
1.167 + sec_to_year = 86400.*365.
1.168 +
1.169 + ymod933 = linint2_points(xm,ym,ymod(0,:,:),True,xo933,yo933,0) * sec_to_year
1.170 + xmod933 = x933
1.171 +
1.172 + ymod933@long_name = "NPP (gC/m2/year)"
1.173 + xmod933@long_name = "SITE_ID"
1.174 +;print (ymod933)
1.175 +;print (xmod933)
1.176 +
1.177 + plot_name = "scatter_model_933"
1.178 + title = "Model "+ model_name +" 933 sites"
1.179 +
1.180 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.181 + res = True ; plot mods desired
1.182 + res@tiMainString = title ; add title
1.183 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.184 + res@xyMarkers = 16 ; choose type of marker
1.185 + res@xyMarkerColor = "red" ; Marker color
1.186 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.187 + res@tmLabelAutoStride = True ; nice tick mark labels
1.188 +
1.189 + plot = gsn_csm_xy (wks,xmod933,ymod933,res) ; create plot
1.190 + frame(wks)
1.191 +
1.192 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.193 + system("rm "+plot_name+"."+plot_type)
1.194 + system("rm "+plot_name+"-1."+plot_type_new)
1.195 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.196 +
1.197 + clear (wks)
1.198 +
1.199 +;(E) scatter model vs ob 81
1.200 +
1.201 + ymod81@long_name = "NPP (gC/m2/year)"
1.202 + y81@long_name = "NPP (gC/m2/year)"
1.203 +;print (ymod81)
1.204 +;print (y81)
1.205 +
1.206 + plot_name = "scatter_model_vs_ob_81"
1.207 + title = model_name +" vs ob 81 sites"
1.208 +
1.209 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.210 + res = True ; plot mods desired
1.211 + res@tiMainString = title ; add title
1.212 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.213 + res@xyMarkers = 16 ; choose type of marker
1.214 + res@xyMarkerColor = "red" ; Marker color
1.215 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.216 + res@tmLabelAutoStride = True ; nice tick mark labels
1.217 +;-------------------------------
1.218 +;compute correlation coef. and M
1.219 + ccr81 = esccr(ymod81,y81,0)
1.220 + print (ccr81)
1.221 +
1.222 +;new eq
1.223 + bias = sum(abs(ymod81-y81)/(ymod81+y81))
1.224 + M81 = (1. - (bias/dimsizes(y81)))*5.
1.225 + print (M81)
1.226 +
1.227 + tRes = True
1.228 + tRes@txFontHeightF = 0.025
1.229 +
1.230 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.231 +
1.232 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.233 +;--------------------------------
1.234 + plot = gsn_csm_xy (wks,y81,ymod81,res) ; create plot
1.235 + frame(wks)
1.236 +
1.237 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.238 + system("rm "+plot_name+"."+plot_type)
1.239 + system("rm "+plot_name+"-1."+plot_type_new)
1.240 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.241 +
1.242 + clear (wks)
1.243 +
1.244 +;(F) scatter model vs ob 933
1.245 +
1.246 + ymod933@long_name = "NPP (gC/m2/year)"
1.247 + y933@long_name = "NPP (gC/m2/year)"
1.248 +;print (ymod33)
1.249 +;print (y933)
1.250 +
1.251 + plot_name = "scatter_model_vs_ob_933"
1.252 + title = model_name +" vs ob 933 sites"
1.253 +
1.254 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.255 + res = True ; plot mods desired
1.256 + res@tiMainString = title ; add title
1.257 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.258 + res@xyMarkers = 16 ; choose type of marker
1.259 + res@xyMarkerColor = "red" ; Marker color
1.260 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.261 + res@tmLabelAutoStride = True ; nice tick mark labels
1.262 +;-------------------------------
1.263 +;compute correlation coef. and M
1.264 + ccr933 = esccr(ymod933,y933,0)
1.265 + print (ccr933)
1.266 +
1.267 +;new eq
1.268 + bias = sum(abs(ymod933-y933)/(ymod933+y933))
1.269 + M933 = (1. - (bias/dimsizes(y933)))*5.
1.270 + print (M933)
1.271 +
1.272 + tRes = True
1.273 + tRes@txFontHeightF = 0.025
1.274 +
1.275 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.276 +
1.277 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.278 +;--------------------------------
1.279 + plot = gsn_csm_xy (wks,y933,ymod933,res) ; create plot
1.280 + frame(wks)
1.281 +
1.282 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.283 + system("rm "+plot_name+"."+plot_type)
1.284 + system("rm "+plot_name+"-1."+plot_type_new)
1.285 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.286 +
1.287 + clear (wks)
1.288 +end