Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ; ***********************************************
2 ; using gsn_table for all
3 ; combine 19.metric_plot.ncl and 24.lines.ncl
4 ; ***********************************************
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
8 ;load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"
9 ;************************************************
11 ;************************************************
12 ; read in data: observed
13 ;************************************************
14 diri = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
15 fili = "co2_globalView_98.nc"
16 g = addfile (diri+fili,"r")
20 sta = chartostring(g->STATION)
28 ;**************************************************************
29 ; get only the lowest level at each station
30 ;**************************************************************
32 lat_tmp@_FillValue = 1.e+36
35 if (.not. ismissing(lat_tmp(n))) then
36 indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
37 if (dimsizes(indexes) .gt. 1) then
38 lat_tmp(indexes(1:)) = lat_tmp@_FillValue
44 indexes = ind(.not. ismissing(lat_tmp))
45 ;print (dimsizes(indexes))
50 val_ob = val(indexes,:)
51 ;printVarSummary (val_ob)
52 ;print (lat_ob +"/"+lon_ob)
54 ;************************************************
56 ;************************************************
57 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
58 ; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"
59 fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"
61 g = addfile(diri2+fili2,"r")
70 ; get the co2 at the lowest level
73 ; change to unit of observed (u mol/mol)
74 ; Model_units [=] kgCO2 / kgDryAir
75 ; 28.966 = molecular weight of dry air
76 ; 44. = molecular weight of CO2
79 factor = (28.966/44.) * 1e6
84 ; y = where(y0 .lt. 287.,y@_FillValue,y)
86 ; print (min(y)+"/"+max(y))
88 ; interpolate into observed station
89 ; note: model is 0-360E, 90S-90N
91 ; to be able to handle observation at (-89.98,-24.80)
95 i = ind(lon_ob .lt. 0.)
96 lon_ob(i) = lon_ob(i) + 360.
98 yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
100 val_model = yo(pts|:,time|:)
101 val_model_0 = val_model
102 ; printVarSummary (val_model)
103 ; print (min(val_model)+"/"+max(val_model))
106 val_model = val_model - conform(val_model,dim_avg(val_model),0)
107 ; print (min(val_model)+"/"+max(val_model))
111 ;*******************************************************************
113 ;*******************************************************************
115 table_header_name = "Zone"
117 ; column (not including header column)
118 col_header_zone = (/"Stations","Amplitude Ratio", \
119 "Correlation Coef","M score","Combined Score"/)
120 ncol_zone = dimsizes(col_header_zone )
122 ; row (not including header row)
123 row_header_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)
124 nrow_zone = dimsizes(row_header_zone)
126 ; arrays to be passed to table.
127 value_zone = new ((/nrow_zone, ncol_zone/),string )
128 ;--------------------------------------------------------------------
130 table_length_zone = 0.4
133 ncr1 = (/1,1/) ; 1 row, 1 column
134 x1 = (/0.005,0.15/) ; Start and end X
135 y1 = (/0.900,0.995/) ; Start and end Y
136 text1 = table_header_name
138 res1@txFontHeightF = 0.03
139 res1@gsFillColor = "CornFlowerBlue"
141 ; Column header (equally space in x2)
142 ncr2 = (/1,ncol_zone/) ; 1 rows, 5 columns
143 x2 = (/x1(1),0.995/) ; start from end of x1
145 text2 = col_header_zone
147 res2@txFontHeightF = 0.015
148 res2@gsFillColor = "Gray"
150 ; Row header (equally space in y2)
151 ncr3 = (/nrow_zone,1/) ; 5 rows, 1 columns
153 y3 = (/1.0-table_length_zone,0.900/) ; end at start of y1
154 text3 = row_header_zone
156 res3@txFontHeightF = 0.02
157 res3@gsFillColor = "Gray"
160 ncr4 = (/nrow_zone,ncol_zone/) ; 5 rows, 5 columns
161 x4 = x2 ; Start and end x
162 y4 = y3 ; Start and end Y
163 text4 = new((/nrow_zone,ncol_zone/),string)
165 color_fill4 = new((/nrow_zone,ncol_zone/),string)
166 color_fill4 = "white"
167 color_fill4(:,ncol_zone-1) = "grey"
169 res4 = True ; Set up resource list
170 ; res4@gsnDebug = True ; Useful to print NDC row,col values used.
171 res4@txFontHeightF = 0.02
172 res4@gsFillColor = color_fill4
176 ;*******************************************************************
177 ; for table -- station
178 ;*******************************************************************
179 ; column (not including header column)
180 col_header_sta = (/"Latitude","Longitude","Amplitude Ratio", \
181 "Correlation Coef","M score","Combined Score"/)
182 ncol_sta = dimsizes(col_header_sta )
183 ;-------------------------------------------------------------------
186 ncr5 = (/1,1/) ; 1 row, 1 column
187 x5 = (/0.005,0.15/) ; Start and end X
188 y5 = (/0.900,0.995/) ; Start and end Y
190 res5@txFontHeightF = 0.02
191 res5@gsFillColor = "CornFlowerBlue"
193 ; Column header (equally space in x2)
194 ncr6 = (/1,ncol_sta/) ; 1 rows, 5 columns
195 x6 = (/x5(1),0.995/) ; start from end of x1
197 text6 = col_header_sta
199 res6@txFontHeightF = 0.012
200 res6@gsFillColor = "Gray"
202 ; Row header (equally space in y2)
205 res7@txFontHeightF = 0.015
206 res7@gsFillColor = "Gray"
207 ;--------------------------------------------------------------
208 ; for station line plot
210 ; for x-axis in xyplot
212 mon@long_name = "month"
215 plot_type_new = "png"
217 res = True ; plot mods desired
218 res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker
219 res@xyLineColors = (/"red","black"/) ; change line color
221 ; Add a boxed legend using the more simple method
223 res@pmLegendDisplayMode = "Always"
224 ; res@pmLegendWidthF = 0.1
225 res@pmLegendWidthF = 0.08
226 res@pmLegendHeightF = 0.06
227 ; res@pmLegendOrthogonalPosF = -1.17
228 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
229 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
231 ; res@pmLegendParallelPosF = 0.18
232 res@pmLegendParallelPosF = 0.23 ;(rightward)
234 ; res@lgPerimOn = False
235 res@lgLabelFontHeightF = 0.015
236 res@xyExplicitLegendLabels = (/"b30.061n","observed"/)
237 ;-------------------------------------------------------------------
242 ; maximum score for the zone, 60N-90N
245 ; index of stations in this zone
246 ind_z = ind(lat_ob .ge. 60.)
248 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
249 ; print (val_ob(ind_z,:))
250 ; print (val_model(ind_z,:))
254 ; maximum score for the zone, 30N-60N
257 ; index of stations in this zone
258 ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
260 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
261 ; print (val_ob(ind_z,:))
262 ; print (val_model(ind_z,:))
266 ; maximum score for the zone, EQ-30N
269 ; index of stations in this zone
270 ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
272 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
273 ; print (val_ob(ind_z,:))
274 ; print (val_model(ind_z,:))
278 ; maximum score for the zone, 90S-EQ
281 ; index of stations in this zone
282 ind_z = ind(lat_ob .lt. 0. )
284 ; print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
285 ; print (val_ob(ind_z,:))
286 ; print (val_model(ind_z,:))
289 npts = dimsizes(ind_z)
291 ;-------------------------------------------------------------------------
292 ; for metric table computation
294 amp_ob = new((/npts/),float)
295 amp_model = new((/npts/),float)
297 amp_ratio_sta = new((/npts/),float)
298 ccr_sta = new((/npts/),float)
299 M_sta = new((/npts/),float)
300 score_sta = new((/npts/),float)
302 ;----------------------
303 ; for table -- station
304 ;----------------------
306 ; row (not including header row)
310 ncr5 = (/1,1/) ; 1 row, 1 column
311 x5 = (/0.005,0.15/) ; Start and end X
312 y5 = (/0.900,0.995/) ; Start and end Y
314 res5@txFontHeightF = 0.02
315 res5@gsFillColor = "CornFlowerBlue"
317 ; Column header (equally space in x2)
318 ncr6 = (/1,ncol_sta/) ; 1 rows, 5 columns
319 x6 = (/x1(1),0.995/) ; start from end of x1
321 text6 = col_header_sta
323 res6@txFontHeightF = 0.012
324 res6@gsFillColor = "Gray"
326 ; Row header (equally space in y2)
327 ncr7 = (/nrow_sta,1/) ; 5 rows, 1 columns
329 ; y7 = (/1.0-table_length_sta,0.900/) ; end at start of y1
330 text7 = new((/nrow_sta/),string)
332 res7@txFontHeightF = 0.015
333 res7@gsFillColor = "Gray"
334 ;-------------------------------------------------------------------------
335 ; for station line plot
341 plot_data = new((/2,12,npts/),float)
342 plot_data_0 = new((/12,npts/),float)
345 plot_data!1 = "month"
347 plot_data@long_name = "CO2 Seasonal"
349 plot_data_0!0 = "month"
350 plot_data_0!1 = "pts"
351 plot_data_0@long_name = "CO2"
352 ;--------------------------------------------------------------------------
354 amp_ob(n) = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))
355 amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
357 amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
358 ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
359 M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
360 score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
362 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))
363 ;----------------------------------------------------------------------
364 ; for station line plot
366 plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
367 plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
369 plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
371 plot_name = sta(ind_z(n))
372 title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"
376 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
377 ;------------------------------------------
380 plot=new(2,graphic) ; create graphic array
381 res@gsnFrame = False ; Do not draw plot
382 res@gsnDraw = False ; Do not advance frame
384 pres = True ; panel plot mods desired
385 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
386 ; indiv. plots in panel
387 pres@gsnMaximize = True ; fill the page
388 ;------------------------------------------
389 res@tiMainString = title ; add title
391 plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res) ; create plot 1
393 plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
395 gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
397 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
398 system("rm "+plot_name+"."+plot_type)
400 ;---------------------------------------------------------------------------
402 ;-------------------------------------------------------------------------
403 ; for line plot in a zone
405 plot_name = "All_"+npts_str
406 title = plot_name + " in "+ zone
410 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
411 ;-----------------------------------------
414 plot=new(2,graphic) ; create graphic array
415 res@gsnFrame = False ; Do not draw plot
416 res@gsnDraw = False ; Do not advance frame
418 pres = True ; panel plot mods desired
419 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
420 ; indiv. plots in panel
421 pres@gsnMaximize = True ; fill the page
422 ;-----------------------------------------
423 res@tiMainString = title ; add title
425 plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res) ; create plot 1
427 plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
429 gsn_panel(wks,plot,(/2,1/),pres) ; create panel plot
431 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
432 system("rm "+plot_name+"."+plot_type)
438 ;---------------------------------------------------------------------------
439 ; for zone table value
441 amp_ratio_zone = avg(amp_ratio_sta)
442 ccr_zone = avg(ccr_sta)
443 M_zone = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)
444 score_zone = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
446 print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)
448 text4(z,0) = sprintf("%5.2f", npts)
449 text4(z,1) = sprintf("%5.2f", amp_ratio_zone)
450 text4(z,2) = sprintf("%5.2f", ccr_zone)
451 text4(z,3) = sprintf("%5.2f", M_zone)
452 text4(z,4) = sprintf("%5.2f", score_zone)
453 ;---------------------------------------------------------------------------
455 ;----------------------------------
456 ; header value for station table
459 ; row value for station table
463 table_length_sta = 0.5
466 table_length_sta = 0.995
469 table_length_sta = 0.8
472 table_length_sta = 0.8
476 y7 = (/1.0-table_length_sta,0.900/) ; end at start of y1
479 ncr8 = (/nrow_sta,ncol_sta/) ; 5 rows, 5 columns
480 x8 = x6 ; Start and end x
481 y8 = y7 ; Start and end Y
482 text8 = new((/nrow_sta,ncol_sta/),string)
485 color_fill8 = "white"
486 color_fill8(:,ncol_sta-1) = "grey"
488 res8 = True ; Set up resource list
489 res8@gsnDebug = True ; Useful to print NDC row,col values used.
490 res8@txFontHeightF = 0.02
491 res8@gsFillColor = color_fill8
495 ; table value for station table
496 text8(:,0) = sprintf("%5.2f", (/lat(ind_z)/))
497 text8(:,1) = sprintf("%5.2f", (/lon(ind_z)/))
498 text8(:,2) = sprintf("%5.2f", (/amp_ratio_sta/))
499 text8(:,3) = sprintf("%5.2f", (/ccr_sta/))
500 text8(:,4) = sprintf("%5.2f", (/M_sta/))
501 text8(:,5) = sprintf("%5.2f", (/score_sta/))
503 plot_name = "table_sta." + zone
505 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
507 gsn_table(wks,ncr5,x5,y5,text5,res5)
508 gsn_table(wks,ncr6,x6,y6,text6,res6)
509 gsn_table(wks,ncr7,x7,y7,text7,res7)
510 gsn_table(wks,ncr8,x8,y8,text8,res8)
514 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
515 system("rm "+plot_name+"."+plot_type)
520 delete (amp_ratio_sta)
530 ;**************************************************
532 ;**************************************************
534 text4(4,0) = sum(stringtofloat(text4(0:3,0)))
538 text4(4,4) = sum(stringtofloat(text4(0:3,4)))
540 plot_name = "table_zone"
541 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
543 gsn_table(wks,ncr1,x1,y1,text1,res1)
544 gsn_table(wks,ncr2,x2,y2,text2,res2)
545 gsn_table(wks,ncr3,x3,y3,text3,res3)
546 gsn_table(wks,ncr4,x4,y4,text4,res4)
550 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
551 system("rm "+plot_name+"."+plot_type)