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