forrest@0
|
1 |
; ****************************************************
|
forrest@0
|
2 |
; combine scatter, histogram, global and zonal plots
|
forrest@0
|
3 |
; compute all correlation coef and M score
|
forrest@0
|
4 |
; ****************************************************
|
forrest@0
|
5 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
|
forrest@0
|
6 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
|
forrest@0
|
7 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
|
forrest@0
|
8 |
;************************************************
|
forrest@0
|
9 |
begin
|
forrest@0
|
10 |
|
forrest@0
|
11 |
plot_type = "ps"
|
forrest@0
|
12 |
plot_type_new = "png"
|
forrest@0
|
13 |
;************************************************
|
forrest@0
|
14 |
; read in model data
|
forrest@0
|
15 |
;************************************************
|
forrest@0
|
16 |
;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
|
forrest@0
|
17 |
;model_name = "b30.061n"
|
forrest@0
|
18 |
;model_grid = "T31"
|
forrest@0
|
19 |
|
forrest@0
|
20 |
;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"
|
forrest@0
|
21 |
;model_name = "newcn"
|
forrest@0
|
22 |
;model_grid = "1.9"
|
forrest@0
|
23 |
;--------------------------------------------
|
forrest@0
|
24 |
film = "i01.10cn_1948-2004_ANN_climo.nc"
|
forrest@0
|
25 |
model_name = "10cn"
|
forrest@0
|
26 |
model_grid = "T42"
|
forrest@0
|
27 |
|
forrest@0
|
28 |
dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
|
forrest@0
|
29 |
fm = addfile (dirm+film,"r")
|
forrest@0
|
30 |
|
forrest@0
|
31 |
;unit: gC/m2
|
forrest@0
|
32 |
data1 = fm->LIVESTEMC
|
forrest@0
|
33 |
data2 = fm->DEADSTEMC
|
forrest@0
|
34 |
data3 = fm->LEAFC
|
forrest@0
|
35 |
datamod0 = data1(0,:,:)
|
forrest@0
|
36 |
datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
|
forrest@0
|
37 |
;-------------------------------------------
|
forrest@0
|
38 |
;film = "i01.10casa_1948-2004_ANN_climo.nc"
|
forrest@0
|
39 |
;model_name = "10casa"
|
forrest@0
|
40 |
;model_grid = "T42"
|
forrest@0
|
41 |
|
forrest@0
|
42 |
;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
|
forrest@0
|
43 |
;fm = addfile (dirm+film,"r")
|
forrest@0
|
44 |
|
forrest@0
|
45 |
;unit: gC/m2
|
forrest@0
|
46 |
;factor_WOODC = 0.7
|
forrest@0
|
47 |
;data1 = fm->WOODC
|
forrest@0
|
48 |
;data2 = fm->LEAFC
|
forrest@0
|
49 |
;datamod0 = data1(0,:,:)
|
forrest@0
|
50 |
;datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
|
forrest@0
|
51 |
;----------------------------------------------------
|
forrest@0
|
52 |
html_name = "table.html." + model_name
|
forrest@0
|
53 |
html_new = html_name +".new"
|
forrest@0
|
54 |
system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
|
forrest@0
|
55 |
"mv -f "+html_new+" "+html_name)
|
forrest@0
|
56 |
|
forrest@0
|
57 |
xm = fm->lon
|
forrest@0
|
58 |
ym = fm->lat
|
forrest@0
|
59 |
;************************************************
|
forrest@0
|
60 |
; read in ob data
|
forrest@0
|
61 |
;************************************************
|
forrest@0
|
62 |
ob_name = "LC15_Amazon_Biomass"
|
forrest@0
|
63 |
|
forrest@0
|
64 |
diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
|
forrest@0
|
65 |
filo = "amazon_biomass_"+model_grid+".nc"
|
forrest@0
|
66 |
fo = addfile (diro+filo,"r")
|
forrest@0
|
67 |
|
forrest@0
|
68 |
dataob = fo->BIOMASS
|
forrest@0
|
69 |
xo = fo->lon
|
forrest@0
|
70 |
yo = fo->lat
|
forrest@0
|
71 |
;************************************************
|
forrest@0
|
72 |
; Units for these variables are:
|
forrest@0
|
73 |
; dataob : MgC/ha
|
forrest@0
|
74 |
; datamod0 : gC/m2
|
forrest@0
|
75 |
; We want to convert these to KgC/m2
|
forrest@0
|
76 |
; ha = 100m*100m = 10,000 m2
|
forrest@0
|
77 |
; MgC/ha*1000/10,000 = KgC/m2
|
forrest@0
|
78 |
|
forrest@0
|
79 |
factor_aboveground = 0.5
|
forrest@0
|
80 |
factor_unit_ob = 0.1
|
forrest@0
|
81 |
factor_unit_mod = 0.001
|
forrest@0
|
82 |
|
forrest@0
|
83 |
dataob = dataob * factor_aboveground * factor_unit_ob
|
forrest@0
|
84 |
datamod0 = datamod0 * factor_unit_mod
|
forrest@0
|
85 |
|
forrest@0
|
86 |
dataob@units = "KgC/m2"
|
forrest@0
|
87 |
datamod0@units = "KgC/m2"
|
forrest@0
|
88 |
|
forrest@0
|
89 |
dataob@long_name = "Amazon Biomass"
|
forrest@0
|
90 |
datamod0@long_name = "Amazon Biomass"
|
forrest@0
|
91 |
;********************************************************
|
forrest@0
|
92 |
; get subset of datamod0 that match dataob
|
forrest@0
|
93 |
|
forrest@0
|
94 |
nlon = dimsizes(xo)
|
forrest@0
|
95 |
nlat = dimsizes(yo)
|
forrest@0
|
96 |
|
forrest@0
|
97 |
ind_lonL = ind(xm .eq. xo(0))
|
forrest@0
|
98 |
ind_lonR = ind(xm .eq. xo(nlon-1))
|
forrest@0
|
99 |
ind_latS = ind(ym .eq. yo(0))
|
forrest@0
|
100 |
ind_latN = ind(ym .eq. yo(nlat-1))
|
forrest@0
|
101 |
|
forrest@0
|
102 |
datamod = dataob
|
forrest@0
|
103 |
datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
|
forrest@0
|
104 |
|
forrest@0
|
105 |
; print (datamod)
|
forrest@0
|
106 |
|
forrest@0
|
107 |
;----------------------------------------------------------------------
|
forrest@0
|
108 |
; global contour
|
forrest@0
|
109 |
|
forrest@0
|
110 |
;res
|
forrest@0
|
111 |
resg = True ; Use plot options
|
forrest@0
|
112 |
resg@cnFillOn = True ; Turn on color fill
|
forrest@0
|
113 |
resg@gsnSpreadColors = True ; use full colormap
|
forrest@0
|
114 |
; resg@cnFillMode = "RasterFill" ; Turn on raster color
|
forrest@0
|
115 |
; resg@lbLabelAutoStride = True
|
forrest@0
|
116 |
resg@cnLinesOn = False ; Turn off contourn lines
|
forrest@0
|
117 |
resg@mpFillOn = False ; Turn off map fill
|
forrest@0
|
118 |
resg@gsnAddCyclic = False
|
forrest@0
|
119 |
|
forrest@0
|
120 |
resg@gsnSpreadColors = True ; use full colormap
|
forrest@0
|
121 |
resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
|
forrest@0
|
122 |
resg@cnMinLevelValF = 0. ; Min level
|
forrest@0
|
123 |
resg@cnMaxLevelValF = 30. ; Max level
|
forrest@0
|
124 |
resg@cnLevelSpacingF = 2. ; interval
|
forrest@0
|
125 |
|
forrest@0
|
126 |
resg@mpMinLatF = -21.1 ; range to zoom in on
|
forrest@0
|
127 |
resg@mpMaxLatF = 13.8
|
forrest@0
|
128 |
resg@mpMinLonF = 277.28
|
forrest@0
|
129 |
resg@mpMaxLonF = 326.43
|
forrest@0
|
130 |
;------------------------------------------------------------------------
|
forrest@0
|
131 |
;global contour ob
|
forrest@0
|
132 |
|
forrest@0
|
133 |
plot_name = "global_ob"
|
forrest@0
|
134 |
title = ob_name+" "+ model_grid
|
forrest@0
|
135 |
resg@tiMainString = title
|
forrest@0
|
136 |
|
forrest@0
|
137 |
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
|
forrest@0
|
138 |
gsn_define_colormap(wks,"gui_default") ; choose colormap
|
forrest@0
|
139 |
|
forrest@0
|
140 |
plot = gsn_csm_contour_map_ce(wks,dataob,resg)
|
forrest@0
|
141 |
frame(wks)
|
forrest@0
|
142 |
|
forrest@0
|
143 |
system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
|
forrest@0
|
144 |
"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
|
forrest@0
|
145 |
"rm "+plot_name+"-*."+plot_type_new+";"+ \
|
forrest@0
|
146 |
"rm "+plot_name+"."+plot_type)
|
forrest@0
|
147 |
|
forrest@0
|
148 |
clear (wks)
|
forrest@0
|
149 |
;------------------------------------------------------------------------
|
forrest@0
|
150 |
;global contour model
|
forrest@0
|
151 |
|
forrest@0
|
152 |
plot_name = "global_model"
|
forrest@0
|
153 |
title = "Model "+ model_name
|
forrest@0
|
154 |
resg@tiMainString = title
|
forrest@0
|
155 |
|
forrest@0
|
156 |
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
|
forrest@0
|
157 |
gsn_define_colormap(wks,"gui_default") ; choose colormap
|
forrest@0
|
158 |
|
forrest@0
|
159 |
plot = gsn_csm_contour_map_ce(wks,datamod,resg)
|
forrest@0
|
160 |
frame(wks)
|
forrest@0
|
161 |
|
forrest@0
|
162 |
system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
|
forrest@0
|
163 |
"mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
|
forrest@0
|
164 |
"rm "+plot_name+"-*."+plot_type_new+";"+ \
|
forrest@0
|
165 |
"rm "+plot_name+"."+plot_type)
|
forrest@0
|
166 |
|
forrest@0
|
167 |
clear (wks)
|
forrest@0
|
168 |
;------------------------------------------------------------------------
|
forrest@0
|
169 |
;global contour model vs ob
|
forrest@0
|
170 |
|
forrest@0
|
171 |
plot_name = "global_model_vs_ob"
|
forrest@0
|
172 |
|
forrest@0
|
173 |
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
|
forrest@0
|
174 |
gsn_define_colormap(wks,"gui_default") ; choose colormap
|
forrest@0
|
175 |
|
forrest@0
|
176 |
delete (plot)
|
forrest@0
|
177 |
plot=new(3,graphic) ; create graphic array
|
forrest@0
|
178 |
|
forrest@0
|
179 |
resg@gsnFrame = False ; Do not draw plot
|
forrest@0
|
180 |
resg@gsnDraw = False ; Do not advance frame
|
forrest@0
|
181 |
|
forrest@0
|
182 |
;(d) compute correlation coef and M score
|
forrest@0
|
183 |
|
forrest@0
|
184 |
uu = ndtooned(datamod)
|
forrest@0
|
185 |
vv = ndtooned(dataob)
|
forrest@0
|
186 |
|
forrest@0
|
187 |
good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
|
forrest@0
|
188 |
|
forrest@0
|
189 |
ug = uu(good)
|
forrest@0
|
190 |
vg = vv(good)
|
forrest@0
|
191 |
|
forrest@0
|
192 |
ccrG = esccr(ug,vg,0)
|
forrest@0
|
193 |
; print (ccrG)
|
forrest@0
|
194 |
|
forrest@0
|
195 |
MG = (ccrG*ccrG)* 5.0
|
forrest@0
|
196 |
; new eq
|
forrest@0
|
197 |
bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
|
forrest@0
|
198 |
MG = (1. - (bias/dimsizes(ug)))*5.
|
forrest@0
|
199 |
|
forrest@0
|
200 |
M_biomass = sprintf("%.2f", MG)
|
forrest@0
|
201 |
system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \
|
forrest@0
|
202 |
"mv -f "+html_new+" "+html_name)
|
forrest@0
|
203 |
print (M_biomass)
|
forrest@0
|
204 |
|
forrest@0
|
205 |
; plot correlation coef
|
forrest@0
|
206 |
|
forrest@0
|
207 |
gRes = True
|
forrest@0
|
208 |
gRes@txFontHeightF = 0.02
|
forrest@0
|
209 |
gRes@txAngleF = 90
|
forrest@0
|
210 |
|
forrest@0
|
211 |
correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
|
forrest@0
|
212 |
|
forrest@0
|
213 |
gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
|
forrest@0
|
214 |
;--------------------------------------------------------------------
|
forrest@0
|
215 |
|
forrest@0
|
216 |
;(a) ob
|
forrest@0
|
217 |
|
forrest@0
|
218 |
title = ob_name+" "+ model_grid
|
forrest@0
|
219 |
resg@tiMainString = title
|
forrest@0
|
220 |
|
forrest@0
|
221 |
plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)
|
forrest@0
|
222 |
|
forrest@0
|
223 |
;(b) model
|
forrest@0
|
224 |
|
forrest@0
|
225 |
title = "Model "+ model_name
|
forrest@0
|
226 |
resg@tiMainString = title
|
forrest@0
|
227 |
|
forrest@0
|
228 |
plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
|
forrest@0
|
229 |
|
forrest@0
|
230 |
;(c) model-ob
|
forrest@0
|
231 |
|
forrest@0
|
232 |
zz = datamod
|
forrest@0
|
233 |
zz = datamod - dataob
|
forrest@0
|
234 |
title = "Model_"+model_name+" - Observed"
|
forrest@0
|
235 |
|
forrest@0
|
236 |
resg@cnMinLevelValF = -10. ; Min level
|
forrest@0
|
237 |
resg@cnMaxLevelValF = 10. ; Max level
|
forrest@0
|
238 |
resg@cnLevelSpacingF = 2. ; interval
|
forrest@0
|
239 |
resg@tiMainString = title
|
forrest@0
|
240 |
|
forrest@0
|
241 |
plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
|
forrest@0
|
242 |
|
forrest@0
|
243 |
pres = True ; panel plot mods desired
|
forrest@0
|
244 |
; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
|
forrest@0
|
245 |
; indiv. plots in panel
|
forrest@0
|
246 |
pres@gsnMaximize = True ; fill the page
|
forrest@0
|
247 |
|
forrest@0
|
248 |
gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
|
forrest@0
|
249 |
|
forrest@0
|
250 |
system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
|
forrest@0
|
251 |
"rm "+plot_name+"."+plot_type)
|
forrest@0
|
252 |
|
forrest@0
|
253 |
clear (wks)
|
forrest@0
|
254 |
delete (plot)
|
forrest@0
|
255 |
;-------------------------------------------------------------------
|
forrest@0
|
256 |
temp_name = "biomass." + model_name
|
forrest@0
|
257 |
system("mkdir -p " + temp_name+";"+ \
|
forrest@0
|
258 |
"cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
|
forrest@0
|
259 |
"mv *.png " + temp_name +";"+ \
|
forrest@0
|
260 |
"tar cf "+ temp_name +".tar " + temp_name)
|
forrest@0
|
261 |
;-------------------------------------------------------------------
|
forrest@0
|
262 |
end
|