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