diff -r 0c6405ab2ff4 -r 4be95183fbcd all/01.npp.ncl
--- a/all/01.npp.ncl Mon Jan 26 22:08:20 2009 -0500
+++ b/all/01.npp.ncl Thu Mar 26 14:02:21 2009 -0400
@@ -193,7 +193,8 @@
;(3) global data, interpolated into model grid
dirg = diro + "npp/"
- filg = "Npp_"+model_grid+"_mean.nc"
+ ;filg = "Npp_"+model_grid+"_mean.nc"
+ filg = "npp_"+model_grid+"_mean_2000-2004.nc"
fglobe = addfile (dirg+filg,"r")
nppglobe0 = fglobe->NPP
@@ -225,11 +226,11 @@
;
; rain81 : mm/year
; rainmod : mm/s
-; npp81 : g C/m^2/year
-; nppmod81: g C/m^2/s
-; nppglobe: g C/m^2/year
+; npp81 : gC/m^2/year
+; nppmod81: gC/m^2/s
+; nppglobe: gC/m^2/year
;
-; We want to convert these to "m/year" and "g C/m^2/year".
+; We want to convert these to "m/year" and "gC/m^2/year".
nsec_per_year = 60*60*24*365
@@ -243,27 +244,27 @@
nppmod = nppmod * nsec_per_year
- npp81@units = "gC/m^2/yr"
- nppmod81@units = "gC/m^2/yr"
- npp933@units = "gC/m^2/yr"
- nppmod933@units = "gC/m^2/yr"
- nppmod@units = "gC/m^2/yr"
- nppglobe@units = "gC/m^2/yr"
- rain81@units = "m/yr"
- rainmod81@units = "m/yr"
- rain933@units = "m/yr"
- rainmod933@units = "m/yr"
+ npp81@units = "gC m~S~-2~N~ y~S~-1~N~"
+ nppmod81@units = "gC m~S~-2~N~ y~S~-1~N~"
+ npp933@units = "gC m~S~-2~N~ y~S~-1~N~"
+ nppmod933@units = "gC m~S~-2~N~ y~S~-1~N~"
+ nppmod@units = "gC m~S~-2~N~ y~S~-1~N~"
+ nppglobe@units = "gC m~S~-2~N~ y~S~-1~N~"
+ rain81@units = "m y~S~-1~N~"
+ rainmod81@units = "m y~S~-1~N~"
+ rain933@units = "m y~S~-1~N~"
+ rainmod933@units = "m y~S~-1~N~"
- npp81@long_name = "Obs. NPP (gC/m2/year)"
- npp933@long_name = "Obs. NPP (gC/m2/year)"
- nppmod81@long_name = "Model NPP (gC/m2/year)"
- nppmod933@long_name = "Model NPP (gC/m2/year)"
- nppmod@long_name = "Model NPP (gC/m2/year)"
- nppglobe@long_name = "NPP (gC/m2/year)"
- rain81@long_name = "PREC (m/year)"
- rain933@long_name = "PREC (m/year)"
- rainmod81@long_name = "PREC (m/year)"
- rainmod933@long_name = "PREC (m/year)"
+ npp81@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ npp933@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ nppmod81@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ nppmod933@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ nppglobe@long_name = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ rain81@long_name = "PREC (m y~S~-1~N~)"
+ rain933@long_name = "PREC (m y~S~-1~N~)"
+ rainmod81@long_name = "PREC (m y~S~-1~N~)"
+ rainmod933@long_name = "PREC (m y~S~-1~N~)"
; change longitude from 0-360 to -180-180
x81 = where(x81 .gt. 180., x81 -360., x81 )
@@ -289,8 +290,8 @@
," | Site ID | " \
," Latitude | " \
," Longitude | " \
- ," NPP(gC/m2/year) | " \
- ," PPT(m/year) | " \
+ ," NPP (gC m-2 y-1) | " \
+ ," PPT (m y-1) | " \
,"" \
/)
table_footer = ""
@@ -354,8 +355,8 @@
," Site ID | " \
," Latitude | " \
," Longitude | " \
- ," NPP(gC/m2/year) | " \
- ," PPT(m/year) | " \
+ ," NPP (gC m-2 y-1) | " \
+ ," PPT (m y-1) | " \
,"" \
/)
table_footer = ""
@@ -425,8 +426,8 @@
," Site ID | " \
," Latitude | " \
," Longitude | " \
- ," NPP(gC/m2.year) | " \
- ," RAIN(m/year) | " \
+ ," NPP (gC m-2 y-1) | " \
+ ," RAIN (m y-1) | " \
,"" \
,"" \
," | observed | " \
@@ -506,8 +507,8 @@
," Site ID | " \
," Latitude | " \
," Longitude | " \
- ," NPP(gC/m2.year) | " \
- ," RAIN(m/year) | " \
+ ," NPP (gC m-2 y-1) | " \
+ ," RAIN (m y-1) | " \
,"
" \
,"" \
," | observed | " \
@@ -569,7 +570,8 @@
;***************************************************************************
plot_name = "scatter_model_vs_ob_81"
- title = model_name +" vs Class A observations (81 sites)"
+ ;title = model_name +" vs Class A observations (81 sites)"
+ title = model_name +" (1975-2000) vs Class A observations (81 sites)"
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
res = True ; plot mods desired
@@ -587,7 +589,7 @@
ccr81 = esccr(nppmod81,npp81,0)
;print (ccr81)
- score_max = 2.5
+ score_max = 1.0
bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
M81s = (1. - (bias/dimsizes(y81)))*score_max
@@ -618,15 +620,16 @@
delete (wks)
delete (dum)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;***************************************************************************
;(A)-6 scatter plot, model vs ob 933
;***************************************************************************
plot_name = "scatter_model_vs_ob_933"
- title = model_name +" vs Class B Observations (933 sites)"
+ title = model_name +" (1975-2000) vs Class B Observations (933 sites)"
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
res = True ; plot mods desired
@@ -644,7 +647,7 @@
ccr933 = esccr(nppmod933,npp933,0)
- score_max = 2.5
+ score_max = 1.0
bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
M933s = (1. - (bias/dimsizes(y933)))*score_max
@@ -675,8 +678,9 @@
delete (wks)
delete (dum)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;**************************************************************************
;(B) histogram 81
@@ -711,7 +715,7 @@
ccr81h = esccr(uu,vv,0)
- score_max = 2.5
+ score_max = 2.0
bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
M81h = (1.- (bias/dimsizes(uu)))*score_max
@@ -742,8 +746,8 @@
resh@trYMinF = min(mn_yvalues) - 10.
resh@trYMaxF = max(mx_yvalues) + 10.
- resh@tiYAxisString = "NPP (g C/m2/year)"
- resh@tiXAxisString = "Precipitation (m/year)"
+ resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ resh@tiXAxisString = "Precipitation (m y~S~-1~N~)"
max_bar = new((/2,nx/),graphic)
min_bar = new((/2,nx/),graphic)
@@ -765,6 +769,9 @@
xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
+ ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
+ print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))
+
;-------------------------------
;Attach the vertical bar and the horizontal cap line
@@ -809,15 +816,16 @@
draw(xy)
delete (wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;****************************************************************************
;(B)-2 histogram plot, model vs ob 81 site
;****************************************************************************
plot_name = "histogram_model_vs_ob_81"
- title = model_name+ " vs Class A Observations (81 sites)"
+ title = model_name+ " (1975-2000) vs Class A Observations (81 sites)"
resh@tiMainString = title
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
@@ -847,6 +855,10 @@
gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
+
+ ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
+ print(model_name + ": " + xvalues(1,:) + " " + yvalues(1,:))
+ ;print("All: " + xvalues(:,:) + " " + yvalues(:,:))
;-------------------------------
;Attach the vertical bar and the horizontal cap line
@@ -883,8 +895,9 @@
draw(xy)
delete(wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
delete (RAIN1_1D)
delete (RAIN2_1D)
@@ -935,7 +948,7 @@
ccr933h = esccr(uu,vv,0)
- score_max = 2.5
+ score_max = 2.0
bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
M933h = (1.- (bias/dimsizes(uu)))*score_max
@@ -967,8 +980,8 @@
resh@trYMinF = min(mn_yvalues) - 10.
resh@trYMaxF = max(mx_yvalues) + 10.
- resh@tiYAxisString = "NPP (g C/m2/year)"
- resh@tiXAxisString = "Precipitation (m/year)"
+ resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
+ resh@tiXAxisString = "Precipitation (m y~S~-1~N~)"
max_bar = new((/2,nx/),graphic)
min_bar = new((/2,nx/),graphic)
@@ -990,6 +1003,9 @@
xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
+ ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
+ print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))
+
;-------------------------------
;Attach the vertical bar and the horizontal cap line
@@ -1027,15 +1043,16 @@
delete (xy)
delete (wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;****************************************************************************
;(B)-4 histogram plot, model vs ob 933 site
;****************************************************************************
plot_name = "histogram_model_vs_ob_933"
- title = model_name+ " vs Class B Observations (933 sites)"
+ title = model_name+ " (1975-2000) vs Class B Observations (933 sites)"
resh@tiMainString = title
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
@@ -1065,6 +1082,9 @@
gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
+
+ ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
+ print("Observations: " + xvalues(1,:) + " " + yvalues(1,:))
;-------------------------------
;Attach the vertical bar and the horizontal cap line
@@ -1102,8 +1122,30 @@
delete(xy)
delete(wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
+
+;***************************************************************************
+; Read 2000-2004 dataset for MODIS comparison
+;***************************************************************************
+
+;------------------------------------------------------
+; read model data
+
+ fm = addfile (dirm+film9,"r")
+
+ nppmod0 = fm->NPP
+ xm = fm->lon
+ ym = fm->lat
+
+ delete (fm)
+
+ nppmod = nppmod0(0,:,:)
+ delete (nppmod0)
+ nppmod = nppmod * nsec_per_year
+ nppmod@units = "gC m~S~-2~N~ y~S~-1~N~"
+ nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
;***************************************************************************
;(C) global contour
@@ -1132,7 +1174,7 @@
nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
plot_name = "global_ob"
- title = "Observed MODIS MOD 17"
+ title = "MODIS MOD17A3 (2000-2004)"
resg@tiMainString = title
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
@@ -1142,14 +1184,15 @@
delete (wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;****************************************************************************
;(C)-2 global contour plot, model
;****************************************************************************
plot_name = "global_model"
- title = "Model "+ model_name
+ title = "Model "+ model_name + " (2000-2004)"
resg@tiMainString = title
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
@@ -1159,8 +1202,9 @@
delete (wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;****************************************************************************
;(C)-3 global contour plot, model vs ob
@@ -1190,7 +1234,7 @@
ccrG = esccr(ug,vg,0)
- score_max = 5.0
+ score_max = 2.0
MG = (ccrG*ccrG)* score_max
M_npp_G = sprintf("%.2f", MG)
@@ -1215,14 +1259,14 @@
;--------------------------------------------------------------------
;(a) ob
- title = "Observed MODIS MOD 17"
+ title = "MODIS MOD17A3 (2000-2004)"
resg@tiMainString = title
plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
;(b) model
- title = "Model "+ model_name
+ title = "Model "+ model_name + " (2000-2004)"
resg@tiMainString = title
plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
@@ -1231,7 +1275,7 @@
zz = nppmod
zz = nppmod - nppglobe
- title = "Model_"+model_name+" - Observed"
+ title = "Model "+model_name+" - MODIS MOD17A3"
resg@cnMinLevelValF = -500 ; Min level
resg@cnMaxLevelValF = 500. ; Max level
@@ -1250,8 +1294,9 @@
delete (wks)
delete (plot)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;***************************************************************************
;(D)-1 zonal line plot, ob
@@ -1261,7 +1306,7 @@
vv@long_name = nppglobe@long_name
plot_name = "zonal_ob"
- title = "Observed MODIS MOD 17"
+ title = "MODIS MOD17A3 (2000-2004)"
resz = True
resz@tiMainString = title
@@ -1272,8 +1317,9 @@
delete (wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;****************************************************************************
;(D)-2 zonal line plot, model vs ob
@@ -1288,7 +1334,7 @@
;-------------------------------------------
; compute correlation coef and M score
- score_max = 5.0
+ score_max = 2.0
ccrZ = esccr(s(1,:), s(0,:),0)
MZ = (ccrZ*ccrZ)* score_max
@@ -1303,7 +1349,7 @@
"mv -f "+html_new+" "+html_name)
;-------------------------------------------
plot_name = "zonal_model_vs_ob"
- title = "Zonal Average"
+ title = "Zonal Average (2000-2004)"
resz@tiMainString = title
wks = gsn_open_wks (plot_type,plot_name)
@@ -1324,16 +1370,16 @@
resz@pmLegendDisplayMode = "Always" ; turn on legend
resz@pmLegendSide = "Top" ; Change location of
; resz@pmLegendParallelPosF = .45 ; move units right
- resz@pmLegendParallelPosF = .82 ; move units right
+ resz@pmLegendParallelPosF = .75 ; move units right
resz@pmLegendOrthogonalPosF = -0.4 ; move units down
resz@pmLegendWidthF = 0.10 ; Change width and
resz@pmLegendHeightF = 0.10 ; height of legend.
- resz@lgLabelFontHeightF = .02 ; change font height
+ resz@lgLabelFontHeightF = .015 ; change font height
; resz@lgTitleOn = True ; turn on legend title
; resz@lgTitleString = "Example" ; create legend title
; resz@lgTitleFontHeightF = .025 ; font of legend title
- resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
+ resz@xyExplicitLegendLabels = (/"MODIS MOD17A3",model_name/) ; explicit labels
;--------------------------------------------------------------------
zRes = True
zRes@txFontHeightF = 0.025
@@ -1347,8 +1393,9 @@
delete (wks)
- system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
"rm "+plot_name+"."+plot_type)
+ ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
;***************************************************************************
; add total score and write to file