1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/99.all.ncl.0 Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1683 @@
1.4 +; ****************************************************
1.5 +; combine scatter, histogram, global and zonal plots
1.6 +; compute all correlation coef and M score
1.7 +; add 1-to-1 line in scatter plots
1.8 +; ****************************************************
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.12 +;************************************************
1.13 +procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
1.14 + ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
1.15 + ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
1.16 + ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
1.17 + ,dx4[1]:numeric \
1.18 + )
1.19 +begin
1.20 +; Calculaee "nice" bins for binning the data in equally spaced ranges.
1.21 +; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
1.22 +; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
1.23 +
1.24 + nbins = 15 ; Number of bins to use.
1.25 +
1.26 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.27 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.28 + range = fspan(nicevals(0),nicevals(1),nvals)
1.29 +
1.30 +; Use this range information to grab all the values in a
1.31 +; particular range, and then take an average.
1.32 +
1.33 + nr = dimsizes(range)
1.34 + nx = nr-1
1.35 +; print (nx)
1.36 +
1.37 +; xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.38 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.39 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.40 + dx4 = dx/4 ; 1/4 of the range
1.41 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.42 +; yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.43 +; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.44 +; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.45 +
1.46 + do nd=0,1
1.47 +
1.48 +; See if we are doing model or observational data.
1.49 +
1.50 + if(nd.eq.0) then
1.51 + data = RAIN1_1D
1.52 + npp_data = NPP1_1D
1.53 + else
1.54 + data = RAIN2_1D
1.55 + npp_data = NPP2_1D
1.56 + end if
1.57 +
1.58 +; Loop through each range and check for values.
1.59 +
1.60 + do i=0,nr-2
1.61 + if (i.ne.(nr-2)) then
1.62 +; print("")
1.63 +; print("In range ["+range(i)+","+range(i+1)+")")
1.64 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.65 + else
1.66 +; print("")
1.67 +; print("In range ["+range(i)+",)")
1.68 + idx = ind(range(i).le.data)
1.69 + end if
1.70 +
1.71 +; Calculate average, and get min and max.
1.72 +
1.73 + if(.not.any(ismissing(idx))) then
1.74 + yvalues(nd,i) = avg(npp_data(idx))
1.75 + mn_yvalues(nd,i) = min(npp_data(idx))
1.76 + mx_yvalues(nd,i) = max(npp_data(idx))
1.77 + count = dimsizes(idx)
1.78 + else
1.79 + count = 0
1.80 + yvalues(nd,i) = yvalues@_FillValue
1.81 + mn_yvalues(nd,i) = yvalues@_FillValue
1.82 + mx_yvalues(nd,i) = yvalues@_FillValue
1.83 + end if
1.84 +
1.85 +; Print out information.
1.86 +
1.87 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.88 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.89 +
1.90 +
1.91 +; Clean up for next time in loop.
1.92 +
1.93 + delete(idx)
1.94 + end do
1.95 + delete(data)
1.96 + delete(npp_data)
1.97 + end do
1.98 +end
1.99 +
1.100 +; Main code.
1.101 +begin
1.102 +
1.103 + plot_type = "ps"
1.104 + plot_type_new = "png"
1.105 +;************************************************
1.106 +; read in model data
1.107 +;************************************************
1.108 +;film = "i01.06cn_1798-2004_ANN_climo.nc"
1.109 +;model_name = "06cn"
1.110 +;model_grid = "T42"
1.111 +
1.112 +;film = "i01.06casa_1798-2004_ANN_climo.nc"
1.113 +;model_name = "06casa"
1.114 +;model_grid = "T42"
1.115 +
1.116 + film = "i01.10casa_1948-2004_ANN_climo.nc"
1.117 + model_name = "10casa"
1.118 + model_grid = "T42"
1.119 +
1.120 +;film = "i01.10cn_1948-2004_ANN_climo.nc"
1.121 +;model_name = "10cn"
1.122 +;model_grid = "T42"
1.123 +;--------------------------------------------------
1.124 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.125 + fm = addfile (dirm+film,"r")
1.126 +
1.127 + nppmod0 = fm->NPP
1.128 + rainmod0 = fm->RAIN
1.129 + xm = fm->lon
1.130 + ym = fm->lat
1.131 +
1.132 +;************************************************
1.133 +; read in ob data
1.134 +;************************************************
1.135 +;(1)
1.136 + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.137 + fil81 = "data.81.nc"
1.138 + f81 = addfile (dir81+fil81,"r")
1.139 +
1.140 + id81 = f81->SITE_ID
1.141 + npp81 = f81->TNPP_C
1.142 + rain81 = tofloat(f81->PREC_ANN)
1.143 + x81 = f81->LONG_DD
1.144 + y81 = f81->LAT_DD
1.145 +
1.146 + id81@long_name = "SITE_ID"
1.147 +
1.148 +;change longitude from (-180,180) to (0,360)
1.149 +;for model data interpolation
1.150 +
1.151 + do i= 0,dimsizes(x81)-1
1.152 + if (x81(i) .lt. 0.) then
1.153 + x81(i) = x81(i)+ 360.
1.154 + end if
1.155 + end do
1.156 +;print (x81)
1.157 +;-------------------------------------
1.158 +;(2)
1.159 + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.160 + fil933 = "data.933.nc"
1.161 + f933 = addfile (dir933+fil933,"r")
1.162 +
1.163 + id933 = f933->SITE_ID
1.164 + npp933 = f933->TNPP_C
1.165 + rain933 = f933->PREC
1.166 + x933 = f933->LONG_DD
1.167 + y933 = f933->LAT_DD
1.168 +
1.169 + id933@long_name = "SITE_ID"
1.170 +
1.171 +;change longitude from (-180,180) to (0,360)
1.172 +;for model data interpolation
1.173 +
1.174 + do i= 0,dimsizes(x933)-1
1.175 + if (x933(i) .lt. 0.) then
1.176 + x933(i) = x933(i)+ 360.
1.177 + end if
1.178 + end do
1.179 +;print (x933)
1.180 +;----------------------------------------
1.181 +;(3)
1.182 + dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.183 + filglobe = "Npp_"+model_grid+"_mean.nc"
1.184 + fglobe = addfile (dirglobe+filglobe,"r")
1.185 +
1.186 + nppglobe0 = fglobe->NPP
1.187 + nppglobe = nppglobe0
1.188 +;***********************************************************************
1.189 +; interpolate model data
1.190 +;***********************************************************************
1.191 + nppmod = nppmod0(0,:,:)
1.192 + rainmod = rainmod0(0,:,:)
1.193 + delete (nppmod0)
1.194 + delete (rainmod0)
1.195 +
1.196 + nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
1.197 +
1.198 + nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
1.199 +
1.200 + rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
1.201 +
1.202 + rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
1.203 +
1.204 +;************************************************
1.205 +; Units for these variables are:
1.206 +;
1.207 +; rain81 : mm/year
1.208 +; rainmod : mm/s
1.209 +; npp81 : g C/m^2/year
1.210 +; nppmod81: g C/m^2/s
1.211 +; nppglobe: g C/m^2/year
1.212 +;
1.213 +; We want to convert these to "m/year" and "g C/m^2/year".
1.214 +;
1.215 + nsec_per_year = 60*60*24*365
1.216 +
1.217 + rain81 = rain81 / 1000.
1.218 + rainmod81 = (rainmod81/ 1000.) * nsec_per_year
1.219 + nppmod81 = nppmod81 * nsec_per_year
1.220 +
1.221 + rain933 = rain933 / 1000.
1.222 + rainmod933 = (rainmod933/ 1000.) * nsec_per_year
1.223 + nppmod933 = nppmod933 * nsec_per_year
1.224 +
1.225 + nppmod = nppmod * nsec_per_year
1.226 +
1.227 + npp81@units = "gC/m^2/yr"
1.228 + nppmod81@units = "gC/m^2/yr"
1.229 + npp933@units = "gC/m^2/yr"
1.230 + nppmod933@units = "gC/m^2/yr"
1.231 + nppmod@units = "gC/m^2/yr"
1.232 + nppglobe@units = "gC/m^2/yr"
1.233 + rain81@units = "m/yr"
1.234 + rainmod81@units = "m/yr"
1.235 + rain933@units = "m/yr"
1.236 + rainmod933@units = "m/yr"
1.237 +
1.238 + npp81@long_name = "NPP (gC/m2/year)"
1.239 + npp933@long_name = "NPP (gC/m2/year)"
1.240 + nppmod81@long_name = "NPP (gC/m2/year)"
1.241 + nppmod933@long_name = "NPP (gC/m2/year)"
1.242 + nppmod@long_name = "NPP (gC/m2/year)"
1.243 + nppglobe@long_name = "NPP (gC/m2/year)"
1.244 + rain81@long_name = "PREC (m/year)"
1.245 + rain933@long_name = "PREC (m/year)"
1.246 + rainmod81@long_name = "PREC (m/year)"
1.247 + rainmod933@long_name = "PREC (m/year)"
1.248 +
1.249 +;(A)-1 plot scatter ob 81
1.250 +
1.251 +;plot_name = "scatter_ob_81"
1.252 +;title = "Observed 81 sites"
1.253 +
1.254 +;wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.255 +;res = True ; plot mods desired
1.256 +;res@tiMainString = title ; add title
1.257 +;res@xyMarkLineModes = "Markers" ; choose which have markers
1.258 +;res@xyMarkers = 16 ; choose type of marker
1.259 +;res@xyMarkerColor = "red" ; Marker color
1.260 +;res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.261 +;res@tmLabelAutoStride = True ; nice tick mark labels
1.262 +
1.263 +;plot = gsn_csm_xy (wks,id81,npp81,res) ; create plot
1.264 +;frame(wks)
1.265 +
1.266 +;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.267 +;system("rm "+plot_name+"."+plot_type)
1.268 +;system("rm "+plot_name+"-1."+plot_type_new)
1.269 +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.270 +
1.271 +;clear (wks)
1.272 +;*******************************************************************
1.273 +;(A)-2 for table of site(1) (the first 41 sites)
1.274 +;*******************************************************************
1.275 +
1.276 + table_length = 0.95
1.277 +
1.278 +; table header name
1.279 + table_header_name = "Site ID"
1.280 +
1.281 +; column (not including header column)
1.282 + col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
1.283 + ncol = dimsizes(col_header)
1.284 +
1.285 +; row (not including header row)
1.286 + nrow = 41
1.287 + row_header = new ((/nrow/),string )
1.288 + row_header(0:nrow-1) = id81(0:nrow-1)
1.289 +
1.290 +; Table header
1.291 + ncr1 = (/1,1/) ; 1 row, 1 column
1.292 + x1 = (/0.005,0.15/) ; Start and end X
1.293 + y1 = (/0.900,0.95/) ; Start and end Y
1.294 + text1 = table_header_name
1.295 + res1 = True
1.296 + res1@txFontHeightF = 0.015
1.297 + res1@gsFillColor = "CornFlowerBlue"
1.298 +
1.299 +; Column header (equally space in x2)
1.300 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.301 + x2 = (/x1(1),0.995/) ; start from end of x1
1.302 + y2 = y1 ; same as y1
1.303 + text2 = col_header
1.304 + res2 = True
1.305 + res2@txFontHeightF = 0.015
1.306 + res2@gsFillColor = "Gray"
1.307 +
1.308 +; Row header (equally space in y2)
1.309 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.310 + x3 = x1 ; same as x1
1.311 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.312 + text3 = row_header
1.313 + res3 = True
1.314 + res3@txFontHeightF = 0.010
1.315 + res3@gsFillColor = "Gray"
1.316 +
1.317 +; Main table body
1.318 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.319 + x4 = x2 ; Start and end x
1.320 + y4 = y3 ; Start and end Y
1.321 + text4 = new((/nrow,ncol/),string)
1.322 +
1.323 + color_fill4 = new((/nrow,ncol/),string)
1.324 + color_fill4 = "white"
1.325 +; color_fill4(:,ncol-1) = "grey"
1.326 +; color_fill4(nrow-1,:) = "green"
1.327 +
1.328 + res4 = True ; Set up resource list
1.329 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.330 + res4@txFontHeightF = 0.015
1.331 + res4@gsFillColor = color_fill4
1.332 +
1.333 + delete (color_fill4)
1.334 +;-------------------------------------------------------------------
1.335 +; for table value
1.336 +
1.337 + do n = 0,nrow-1
1.338 + text4(n,0) = sprintf("%5.2f", y81(n))
1.339 + text4(n,1) = sprintf("%5.2f", x81(n))
1.340 + text4(n,2) = sprintf("%5.2f", npp81(n))
1.341 + text4(n,3) = sprintf("%5.2f", rain81(n))
1.342 + end do
1.343 +;---------------------------------------------------------------------------
1.344 +
1.345 + plot_name = "table_site81_ob1"
1.346 +
1.347 + wks = gsn_open_wks (plot_type,plot_name)
1.348 +;------------------------------------------
1.349 +; for table title
1.350 +
1.351 + gRes = True
1.352 + gRes@txFontHeightF = 0.02
1.353 +; gRes@txAngleF = 90
1.354 +
1.355 + title_text = "Observation at 81 Sites (1)"
1.356 +
1.357 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.358 +;------------------------------------------
1.359 +
1.360 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.361 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.362 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.363 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.364 +
1.365 + frame(wks)
1.366 + clear (wks)
1.367 +
1.368 + delete (row_header)
1.369 + delete (res3)
1.370 + delete (ncr3)
1.371 + delete (text3)
1.372 + delete (res4)
1.373 + delete (ncr4)
1.374 + delete (text4)
1.375 +
1.376 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.377 + system("rm "+plot_name+"."+plot_type)
1.378 + system("rm "+plot_name+"-1."+plot_type_new)
1.379 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.380 +
1.381 +;*******************************************************************
1.382 +;(A)-3 for table of site(2) (the second 40 sites)
1.383 +;*******************************************************************
1.384 +
1.385 + table_length = 0.95
1.386 +
1.387 +; table header name
1.388 + table_header_name = "Site ID"
1.389 +
1.390 +; column (not including header column)
1.391 + col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
1.392 + ncol = dimsizes(col_header)
1.393 +
1.394 +; row (not including header row)
1.395 + nrow = 40
1.396 + row_header = new ((/nrow/),string )
1.397 + row_header(0:nrow-1) = id81(41:41+nrow-1)
1.398 +
1.399 +; Table header
1.400 + ncr1 = (/1,1/) ; 1 row, 1 column
1.401 + x1 = (/0.005,0.15/) ; Start and end X
1.402 + y1 = (/0.900,0.95/) ; Start and end Y
1.403 + text1 = table_header_name
1.404 + res1 = True
1.405 + res1@txFontHeightF = 0.015
1.406 + res1@gsFillColor = "CornFlowerBlue"
1.407 +
1.408 +; Column header (equally space in x2)
1.409 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.410 + x2 = (/x1(1),0.995/) ; start from end of x1
1.411 + y2 = y1 ; same as y1
1.412 + text2 = col_header
1.413 + res2 = True
1.414 + res2@txFontHeightF = 0.015
1.415 + res2@gsFillColor = "Gray"
1.416 +
1.417 +; Row header (equally space in y2)
1.418 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.419 + x3 = x1 ; same as x1
1.420 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.421 + text3 = row_header
1.422 + res3 = True
1.423 + res3@txFontHeightF = 0.010
1.424 + res3@gsFillColor = "Gray"
1.425 +
1.426 +; Main table body
1.427 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.428 + x4 = x2 ; Start and end x
1.429 + y4 = y3 ; Start and end Y
1.430 + text4 = new((/nrow,ncol/),string)
1.431 +
1.432 + color_fill4 = new((/nrow,ncol/),string)
1.433 + color_fill4 = "white"
1.434 +; color_fill4(:,ncol-1) = "grey"
1.435 +; color_fill4(nrow-1,:) = "green"
1.436 +
1.437 + res4 = True ; Set up resource list
1.438 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.439 + res4@txFontHeightF = 0.015
1.440 + res4@gsFillColor = color_fill4
1.441 +
1.442 + delete (color_fill4)
1.443 +;-------------------------------------------------------------------
1.444 +; for table value
1.445 +
1.446 + do n = 0,nrow-1
1.447 + text4(n,0) = sprintf("%5.2f", y81(n+41))
1.448 + text4(n,1) = sprintf("%5.2f", x81(n+41))
1.449 + text4(n,2) = sprintf("%5.2f", npp81(n+41))
1.450 + text4(n,3) = sprintf("%5.2f", rain81(n+41))
1.451 + end do
1.452 +;---------------------------------------------------------------------------
1.453 +
1.454 + plot_name = "table_site81_ob2"
1.455 +
1.456 + wks = gsn_open_wks (plot_type,plot_name)
1.457 +;------------------------------------------
1.458 +; for table title
1.459 +
1.460 + gRes = True
1.461 + gRes@txFontHeightF = 0.02
1.462 +; gRes@txAngleF = 90
1.463 +
1.464 + title_text = "Observation at 81 Sites (2)"
1.465 +
1.466 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.467 +;------------------------------------------
1.468 +
1.469 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.470 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.471 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.472 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.473 +
1.474 + frame(wks)
1.475 + clear (wks)
1.476 +
1.477 + delete (row_header)
1.478 + delete (res3)
1.479 + delete (ncr3)
1.480 + delete (text3)
1.481 + delete (res4)
1.482 + delete (ncr4)
1.483 + delete (text4)
1.484 +
1.485 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.486 + system("rm "+plot_name+"."+plot_type)
1.487 + system("rm "+plot_name+"-1."+plot_type_new)
1.488 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.489 +;**************************************************************************
1.490 +;(A)-4 plot scatter ob 933
1.491 +
1.492 + plot_name = "scatter_ob_933"
1.493 + title = "Observed 933 sites"
1.494 +
1.495 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.496 + res = True ; plot mods desired
1.497 + res@tiMainString = title ; add title
1.498 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.499 + res@xyMarkers = 16 ; choose type of marker
1.500 + res@xyMarkerColor = "red" ; Marker color
1.501 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.502 + res@tmLabelAutoStride = True ; nice tick mark labels
1.503 +
1.504 + plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
1.505 + frame(wks)
1.506 +
1.507 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.508 + system("rm "+plot_name+"."+plot_type)
1.509 + system("rm "+plot_name+"-1."+plot_type_new)
1.510 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.511 +
1.512 + clear (wks)
1.513 +
1.514 +;(A)-5 plot scatter model 81
1.515 +
1.516 +;plot_name = "scatter_model_81"
1.517 +;title = "Model "+ model_name +" 81 sites"
1.518 +
1.519 +;wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.520 +;res = True ; plot mods desired
1.521 +;res@tiMainString = title ; add title
1.522 +;res@xyMarkLineModes = "Markers" ; choose which have markers
1.523 +;res@xyMarkers = 16 ; choose type of marker
1.524 +;res@xyMarkerColor = "red" ; Marker color
1.525 +;res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.526 +;res@tmLabelAutoStride = True ; nice tick mark labels
1.527 +
1.528 +;plot = gsn_csm_xy (wks,id81,nppmod81,res) ; create plot
1.529 +;frame(wks)
1.530 +
1.531 +;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.532 +;system("rm "+plot_name+"."+plot_type)
1.533 +;system("rm "+plot_name+"-1."+plot_type_new)
1.534 +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.535 +
1.536 +;clear (wks)
1.537 +
1.538 +;(A)-6 plot scatter model 933
1.539 +
1.540 + plot_name = "scatter_model_933"
1.541 + title = "Model "+ model_name +" 933 sites"
1.542 +
1.543 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.544 + res = True ; plot mods desired
1.545 + res@tiMainString = title ; add title
1.546 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.547 + res@xyMarkers = 16 ; choose type of marker
1.548 + res@xyMarkerColor = "red" ; Marker color
1.549 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.550 + res@tmLabelAutoStride = True ; nice tick mark labels
1.551 +
1.552 + plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
1.553 + frame(wks)
1.554 +
1.555 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.556 + system("rm "+plot_name+"."+plot_type)
1.557 + system("rm "+plot_name+"-1."+plot_type_new)
1.558 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.559 +
1.560 + clear (wks)
1.561 +
1.562 +;(A)-7 scatter model vs ob 81
1.563 +
1.564 + plot_name = "scatter_model_vs_ob_81"
1.565 + title = model_name +" vs ob 81 sites"
1.566 +
1.567 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.568 + res = True ; plot mods desired
1.569 + res@tiMainString = title ; add title
1.570 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.571 + res@xyMarkers = 16 ; choose type of marker
1.572 + res@xyMarkerColor = "red" ; Marker color
1.573 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.574 + res@tmLabelAutoStride = True ; nice tick mark labels
1.575 +
1.576 + res@gsnDraw = False
1.577 + res@gsnFrame = False ; don't advance frame yet
1.578 +;-------------------------------
1.579 +;compute correlation coef. and M
1.580 + ccr81 = esccr(nppmod81,npp81,0)
1.581 +;print (ccr81)
1.582 +
1.583 +;new eq
1.584 + bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
1.585 + M81 = (1. - (bias/dimsizes(y81)))*5.
1.586 + print (M81)
1.587 +
1.588 + tRes = True
1.589 + tRes@txFontHeightF = 0.025
1.590 +
1.591 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.592 +
1.593 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.594 +;--------------------------------
1.595 + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
1.596 +;-------------------------------
1.597 +; add polyline
1.598 +
1.599 + dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
1.600 +;-------------------------------
1.601 + draw (plot)
1.602 + frame(wks)
1.603 +
1.604 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.605 + system("rm "+plot_name+"."+plot_type)
1.606 +;system("rm "+plot_name+"-1."+plot_type_new)
1.607 +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.608 +
1.609 + clear (wks)
1.610 + delete (dum)
1.611 +
1.612 +;(A)-8 scatter model vs ob 933
1.613 +
1.614 + plot_name = "scatter_model_vs_ob_933"
1.615 + title = model_name +" vs ob 933 sites"
1.616 +
1.617 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.618 + res = True ; plot mods desired
1.619 + res@tiMainString = title ; add title
1.620 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.621 + res@xyMarkers = 16 ; choose type of marker
1.622 + res@xyMarkerColor = "red" ; Marker color
1.623 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.624 + res@tmLabelAutoStride = True ; nice tick mark labels
1.625 +
1.626 + res@gsnDraw = False
1.627 + res@gsnFrame = False ; don't advance frame yet
1.628 +;-------------------------------
1.629 +;compute correlation coef. and M
1.630 + ccr933 = esccr(nppmod933,npp933,0)
1.631 +;print (ccr933)
1.632 +
1.633 +;new eq
1.634 + bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
1.635 + M933 = (1. - (bias/dimsizes(y933)))*5.
1.636 + print (M933)
1.637 +
1.638 + tRes = True
1.639 + tRes@txFontHeightF = 0.025
1.640 +
1.641 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.642 +
1.643 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.644 +;--------------------------------
1.645 + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
1.646 +;-------------------------------
1.647 +; add polyline
1.648 +
1.649 + dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
1.650 +;-------------------------------
1.651 + draw (plot)
1.652 + frame(wks)
1.653 +
1.654 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.655 + system("rm "+plot_name+"."+plot_type)
1.656 +;system("rm "+plot_name+"-1."+plot_type_new)
1.657 +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.658 +
1.659 + clear (wks)
1.660 + delete (dum)
1.661 +
1.662 +;*******************************************************************
1.663 +;(A)-9 for table of site(1) (the first 41 sites), model vs ob
1.664 +;*******************************************************************
1.665 +
1.666 + table_length = 0.95
1.667 +
1.668 +; table header name
1.669 + table_header_name = "Site ID"
1.670 +
1.671 +; row (not including header row)
1.672 + nrow = 41
1.673 + row_header = new ((/nrow/),string )
1.674 + row_header(0:nrow-1) = id81(0:nrow-1)
1.675 +
1.676 +; Table header
1.677 + ncr1 = (/1,1/) ; 1 row, 1 column
1.678 + x1 = (/0.005,0.15/) ; Start and end X
1.679 + y1 = (/0.85 ,0.95/) ; Start and end Y
1.680 + text1 = table_header_name
1.681 + res1 = True
1.682 + res1@txFontHeightF = 0.015
1.683 + res1@gsFillColor = "CornFlowerBlue"
1.684 +
1.685 +; Column header1 (equally space in x2)
1.686 +
1.687 + delete (col_header)
1.688 + delete (ncr2)
1.689 + delete (text2)
1.690 +
1.691 + col_header = (/"Latitude","Longitude"/)
1.692 + ncol = dimsizes(col_header)
1.693 +
1.694 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.695 + x2 = (/x1(1),0.35/) ; start from end of x1
1.696 + y2 = y1 ; same as y1
1.697 + text2 = col_header
1.698 + res2 = True
1.699 + res2@txFontHeightF = 0.015
1.700 + res2@gsFillColor = "Gray"
1.701 +
1.702 +; Column header2 (equally space in x2)
1.703 +
1.704 + col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
1.705 + col_header2 = (/"ob",model_name \
1.706 + ,"ob",model_name \
1.707 + /)
1.708 +
1.709 + ncol1 = dimsizes(col_header1)
1.710 + ncol2 = dimsizes(col_header2)
1.711 +
1.712 + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
1.713 + x21 = (/x2(1),0.995/) ; start from end of x2
1.714 + yhalf = (y1(0)+y1(1))*0.5
1.715 + y21 = (/yhalf,y1(1)/) ; top half of y1
1.716 + text21 = col_header1
1.717 + res21 = True
1.718 + res21@txFontHeightF = 0.015
1.719 + res21@gsFillColor = "Gray"
1.720 +
1.721 + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
1.722 + x22 = x21 ; start from end of x1
1.723 + y22 = (/y1(0),yhalf/) ; bottom half of y1
1.724 + text22 = col_header2
1.725 + res22 = True
1.726 + res22@txFontHeightF = 0.015
1.727 + res22@gsFillColor = "Gray"
1.728 +
1.729 +; Row header (equally space in y2)
1.730 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.731 + x3 = x1 ; same as x1
1.732 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.733 + text3 = row_header
1.734 + res3 = True
1.735 + res3@txFontHeightF = 0.010
1.736 + res3@gsFillColor = "Gray"
1.737 +
1.738 +; Main table body-1
1.739 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.740 + x4 = x2 ; Start and end x
1.741 + y4 = y3 ; Start and end Y
1.742 + text4 = new((/nrow,ncol/),string)
1.743 +
1.744 + color_fill4 = new((/nrow,ncol/),string)
1.745 + color_fill4 = "white"
1.746 +; color_fill4(:,ncol-1) = "grey"
1.747 +; color_fill4(nrow-1,:) = "green"
1.748 +
1.749 + res4 = True ; Set up resource list
1.750 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.751 + res4@txFontHeightF = 0.015
1.752 + res4@gsFillColor = color_fill4
1.753 +
1.754 + delete (color_fill4)
1.755 +
1.756 +; Main table body-2
1.757 + ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
1.758 + x42 = x21 ; Start and end x
1.759 + y42 = y3 ; Start and end Y
1.760 + text42 = new((/nrow,ncol2/),string)
1.761 +
1.762 + color_fill42 = new((/nrow,ncol2/),string)
1.763 + color_fill42 = "white"
1.764 +; color_fill42(:,ncol-1) = "grey"
1.765 +; color_fill42(nrow-1,:) = "green"
1.766 +
1.767 + res42 = True ; Set up resource list
1.768 +; res42@gsnDebug = True ; Useful to print NDC row,col values used.
1.769 + res42@txFontHeightF = 0.015
1.770 + res42@gsFillColor = color_fill42
1.771 +
1.772 + delete (color_fill42)
1.773 +;-------------------------------------------------------------------
1.774 +; for table value
1.775 +
1.776 + do n = 0,nrow-1
1.777 + text4(n,0) = sprintf("%5.2f", y81(n))
1.778 + text4(n,1) = sprintf("%5.2f", x81(n))
1.779 + end do
1.780 +
1.781 + do n = 0,nrow-1
1.782 + text42(n,0) = sprintf("%5.2f", npp81(n))
1.783 + text42(n,1) = sprintf("%5.2f", nppmod81(n))
1.784 + text42(n,2) = sprintf("%5.2f", rain81(n))
1.785 + text42(n,3) = sprintf("%5.2f", rainmod81(n))
1.786 + end do
1.787 +;---------------------------------------------------------------------------
1.788 +
1.789 + plot_name = "table_site81_model_vs_ob1"
1.790 +
1.791 + wks = gsn_open_wks (plot_type,plot_name)
1.792 +;------------------------------------------
1.793 +; for table title
1.794 +
1.795 + gRes = True
1.796 + gRes@txFontHeightF = 0.02
1.797 +; gRes@txAngleF = 90
1.798 +
1.799 + title_text = "Model vs Observation at 81 Sites (1)"
1.800 +
1.801 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.802 +;------------------------------------------
1.803 +
1.804 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.805 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.806 + gsn_table(wks,ncr21,x21,y21,text21,res21)
1.807 + gsn_table(wks,ncr22,x22,y22,text22,res22)
1.808 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.809 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.810 + gsn_table(wks,ncr42,x42,y42,text42,res42)
1.811 +
1.812 + frame(wks)
1.813 + clear (wks)
1.814 +
1.815 + delete (col_header)
1.816 + delete (row_header)
1.817 + delete (res1)
1.818 + delete (ncr1)
1.819 + delete (text1)
1.820 + delete (res2)
1.821 + delete (ncr2)
1.822 + delete (text2)
1.823 + delete (res21)
1.824 + delete (ncr21)
1.825 + delete (text21)
1.826 + delete (res22)
1.827 + delete (ncr22)
1.828 + delete (text22)
1.829 + delete (res3)
1.830 + delete (ncr3)
1.831 + delete (text3)
1.832 + delete (res4)
1.833 + delete (ncr4)
1.834 + delete (text4)
1.835 + delete (res42)
1.836 + delete (ncr42)
1.837 + delete (text42)
1.838 +
1.839 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.840 + system("rm "+plot_name+"."+plot_type)
1.841 + system("rm "+plot_name+"-1."+plot_type_new)
1.842 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.843 +
1.844 +;*******************************************************************
1.845 +;(A)-10 for table of site(2) (the last 40 sites), model vs ob
1.846 +;*******************************************************************
1.847 +
1.848 + table_length = 0.95
1.849 +
1.850 +; table header name
1.851 + table_header_name = "Site ID"
1.852 +
1.853 +; row (not including header row)
1.854 + nrow = 40
1.855 + row_header = new ((/nrow/),string )
1.856 + row_header(0:nrow-1) = id81(41:41+nrow-1)
1.857 +
1.858 +; Table header
1.859 + ncr1 = (/1,1/) ; 1 row, 1 column
1.860 + x1 = (/0.005,0.15/) ; Start and end X
1.861 + y1 = (/0.85 ,0.95/) ; Start and end Y
1.862 + text1 = table_header_name
1.863 + res1 = True
1.864 + res1@txFontHeightF = 0.015
1.865 + res1@gsFillColor = "CornFlowerBlue"
1.866 +
1.867 +; Column header1 (equally space in x2)
1.868 +
1.869 + col_header = (/"Latitude","Longitude"/)
1.870 + ncol = dimsizes(col_header)
1.871 +
1.872 + ncr2 = (/1,ncol/) ; 1 rows, ncol columns
1.873 + x2 = (/x1(1),0.35/) ; start from end of x1
1.874 + y2 = y1 ; same as y1
1.875 + text2 = col_header
1.876 + res2 = True
1.877 + res2@txFontHeightF = 0.015
1.878 + res2@gsFillColor = "Gray"
1.879 +
1.880 +; Column header2 (equally space in x2)
1.881 +
1.882 +; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
1.883 + col_header2 = (/"ob",model_name \
1.884 + ,"ob",model_name \
1.885 + /)
1.886 +
1.887 + ncol1 = dimsizes(col_header1)
1.888 + ncol2 = dimsizes(col_header2)
1.889 +
1.890 + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
1.891 + x21 = (/x2(1),0.995/) ; start from end of x2
1.892 + yhalf = (y1(0)+y1(1))*0.5
1.893 + y21 = (/yhalf,y1(1)/) ; top half of y1
1.894 + text21 = col_header1
1.895 + res21 = True
1.896 + res21@txFontHeightF = 0.015
1.897 + res21@gsFillColor = "Gray"
1.898 +
1.899 + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
1.900 + x22 = x21 ; start from end of x1
1.901 + y22 = (/y1(0),yhalf/) ; bottom half of y1
1.902 + text22 = col_header2
1.903 + res22 = True
1.904 + res22@txFontHeightF = 0.015
1.905 + res22@gsFillColor = "Gray"
1.906 +
1.907 +; Row header (equally space in y2)
1.908 + ncr3 = (/nrow,1/) ; nrow rows, 1 columns
1.909 + x3 = x1 ; same as x1
1.910 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
1.911 + text3 = row_header
1.912 + res3 = True
1.913 + res3@txFontHeightF = 0.010
1.914 + res3@gsFillColor = "Gray"
1.915 +
1.916 +; Main table body-1
1.917 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
1.918 + x4 = x2 ; Start and end x
1.919 + y4 = y3 ; Start and end Y
1.920 + text4 = new((/nrow,ncol/),string)
1.921 +
1.922 + color_fill4 = new((/nrow,ncol/),string)
1.923 + color_fill4 = "white"
1.924 +; color_fill4(:,ncol-1) = "grey"
1.925 +; color_fill4(nrow-1,:) = "green"
1.926 +
1.927 + res4 = True ; Set up resource list
1.928 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.929 + res4@txFontHeightF = 0.015
1.930 + res4@gsFillColor = color_fill4
1.931 +
1.932 + delete (color_fill4)
1.933 +
1.934 +; Main table body-2
1.935 + ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
1.936 + x42 = x21 ; Start and end x
1.937 + y42 = y3 ; Start and end Y
1.938 + text42 = new((/nrow,ncol2/),string)
1.939 +
1.940 + color_fill42 = new((/nrow,ncol2/),string)
1.941 + color_fill42 = "white"
1.942 +; color_fill42(:,ncol-1) = "grey"
1.943 +; color_fill42(nrow-1,:) = "green"
1.944 +
1.945 + res42 = True ; Set up resource list
1.946 +; res42@gsnDebug = True ; Useful to print NDC row,col values used.
1.947 + res42@txFontHeightF = 0.015
1.948 + res42@gsFillColor = color_fill42
1.949 +
1.950 + delete (color_fill42)
1.951 +;-------------------------------------------------------------------
1.952 +; for table value
1.953 +
1.954 + do n = 0,nrow-1
1.955 + text4(n,0) = sprintf("%5.2f", y81(n+41))
1.956 + text4(n,1) = sprintf("%5.2f", x81(n+41))
1.957 + end do
1.958 +
1.959 + do n = 0,nrow-1
1.960 + text42(n,0) = sprintf("%5.2f", npp81(n+41))
1.961 + text42(n,1) = sprintf("%5.2f", nppmod81(n+41))
1.962 + text42(n,2) = sprintf("%5.2f", rain81(n+41))
1.963 + text42(n,3) = sprintf("%5.2f", rainmod81(n+41))
1.964 + end do
1.965 +;---------------------------------------------------------------------------
1.966 +
1.967 + plot_name = "table_site81_model_vs_ob2"
1.968 +
1.969 + wks = gsn_open_wks (plot_type,plot_name)
1.970 +;------------------------------------------
1.971 +; for table title
1.972 +
1.973 + gRes = True
1.974 + gRes@txFontHeightF = 0.02
1.975 +; gRes@txAngleF = 90
1.976 +
1.977 + title_text = "Model vs Observation at 81 Sites (2)"
1.978 +
1.979 + gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
1.980 +;------------------------------------------
1.981 +
1.982 + gsn_table(wks,ncr1,x1,y1,text1,res1)
1.983 + gsn_table(wks,ncr2,x2,y2,text2,res2)
1.984 + gsn_table(wks,ncr21,x21,y21,text21,res21)
1.985 + gsn_table(wks,ncr22,x22,y22,text22,res22)
1.986 + gsn_table(wks,ncr3,x3,y3,text3,res3)
1.987 + gsn_table(wks,ncr4,x4,y4,text4,res4)
1.988 + gsn_table(wks,ncr42,x42,y42,text42,res42)
1.989 +
1.990 + frame(wks)
1.991 + clear (wks)
1.992 +
1.993 + delete (col_header)
1.994 + delete (row_header)
1.995 + delete (res1)
1.996 + delete (ncr1)
1.997 + delete (text1)
1.998 + delete (res2)
1.999 + delete (ncr2)
1.1000 + delete (text2)
1.1001 + delete (res21)
1.1002 + delete (ncr21)
1.1003 + delete (text21)
1.1004 + delete (res22)
1.1005 + delete (ncr22)
1.1006 + delete (text22)
1.1007 + delete (res3)
1.1008 + delete (ncr3)
1.1009 + delete (text3)
1.1010 + delete (res4)
1.1011 + delete (ncr4)
1.1012 + delete (text4)
1.1013 + delete (res42)
1.1014 + delete (ncr42)
1.1015 + delete (text42)
1.1016 +
1.1017 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1018 + system("rm "+plot_name+"."+plot_type)
1.1019 + system("rm "+plot_name+"-1."+plot_type_new)
1.1020 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1021 +
1.1022 +;**************************************************************************
1.1023 +;(B) histogram 81
1.1024 +;***********************************************************************-
1.1025 +;get data
1.1026 +
1.1027 + RAIN1_1D = ndtooned(rain81)
1.1028 + RAIN2_1D = ndtooned(rainmod81)
1.1029 + NPP1_1D = ndtooned(npp81)
1.1030 + NPP2_1D = ndtooned(nppmod81)
1.1031 +
1.1032 + nx = 8
1.1033 +
1.1034 + xvalues = new((/2,nx/),float)
1.1035 + yvalues = new((/2,nx/),float)
1.1036 + mn_yvalues = new((/2,nx/),float)
1.1037 + mx_yvalues = new((/2,nx/),float)
1.1038 + dx4 = new((/1/),float)
1.1039 +
1.1040 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.1041 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.1042 +
1.1043 +;----------------------------------------
1.1044 +;compute correlation coeff and M score
1.1045 +
1.1046 + u = yvalues(0,:)
1.1047 + v = yvalues(1,:)
1.1048 +
1.1049 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.1050 + uu = u(good)
1.1051 + vv = v(good)
1.1052 +
1.1053 + ccr81h = esccr(uu,vv,0)
1.1054 +;print (ccr81h)
1.1055 +
1.1056 +;new eq
1.1057 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.1058 + M81h = (1.- (bias/dimsizes(uu)))*5.
1.1059 + print (M81h)
1.1060 +
1.1061 + delete (u)
1.1062 + delete (v)
1.1063 + delete (uu)
1.1064 + delete (vv)
1.1065 +;----------------------------------------------------------------------
1.1066 +; histogram res
1.1067 +
1.1068 + resh = True
1.1069 + resh@gsnMaximize = True
1.1070 + resh@gsnDraw = False
1.1071 + resh@gsnFrame = False
1.1072 + resh@xyMarkLineMode = "Markers"
1.1073 + resh@xyMarkerSizeF = 0.014
1.1074 + resh@xyMarker = 16
1.1075 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.1076 + resh@trYMinF = min(mn_yvalues) - 10.
1.1077 + resh@trYMaxF = max(mx_yvalues) + 10.
1.1078 +
1.1079 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.1080 + resh@tiXAxisString = "Precipitation (m/year)"
1.1081 +
1.1082 + max_bar = new((/2,nx/),graphic)
1.1083 + min_bar = new((/2,nx/),graphic)
1.1084 + max_cap = new((/2,nx/),graphic)
1.1085 + min_cap = new((/2,nx/),graphic)
1.1086 +
1.1087 + lnres = True
1.1088 + line_colors = (/"brown","blue"/)
1.1089 +;=================================================================
1.1090 +;(B)-1 histogram ob 81 site only
1.1091 +;
1.1092 + plot_name = "histogram_ob_81"
1.1093 + title = "Observed 81 site"
1.1094 + resh@tiMainString = title
1.1095 +
1.1096 + wks = gsn_open_wks (plot_type,plot_name)
1.1097 +
1.1098 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.1099 +
1.1100 +;-------------------------------
1.1101 +;Attach the vertical bar and the horizontal cap line
1.1102 +
1.1103 + delete (x1)
1.1104 + delete (x2)
1.1105 + delete (y1)
1.1106 + delete (y2)
1.1107 +
1.1108 + do nd=0,0
1.1109 + lnres@gsLineColor = line_colors(nd)
1.1110 + do i=0,nx-1
1.1111 +
1.1112 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1113 + .not.ismissing(mx_yvalues(nd,i))) then
1.1114 +;
1.1115 +; Attach the vertical bar, both above and below the marker.
1.1116 +;
1.1117 + x1 = xvalues(nd,i)
1.1118 + y1 = yvalues(nd,i)
1.1119 + y2 = mn_yvalues(nd,i)
1.1120 + plx = (/x1,x1/)
1.1121 + ply = (/y1,y2/)
1.1122 + min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1123 +
1.1124 + y2 = mx_yvalues(nd,i)
1.1125 + plx = (/x1,x1/)
1.1126 + ply = (/y1,y2/)
1.1127 + max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1128 +;
1.1129 +; Attach the horizontal cap line, both above and below the marker.
1.1130 +;
1.1131 + x1 = xvalues(nd,i) - dx4
1.1132 + x2 = xvalues(nd,i) + dx4
1.1133 + y1 = mn_yvalues(nd,i)
1.1134 + plx = (/x1,x2/)
1.1135 + ply = (/y1,y1/)
1.1136 + min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1137 +
1.1138 + y1 = mx_yvalues(nd,i)
1.1139 + plx = (/x1,x2/)
1.1140 + ply = (/y1,y1/)
1.1141 + max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.1142 + end if
1.1143 + end do
1.1144 + end do
1.1145 +
1.1146 + draw(xy)
1.1147 + frame(wks)
1.1148 +
1.1149 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1150 + system("rm "+plot_name+"."+plot_type)
1.1151 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1152 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1153 +
1.1154 + clear (wks)
1.1155 +;===========================================================================
1.1156 +;(B)-2 histogram model vs ob 81 site
1.1157 +
1.1158 + plot_name = "histogram_mod-ob_81"
1.1159 + title = model_name+ " vs Observed 81 site"
1.1160 + resh@tiMainString = title
1.1161 +
1.1162 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1163 +
1.1164 +;-----------------------------
1.1165 +; Add a boxed legend using the more simple method, which won't have
1.1166 +; vertical lines going through the markers.
1.1167 +
1.1168 + resh@pmLegendDisplayMode = "Always"
1.1169 +; resh@pmLegendWidthF = 0.1
1.1170 + resh@pmLegendWidthF = 0.08
1.1171 + resh@pmLegendHeightF = 0.05
1.1172 + resh@pmLegendOrthogonalPosF = -1.17
1.1173 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1174 +; resh@pmLegendParallelPosF = 0.18
1.1175 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.1176 +
1.1177 +; resh@lgPerimOn = False
1.1178 + resh@lgLabelFontHeightF = 0.015
1.1179 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.1180 +;-----------------------------
1.1181 + tRes = True
1.1182 + tRes@txFontHeightF = 0.025
1.1183 +
1.1184 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1.1185 +
1.1186 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1187 +
1.1188 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.1189 +;-------------------------------
1.1190 +;Attach the vertical bar and the horizontal cap line
1.1191 +
1.1192 + do nd=0,1
1.1193 + lnres@gsLineColor = line_colors(nd)
1.1194 + do i=0,nx-1
1.1195 +
1.1196 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1197 + .not.ismissing(mx_yvalues(nd,i))) then
1.1198 +;
1.1199 +; Attach the vertical bar, both above and below the marker.
1.1200 +;
1.1201 + x1 = xvalues(nd,i)
1.1202 + y1 = yvalues(nd,i)
1.1203 + y2 = mn_yvalues(nd,i)
1.1204 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1205 +
1.1206 + y2 = mx_yvalues(nd,i)
1.1207 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1208 +;
1.1209 +; Attach the horizontal cap line, both above and below the marker.
1.1210 +;
1.1211 + x1 = xvalues(nd,i) - dx4
1.1212 + x2 = xvalues(nd,i) + dx4
1.1213 + y1 = mn_yvalues(nd,i)
1.1214 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1215 +
1.1216 + y1 = mx_yvalues(nd,i)
1.1217 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1218 + end if
1.1219 + end do
1.1220 + end do
1.1221 + draw(xy)
1.1222 + frame(wks)
1.1223 +
1.1224 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1225 + system("rm "+plot_name+"."+plot_type)
1.1226 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1227 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1228 +
1.1229 + clear (wks)
1.1230 +
1.1231 + delete (RAIN1_1D)
1.1232 + delete (RAIN2_1D)
1.1233 + delete (NPP1_1D)
1.1234 + delete (NPP2_1D)
1.1235 +;delete (range)
1.1236 + delete (xvalues)
1.1237 + delete (yvalues)
1.1238 + delete (mn_yvalues)
1.1239 + delete (mx_yvalues)
1.1240 + delete (good)
1.1241 + delete (max_bar)
1.1242 + delete (min_bar)
1.1243 + delete (max_cap)
1.1244 + delete (min_cap)
1.1245 +;**************************************************************************
1.1246 +;(B) histogram 933
1.1247 +
1.1248 +;--------------------------------------------------------------------
1.1249 +;get data
1.1250 +
1.1251 + RAIN1_1D = ndtooned(rain933)
1.1252 + RAIN2_1D = ndtooned(rainmod933)
1.1253 + NPP1_1D = ndtooned(npp933)
1.1254 + NPP2_1D = ndtooned(nppmod933)
1.1255 +
1.1256 + nx = 10
1.1257 + xvalues = new((/2,nx/),float)
1.1258 + yvalues = new((/2,nx/),float)
1.1259 + mn_yvalues = new((/2,nx/),float)
1.1260 + mx_yvalues = new((/2,nx/),float)
1.1261 + dx4 = new((/1/),float)
1.1262 +
1.1263 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.1264 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.1265 +
1.1266 +;----------------------------------------
1.1267 +;compute correlation coeff and M score
1.1268 +
1.1269 + u = yvalues(0,:)
1.1270 + v = yvalues(1,:)
1.1271 +
1.1272 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.1273 + uu = u(good)
1.1274 + vv = v(good)
1.1275 +
1.1276 + ccr933h = esccr(uu,vv,0)
1.1277 +;print (ccr933h)
1.1278 +
1.1279 +;new eq
1.1280 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.1281 + M933h = (1.- (bias/dimsizes(uu)))*5.
1.1282 + print (M933h)
1.1283 +
1.1284 + delete (u)
1.1285 + delete (v)
1.1286 + delete (uu)
1.1287 + delete (vv)
1.1288 +;----------------------------------------------------------------------
1.1289 +; histogram res
1.1290 +
1.1291 + delete (resh)
1.1292 + resh = True
1.1293 + resh@gsnMaximize = True
1.1294 + resh@gsnDraw = False
1.1295 + resh@gsnFrame = False
1.1296 + resh@xyMarkLineMode = "Markers"
1.1297 + resh@xyMarkerSizeF = 0.014
1.1298 + resh@xyMarker = 16
1.1299 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.1300 + resh@trYMinF = min(mn_yvalues) - 10.
1.1301 + resh@trYMaxF = max(mx_yvalues) + 10.
1.1302 +
1.1303 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.1304 + resh@tiXAxisString = "Precipitation (m/year)"
1.1305 +
1.1306 + max_bar = new((/2,nx/),graphic)
1.1307 + min_bar = new((/2,nx/),graphic)
1.1308 + max_cap = new((/2,nx/),graphic)
1.1309 + min_cap = new((/2,nx/),graphic)
1.1310 +
1.1311 + lnres = True
1.1312 + line_colors = (/"brown","blue"/)
1.1313 +;=================================================================
1.1314 +;(B)-3 histogram ob 933 site only
1.1315 +
1.1316 + plot_name = "histogram_ob_933"
1.1317 + title = "Observed 933 site"
1.1318 + resh@tiMainString = title
1.1319 +
1.1320 + wks = gsn_open_wks (plot_type,plot_name)
1.1321 +
1.1322 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.1323 +
1.1324 +;-------------------------------
1.1325 +;Attach the vertical bar and the horizontal cap line
1.1326 +
1.1327 + do nd=0,0
1.1328 + lnres@gsLineColor = line_colors(nd)
1.1329 + do i=0,nx-1
1.1330 +
1.1331 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1332 + .not.ismissing(mx_yvalues(nd,i))) then
1.1333 +
1.1334 +; Attach the vertical bar, both above and below the marker.
1.1335 +
1.1336 + x1 = xvalues(nd,i)
1.1337 + y1 = yvalues(nd,i)
1.1338 + y2 = mn_yvalues(nd,i)
1.1339 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1340 +
1.1341 + y2 = mx_yvalues(nd,i)
1.1342 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1343 +
1.1344 +; Attach the horizontal cap line, both above and below the marker.
1.1345 +
1.1346 + x1 = xvalues(nd,i) - dx4
1.1347 + x2 = xvalues(nd,i) + dx4
1.1348 + y1 = mn_yvalues(nd,i)
1.1349 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1350 +
1.1351 + y1 = mx_yvalues(nd,i)
1.1352 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1353 + end if
1.1354 + end do
1.1355 + end do
1.1356 +
1.1357 + draw(xy)
1.1358 + frame(wks)
1.1359 +
1.1360 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1361 + system("rm "+plot_name+"."+plot_type)
1.1362 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1363 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1364 +
1.1365 + delete (xy)
1.1366 + clear (wks)
1.1367 +
1.1368 +;===========================================================================
1.1369 +;(B)-4 histogram model vs ob 933 site
1.1370 +
1.1371 + plot_name = "histogram_mod-ob_933"
1.1372 + title = model_name+ " vs Observed 933 site"
1.1373 + resh@tiMainString = title
1.1374 +
1.1375 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1376 +
1.1377 +;-----------------------------
1.1378 +; Add a boxed legend using the more simple method, which won't have
1.1379 +; vertical lines going through the markers.
1.1380 +
1.1381 + resh@pmLegendDisplayMode = "Always"
1.1382 +; resh@pmLegendWidthF = 0.1
1.1383 + resh@pmLegendWidthF = 0.08
1.1384 + resh@pmLegendHeightF = 0.05
1.1385 + resh@pmLegendOrthogonalPosF = -1.17
1.1386 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1387 +; resh@pmLegendParallelPosF = 0.18
1.1388 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.1389 +
1.1390 +; resh@lgPerimOn = False
1.1391 + resh@lgLabelFontHeightF = 0.015
1.1392 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.1393 +;-----------------------------
1.1394 + tRes = True
1.1395 + tRes@txFontHeightF = 0.025
1.1396 +
1.1397 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1.1398 +
1.1399 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1400 +
1.1401 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.1402 +;-------------------------------
1.1403 +;Attach the vertical bar and the horizontal cap line
1.1404 +
1.1405 + do nd=0,1
1.1406 + lnres@gsLineColor = line_colors(nd)
1.1407 + do i=0,nx-1
1.1408 +
1.1409 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1410 + .not.ismissing(mx_yvalues(nd,i))) then
1.1411 +;
1.1412 +; Attach the vertical bar, both above and below the marker.
1.1413 +;
1.1414 + x1 = xvalues(nd,i)
1.1415 + y1 = yvalues(nd,i)
1.1416 + y2 = mn_yvalues(nd,i)
1.1417 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1418 +
1.1419 + y2 = mx_yvalues(nd,i)
1.1420 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1421 +;
1.1422 +; Attach the horizontal cap line, both above and below the marker.
1.1423 +;
1.1424 + x1 = xvalues(nd,i) - dx4
1.1425 + x2 = xvalues(nd,i) + dx4
1.1426 + y1 = mn_yvalues(nd,i)
1.1427 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1428 +
1.1429 + y1 = mx_yvalues(nd,i)
1.1430 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1431 + end if
1.1432 + end do
1.1433 + end do
1.1434 + draw(xy)
1.1435 + frame(wks)
1.1436 +
1.1437 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1438 + system("rm "+plot_name+"."+plot_type)
1.1439 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1440 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1441 +
1.1442 + delete (xy)
1.1443 + clear (wks)
1.1444 +;***********************************************************************-
1.1445 +; global contour
1.1446 +
1.1447 +;res
1.1448 + resg = True ; Use plot options
1.1449 + resg@cnFillOn = True ; Turn on color fill
1.1450 + resg@gsnSpreadColors = True ; use full colormap
1.1451 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1452 +; resg@lbLabelAutoStride = True
1.1453 + resg@cnLinesOn = False ; Turn off contourn lines
1.1454 + resg@mpFillOn = False ; Turn off map fill
1.1455 +
1.1456 + resg@gsnSpreadColors = True ; use full colormap
1.1457 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1458 + resg@cnMinLevelValF = 0. ; Min level
1.1459 + resg@cnMaxLevelValF = 2200. ; Max level
1.1460 + resg@cnLevelSpacingF = 200. ; interval
1.1461 +;------------------------------------------------------------------------
1.1462 +;(C)-1 global contour ob
1.1463 +
1.1464 + delta = 0.00000000001
1.1465 + nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1.1466 +
1.1467 + plot_name = "global_ob"
1.1468 + title = "Observed MODIS MOD 17"
1.1469 + resg@tiMainString = title
1.1470 +
1.1471 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1472 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1473 +
1.1474 + plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1475 + frame(wks)
1.1476 +
1.1477 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1478 + system("rm "+plot_name+"."+plot_type)
1.1479 + system("rm "+plot_name+"-1."+plot_type_new)
1.1480 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1481 +
1.1482 + clear (wks)
1.1483 +;------------------------------------------------------------------------
1.1484 +;(C)-2 global contour model
1.1485 +
1.1486 + plot_name = "global_model"
1.1487 + title = "Model "+ model_name
1.1488 + resg@tiMainString = title
1.1489 +
1.1490 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1491 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1492 +
1.1493 + plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1494 + frame(wks)
1.1495 +
1.1496 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1497 + system("rm "+plot_name+"."+plot_type)
1.1498 + system("rm "+plot_name+"-1."+plot_type_new)
1.1499 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1500 +
1.1501 + clear (wks)
1.1502 +;------------------------------------------------------------------------
1.1503 +;(C)-3 global contour model vs ob
1.1504 +
1.1505 + plot_name = "global_model_vs_ob"
1.1506 +
1.1507 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1508 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1509 +
1.1510 + delete (plot)
1.1511 + plot=new(3,graphic) ; create graphic array
1.1512 +
1.1513 + resg@gsnFrame = False ; Do not draw plot
1.1514 + resg@gsnDraw = False ; Do not advance frame
1.1515 +
1.1516 +;(d) compute correlation coef and M score
1.1517 +
1.1518 + uu1 = ndtooned(nppmod)
1.1519 + vv1 = ndtooned(nppglobe)
1.1520 +
1.1521 + delete (good)
1.1522 + good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1.1523 +
1.1524 + ug = uu1(good)
1.1525 + vg = vv1(good)
1.1526 +
1.1527 + ccrG = esccr(ug,vg,0)
1.1528 +; print (ccrG)
1.1529 +
1.1530 + MG = (ccrG*ccrG)* 5.0
1.1531 + print (MG)
1.1532 +
1.1533 +; plot correlation coef
1.1534 +
1.1535 + gRes = True
1.1536 + gRes@txFontHeightF = 0.02
1.1537 + gRes@txAngleF = 90
1.1538 +
1.1539 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.1540 +
1.1541 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1542 +;--------------------------------------------------------------------
1.1543 +
1.1544 +;(a) ob
1.1545 +
1.1546 + title = "Observed MODIS MOD 17"
1.1547 + resg@tiMainString = title
1.1548 +
1.1549 + plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1550 +
1.1551 +;(b) model
1.1552 +
1.1553 + title = "Model "+ model_name
1.1554 + resg@tiMainString = title
1.1555 +
1.1556 + plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1557 +
1.1558 +;(c) model-ob
1.1559 +
1.1560 + zz = nppmod
1.1561 + zz = nppmod - nppglobe
1.1562 + title = "Model_"+model_name+" - Observed"
1.1563 +
1.1564 + resg@cnMinLevelValF = -500 ; Min level
1.1565 + resg@cnMaxLevelValF = 500. ; Max level
1.1566 + resg@cnLevelSpacingF = 50. ; interval
1.1567 + resg@tiMainString = title
1.1568 +
1.1569 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1570 +
1.1571 + pres = True ; panel plot mods desired
1.1572 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1573 + ; indiv. plots in panel
1.1574 + pres@gsnMaximize = True ; fill the page
1.1575 +
1.1576 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1577 +
1.1578 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1579 + system("rm "+plot_name+"."+plot_type)
1.1580 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1581 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1582 +
1.1583 + frame (wks)
1.1584 + clear (wks)
1.1585 +
1.1586 + delete (plot)
1.1587 +;**********************************************************************
1.1588 +;(D)-1 zonal line plot ob
1.1589 +
1.1590 + resz = True
1.1591 +
1.1592 + vv = zonalAve(nppglobe)
1.1593 + vv@long_name = nppglobe@long_name
1.1594 +
1.1595 + plot_name = "zonal_ob"
1.1596 + title = "Observed MODIS MOD 17"
1.1597 + resz@tiMainString = title
1.1598 +
1.1599 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1600 +
1.1601 + plot = gsn_csm_xy (wks,ym,vv,resz)
1.1602 + frame(wks)
1.1603 +
1.1604 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1605 + system("rm "+plot_name+"."+plot_type)
1.1606 + system("rm "+plot_name+"-1."+plot_type_new)
1.1607 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1608 +
1.1609 + clear (wks)
1.1610 +;---------------------------------------------------------------------
1.1611 +;(D)-2 zonal line plot model vs ob
1.1612 +
1.1613 + s = new ((/2,dimsizes(ym)/), float)
1.1614 +
1.1615 + s(0,:) = zonalAve(nppglobe)
1.1616 + s(1,:) = zonalAve(nppmod)
1.1617 +
1.1618 + s@long_name = nppglobe@long_name
1.1619 +;-------------------------------------------
1.1620 +;(d) compute correlation coef and M score
1.1621 +
1.1622 + ccrZ = esccr(s(1,:), s(0,:),0)
1.1623 +; print (ccrZ)
1.1624 +
1.1625 + MZ = (ccrZ*ccrZ)* 5.0
1.1626 + print (MZ)
1.1627 +;-------------------------------------------
1.1628 + plot_name = "zonal_model_vs_ob"
1.1629 + title = "Zonal Average"
1.1630 + resz@tiMainString = title
1.1631 +
1.1632 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1633 +
1.1634 +; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1.1635 +; resz@vpWidthF = 0.7
1.1636 +
1.1637 + resz@xyMonoLineColor = "False" ; want colored lines
1.1638 + resz@xyLineColors = (/"black","red"/) ; colors chosen
1.1639 +; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1.1640 + resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1.1641 + resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1.1642 +
1.1643 + resz@tiYAxisString = s@long_name ; add a axis title
1.1644 + resz@txFontHeightF = 0.0195 ; change title font heights
1.1645 +
1.1646 +; Legent
1.1647 + resz@pmLegendDisplayMode = "Always" ; turn on legend
1.1648 + resz@pmLegendSide = "Top" ; Change location of
1.1649 +; resz@pmLegendParallelPosF = .45 ; move units right
1.1650 + resz@pmLegendParallelPosF = .82 ; move units right
1.1651 + resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1.1652 +
1.1653 + resz@pmLegendWidthF = 0.10 ; Change width and
1.1654 + resz@pmLegendHeightF = 0.10 ; height of legend.
1.1655 + resz@lgLabelFontHeightF = .02 ; change font height
1.1656 +; resz@lgTitleOn = True ; turn on legend title
1.1657 +; resz@lgTitleString = "Example" ; create legend title
1.1658 +; resz@lgTitleFontHeightF = .025 ; font of legend title
1.1659 + resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
1.1660 +;--------------------------------------------------------------------
1.1661 + zRes = True
1.1662 + zRes@txFontHeightF = 0.025
1.1663 +
1.1664 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1.1665 +
1.1666 + gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1.1667 +;--------------------------------------------------------------------
1.1668 +
1.1669 + plot = gsn_csm_xy (wks,ym,s,resz) ; create plot
1.1670 +
1.1671 + frame(wks) ; advance frame
1.1672 +
1.1673 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1674 + system("rm "+plot_name+"."+plot_type)
1.1675 + system("rm "+plot_name+"-1."+plot_type_new)
1.1676 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1677 +
1.1678 + clear (wks)
1.1679 +;------------------------------------------------------------------------
1.1680 + temp_name = "temp." + model_name
1.1681 + system("mkdir -p " + temp_name)
1.1682 + system("mv *.png " + temp_name)
1.1683 + system("tar cf "+ temp_name +".tar " + temp_name)
1.1684 +;------------------------------------------------------------------------
1.1685 +end
1.1686 +