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