1 ; ***********************************************
2 ; combine 19.metric_plot.ncl and 24.lines.ncl
3 ; ***********************************************
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
7 load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"
8 ;************************************************
10 ;************************************************
11 ; read in data: observed
12 ;************************************************
13 diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
14 fili = "co2_globalView_98.nc"
15 g = addfile (diri+fili,"r")
19 sta = chartostring(g->STATION)
27 ;**************************************************************
28 ; get only the lowest level at each station
29 ;**************************************************************
31 lat_tmp@_FillValue = 1.e+36
34 if (.not. ismissing(lat_tmp(n))) then
35 indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
36 if (dimsizes(indexes) .gt. 1) then
37 lat_tmp(indexes(1:)) = lat_tmp@_FillValue
43 indexes = ind(.not. ismissing(lat_tmp))
44 ;print (dimsizes(indexes))
49 val_ob = val(indexes,:)
50 ;printVarSummary (val_ob)
51 ;print (lat_ob +"/"+lon_ob)
53 ;************************************************
55 ;************************************************
56 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
57 ; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"
58 fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"
60 g = addfile(diri2+fili2,"r")
69 ; get the co2 at the lowest level
72 ; change to unit of observed (u mol/mol)
73 ; Model_units [=] kgCO2 / kgDryAir
74 ; 28.966 = molecular weight of dry air
75 ; 44. = molecular weight of CO2
78 factor = (28.966/44.) * 1e6
83 ; y = where(y0 .lt. 287.,y@_FillValue,y)
85 ; print (min(y)+"/"+max(y))
87 ; interpolate into observed station
88 ; note: model is 0-360E, 90S-90N
90 ; to be able to handle observation at (-89.98,-24.80)
94 i = ind(lon_ob .lt. 0.)
95 lon_ob(i) = lon_ob(i) + 360.
97 yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
99 val_model = yo(pts|:,time|:)
100 val_model_0 = val_model
101 ; printVarSummary (val_model)
102 ; print (min(val_model)+"/"+max(val_model))
105 val_model = val_model - conform(val_model,dim_avg(val_model),0)
106 ; print (min(val_model)+"/"+max(val_model))
109 ;--------------------------------------------------------------
110 ; for metric table plots
112 case_zone = (/"Stations","Amplitude Ratio", \
113 "Correlation Coef","M score","Combined Score"/)
114 nCase_zone = dimsizes(case_zone )
117 var_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)
118 nVar_zone = dimsizes(var_zone)
120 ; arrays to be passed to diagram.
121 case_value_zone = new ((/nCase_zone, nVar_zone/),float )
122 ;--------------------------------------------------------------
123 ; for station line plot
125 ; for x-axis in xyplot
127 mon@long_name = "month"
130 plot_type_new = "png"
132 res = True ; plot mods desired
133 res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker
134 res@xyLineColors = (/"red","black"/) ; change line color
136 ; Add a boxed legend using the more simple method
138 res@pmLegendDisplayMode = "Always"
139 ; res@pmLegendWidthF = 0.1
140 res@pmLegendWidthF = 0.08
141 res@pmLegendHeightF = 0.06
142 ; res@pmLegendOrthogonalPosF = -1.17
143 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
144 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
146 ; res@pmLegendParallelPosF = 0.18
147 res@pmLegendParallelPosF = 0.23 ;(rightward)
149 ; res@lgPerimOn = False
150 res@lgLabelFontHeightF = 0.015
151 res@xyExplicitLegendLabels = (/"b30.061n","observed"/)
152 ;-------------------------------------------------------------------
157 ; maximum score for the zone, 60N-90N
160 ; index of stations in this zone
161 ind_z = ind(lat_ob .ge. 60.)
163 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
164 ; print (val_ob(ind_z,:))
165 ; print (val_model(ind_z,:))
169 ; maximum score for the zone, 30N-60N
172 ; index of stations in this zone
173 ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
175 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
176 ; print (val_ob(ind_z,:))
177 ; print (val_model(ind_z,:))
181 ; maximum score for the zone, EQ-30N
184 ; index of stations in this zone
185 ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
187 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
188 ; print (val_ob(ind_z,:))
189 ; print (val_model(ind_z,:))
193 ; maximum score for the zone, 90S-EQ
196 ; index of stations in this zone
197 ind_z = ind(lat_ob .lt. 0. )
199 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
200 ; print (val_ob(ind_z,:))
201 ; print (val_model(ind_z,:))
204 npts = dimsizes(ind_z)
206 ;-------------------------------------------------------------------------
207 ; for metric table plot
209 amp_ob = new((/npts/),float)
210 amp_model = new((/npts/),float)
212 amp_ratio_sta = new((/npts/),float)
213 ccr_sta = new((/npts/),float)
214 M_sta = new((/npts/),float)
215 score_sta = new((/npts/),float)
216 ;-------------------------------------------------------------------------
217 ; for station line plot
223 plot_data = new((/2,12,npts/),float)
224 plot_data_0 = new((/12,npts/),float)
227 plot_data!1 = "month"
229 plot_data@long_name = "CO2 Seasonal"
231 plot_data_0!0 = "month"
232 plot_data_0!1 = "pts"
233 plot_data_0@long_name = "CO2"
234 ;--------------------------------------------------------------------------
236 amp_ob(n) = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))
237 amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
239 amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
240 ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
241 M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
242 score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
244 print (sta(ind_z(n))+"/"+lat(ind_z(n))+"/"+lon(ind_z(n))+"/"+amp_ratio_sta(n)+"/"+ccr_sta(n)+"/"+M_sta(n)+"/"+score_sta(n))
245 ;----------------------------------------------------------------------
246 ; for station line plot
248 plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
249 plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
251 plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
253 plot_name = sta(ind_z(n))
254 title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"
258 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
259 ;------------------------------------------
262 plot=new(2,graphic) ; create graphic array
263 res@gsnFrame = False ; Do not draw plot
264 res@gsnDraw = False ; Do not advance frame
266 pres = True ; panel plot mods desired
267 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
268 ; indiv. plots in panel
269 pres@gsnMaximize = True ; fill the page
270 ;------------------------------------------
271 res@tiMainString = title ; add title
273 plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res) ; create plot 1
275 plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
277 gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
279 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
280 system("rm "+plot_name+"."+plot_type)
282 ;---------------------------------------------------------------------------
284 ;-------------------------------------------------------------------------
285 ; fo line plot in a zone
287 plot_name = "All_"+npts_str
288 title = plot_name + " in "+ zone
292 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
293 ;-----------------------------------------
296 plot=new(2,graphic) ; create graphic array
297 res@gsnFrame = False ; Do not draw plot
298 res@gsnDraw = False ; Do not advance frame
300 pres = True ; panel plot mods desired
301 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
302 ; indiv. plots in panel
303 pres@gsnMaximize = True ; fill the page
304 ;-----------------------------------------
305 res@tiMainString = title ; add title
307 plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res) ; create plot 1
309 plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
311 gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
313 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
314 system("rm "+plot_name+"."+plot_type)
320 ;---------------------------------------------------------------------------
321 ;---------------------------------------------------------------------------
322 ; metric table in a zone
324 amp_ratio_zone = avg(amp_ratio_sta)
325 ccr_zone = avg(ccr_sta)
326 M_zone = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)
327 score_zone = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
329 print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)
331 case_value_zone(0,z) = npts
332 case_value_zone(1,z) = (/amp_ratio_zone/)
333 case_value_zone(2,z) = (/ccr_zone/)
334 case_value_zone(3,z) = (/M_zone/)
335 case_value_zone(4,z) = (/score_zone/)
337 ; column for station table
338 case_sta = (/"Latitude","Longitude","Amplitude Ratio", \
339 "Correlation Coef","M score","Combined Score"/)
340 nCase_sta = dimsizes(case_sta )
342 ; row for station table
344 nVar_sta = dimsizes(var_sta)
346 ; arrays to be passed to diagram.
347 case_value_sta = new ((/nCase_sta, nVar_sta/),float )
349 case_value_sta(0,:) = (/lat(ind_z)/)
350 case_value_sta(1,:) = (/lon(ind_z)/)
351 case_value_sta(2,:) = (/amp_ratio_sta/)
352 case_value_sta(3,:) = (/ccr_sta/)
353 case_value_sta(4,:) = (/M_sta/)
354 case_value_sta(5,:) = (/score_sta/)
356 ;**************************************************
358 ;**************************************************
360 tt_opt@pltType= "png" ; "eps" [default], "pdf", "ps"
361 ; "png", "gif" [if you have ImageMajik 'convert']
363 tt_opt@tableTitle = "Zone "+ zone
364 plot_name = "table_sta." + zone
366 metrics_table(plot_name, var_sta, case_sta , case_value_sta, tt_opt)
371 delete (amp_ratio_sta)
376 delete (case_value_sta)
378 ;**************************************************
380 ;**************************************************
381 case_value_zone(0,4) = case_value_zone(0,0)+case_value_zone(0,1)+case_value_zone(0,2)+case_value_zone(0,3)
382 case_value_zone(1,4) = 0.
383 case_value_zone(2,4) = 0.
384 case_value_zone(3,4) = 0.
385 case_value_zone(4,4) = case_value_zone(4,0)+case_value_zone(4,1)+case_value_zone(4,2)+case_value_zone(4,3)
388 tt_opt@pltType= "png" ; "eps" [default], "pdf", "ps"
389 ; "png", "gif" [if you have ImageMajik 'convert']
391 tt_opt@tableTitle = "Zone"
392 plot_name = "table_zone"
394 metrics_table(plot_name, var_zone, case_zone , case_value_zone, tt_opt)