forrest@0: ; *********************************************** forrest@0: ; combine all scatter 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: begin forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;************************************************ forrest@0: ; read in ob data forrest@0: ;************************************************ forrest@0: forrest@0: ;(A) plot ob scatter_81 forrest@0: forrest@0: dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" forrest@0: fil81 = "data.81.nc" forrest@0: f81 = addfile (dir81+fil81,"r") forrest@0: forrest@0: x81 = f81->SITE_ID forrest@0: y81 = f81->TNPP_C forrest@0: xo81 = f81->LONG_DD forrest@0: yo81 = f81->LAT_DD forrest@0: forrest@0: x81@long_name = "SITE_ID" forrest@0: y81@long_name = "NPP (gC/m2/year)" forrest@0: forrest@0: plot_name = "scatter_ob_81" forrest@0: title = "Observed 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: plot = gsn_csm_xy (wks,x81,y81,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(B) plot ob scatter_933 forrest@0: forrest@0: dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" forrest@0: fil933 = "data.933.nc" forrest@0: f933 = addfile (dir933+fil933,"r") forrest@0: forrest@0: x933 = f933->SITE_ID forrest@0: y933 = f933->TNPP_C forrest@0: xo933 = f933->LONG_DD forrest@0: yo933 = f933->LAT_DD forrest@0: forrest@0: x933@long_name = "SITE_ID" forrest@0: y933@long_name = "NPP (gC/m2/year)" forrest@0: forrest@0: plot_name = "scatter_ob_933" forrest@0: title = "Observed 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: plot = gsn_csm_xy (wks,x933,y933,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;************************************************ forrest@0: ; read in model data forrest@0: ;************************************************ forrest@0: forrest@0: ;(C) plot model scatter_81 forrest@0: forrest@0: model_name = "b30.061n" forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = "b30.061n_1995-2004_ANN_climo_lnd.nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: ymod = fm->NPP forrest@0: xm = fm->lon forrest@0: ym = fm->lat forrest@0: forrest@0: nx = dimsizes(xo81) forrest@0: do i= 0,nx-1 forrest@0: if (xo81(i) .lt. 0.) then forrest@0: xo81(i) = xo81(i)+ 360. forrest@0: end if forrest@0: end do forrest@0: print (xo81) forrest@0: forrest@0: sec_to_year = 86400.*365. forrest@0: forrest@0: ymod81 = linint2_points(xm,ym,ymod(0,:,:),True,xo81,yo81,0) * sec_to_year forrest@0: xmod81 = x81 forrest@0: forrest@0: ymod81@long_name = "NPP (gC/m2/year)" forrest@0: xmod81@long_name = "SITE_ID" forrest@0: print (ymod81) forrest@0: print (xmod81) forrest@0: forrest@0: plot_name = "scatter_model_81" forrest@0: title = "Model "+ model_name +" 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: plot = gsn_csm_xy (wks,xmod81,ymod81,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(D) plot model scatter_933 forrest@0: forrest@0: model_name = "b30.061n" forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = "b30.061n_1995-2004_ANN_climo_lnd.nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: ymod = fm->NPP forrest@0: xm = fm->lon forrest@0: ym = fm->lat forrest@0: forrest@0: nx = dimsizes(xo933) forrest@0: do i= 0,nx-1 forrest@0: if (xo933(i) .lt. 0.) then forrest@0: xo933(i) = xo933(i)+ 360. forrest@0: end if forrest@0: end do forrest@0: print (xo933) forrest@0: forrest@0: sec_to_year = 86400.*365. forrest@0: forrest@0: ymod933 = linint2_points(xm,ym,ymod(0,:,:),True,xo933,yo933,0) * sec_to_year forrest@0: xmod933 = x933 forrest@0: forrest@0: ymod933@long_name = "NPP (gC/m2/year)" forrest@0: xmod933@long_name = "SITE_ID" forrest@0: ;print (ymod933) forrest@0: ;print (xmod933) forrest@0: forrest@0: plot_name = "scatter_model_933" forrest@0: title = "Model "+ model_name +" 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: plot = gsn_csm_xy (wks,xmod933,ymod933,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(E) scatter model vs ob 81 forrest@0: forrest@0: ymod81@long_name = "NPP (gC/m2/year)" forrest@0: y81@long_name = "NPP (gC/m2/year)" forrest@0: ;print (ymod81) forrest@0: ;print (y81) forrest@0: forrest@0: plot_name = "scatter_model_vs_ob_81" forrest@0: title = model_name +" vs ob 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: ;compute correlation coef. and M forrest@0: ccr81 = esccr(ymod81,y81,0) forrest@0: print (ccr81) forrest@0: forrest@0: ;new eq forrest@0: bias = sum(abs(ymod81-y81)/(ymod81+y81)) forrest@0: M81 = (1. - (bias/dimsizes(y81)))*5. forrest@0: print (M81) 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,y81,ymod81,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(F) scatter model vs ob 933 forrest@0: forrest@0: ymod933@long_name = "NPP (gC/m2/year)" forrest@0: y933@long_name = "NPP (gC/m2/year)" forrest@0: ;print (ymod33) forrest@0: ;print (y933) forrest@0: forrest@0: plot_name = "scatter_model_vs_ob_933" forrest@0: title = model_name +" vs ob 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: ;compute correlation coef. and M forrest@0: ccr933 = esccr(ymod933,y933,0) forrest@0: print (ccr933) forrest@0: forrest@0: ;new eq forrest@0: bias = sum(abs(ymod933-y933)/(ymod933+y933)) forrest@0: M933 = (1. - (bias/dimsizes(y933)))*5. forrest@0: print (M933) 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,y933,ymod933,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: end