1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/99.all.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1358 @@
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 set_line(lines:string,nline:integer,newlines:string)
1.14 +begin
1.15 +; add line to ascci/html file
1.16 +
1.17 + nnewlines = dimsizes(newlines)
1.18 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.19 + print("set_line: bad index, not setting anything.")
1.20 + return
1.21 + end if
1.22 + lines(nline:nline+nnewlines-1) = newlines
1.23 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.24 + nline = nline + nnewlines
1.25 + return
1.26 +end
1.27 +;****************************************************************************
1.28 +procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
1.29 + ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
1.30 + ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
1.31 + ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
1.32 + ,dx4[1]:numeric \
1.33 + )
1.34 +begin
1.35 +; Calculaee "nice" bins for binning the data in equally spaced ranges.
1.36 +; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
1.37 +; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
1.38 +
1.39 + nbins = 15 ; Number of bins to use.
1.40 +
1.41 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.42 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.43 + range = fspan(nicevals(0),nicevals(1),nvals)
1.44 +
1.45 +; Use this range information to grab all the values in a
1.46 +; particular range, and then take an average.
1.47 +
1.48 + nr = dimsizes(range)
1.49 + nx = nr-1
1.50 +; print (nx)
1.51 +
1.52 +; xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.53 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.54 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.55 + dx4 = dx/4 ; 1/4 of the range
1.56 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.57 +; yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.58 +; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.59 +; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.60 +
1.61 + do nd=0,1
1.62 +
1.63 +; See if we are doing model or observational data.
1.64 +
1.65 + if(nd.eq.0) then
1.66 + data = RAIN1_1D
1.67 + npp_data = NPP1_1D
1.68 + else
1.69 + data = RAIN2_1D
1.70 + npp_data = NPP2_1D
1.71 + end if
1.72 +
1.73 +; Loop through each range and check for values.
1.74 +
1.75 + do i=0,nr-2
1.76 + if (i.ne.(nr-2)) then
1.77 +; print("")
1.78 +; print("In range ["+range(i)+","+range(i+1)+")")
1.79 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.80 + else
1.81 +; print("")
1.82 +; print("In range ["+range(i)+",)")
1.83 + idx = ind(range(i).le.data)
1.84 + end if
1.85 +
1.86 +; Calculate average, and get min and max.
1.87 +
1.88 + if(.not.any(ismissing(idx))) then
1.89 + yvalues(nd,i) = avg(npp_data(idx))
1.90 + mn_yvalues(nd,i) = min(npp_data(idx))
1.91 + mx_yvalues(nd,i) = max(npp_data(idx))
1.92 + count = dimsizes(idx)
1.93 + else
1.94 + count = 0
1.95 + yvalues(nd,i) = yvalues@_FillValue
1.96 + mn_yvalues(nd,i) = yvalues@_FillValue
1.97 + mx_yvalues(nd,i) = yvalues@_FillValue
1.98 + end if
1.99 +
1.100 +; Print out information.
1.101 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.102 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.103 +
1.104 +; Clean up for next time in loop.
1.105 +
1.106 + delete(idx)
1.107 + end do
1.108 + delete(data)
1.109 + delete(npp_data)
1.110 + end do
1.111 +end
1.112 +;****************************************************************************
1.113 +; Main code.
1.114 +;****************************************************************************
1.115 +begin
1.116 +
1.117 + plot_type = "ps"
1.118 + plot_type_new = "png"
1.119 +
1.120 +;************************************************
1.121 +; read model data
1.122 +;************************************************
1.123 +
1.124 +;film = "i01.06cn_1798-2004_ANN_climo.nc"
1.125 +;model_name = "06cn"
1.126 +;model_grid = "T42"
1.127 +
1.128 +;film = "i01.06casa_1798-2004_ANN_climo.nc"
1.129 +;model_name = "06casa"
1.130 +;model_grid = "T42"
1.131 +
1.132 +;film = "i01.10casa_1948-2004_ANN_climo.nc"
1.133 +;model_name = "10casa"
1.134 +;model_grid = "T42"
1.135 +
1.136 + film = "i01.10cn_1948-2004_ANN_climo.nc"
1.137 + model_name = "10cn"
1.138 + model_grid = "T42"
1.139 +
1.140 + html_name = "table.html." + model_name
1.141 + html_new = html_name +".new"
1.142 + system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
1.143 + "mv -f "+html_new+" "+html_name)
1.144 +;--------------------------------------------------
1.145 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.146 + fm = addfile (dirm+film,"r")
1.147 +
1.148 + nppmod0 = fm->NPP
1.149 + rainmod0 = fm->RAIN
1.150 + xm = fm->lon
1.151 + ym = fm->lat
1.152 +
1.153 +;************************************************
1.154 +; read ob data
1.155 +;************************************************
1.156 +
1.157 +;(1) data at 81 sites
1.158 + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.159 + fil81 = "data.81.nc"
1.160 + f81 = addfile (dir81+fil81,"r")
1.161 +
1.162 + id81 = f81->SITE_ID
1.163 + npp81 = f81->TNPP_C
1.164 + rain81 = tofloat(f81->PREC_ANN)
1.165 + x81 = f81->LONG_DD
1.166 + y81 = f81->LAT_DD
1.167 +
1.168 + id81@long_name = "SITE_ID"
1.169 +
1.170 +;change longitude from (-180,180) to (0,360)
1.171 +;for model data interpolation
1.172 +
1.173 + do i= 0,dimsizes(x81)-1
1.174 + if (x81(i) .lt. 0.) then
1.175 + x81(i) = x81(i)+ 360.
1.176 + end if
1.177 + end do
1.178 +;print (x81)
1.179 +;-------------------------------------
1.180 +;(2) data at 933 sites
1.181 + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.182 + fil933 = "data.933.nc"
1.183 + f933 = addfile (dir933+fil933,"r")
1.184 +
1.185 + id933 = f933->SITE_ID
1.186 + npp933 = f933->TNPP_C
1.187 + rain933 = f933->PREC
1.188 + x933 = f933->LONG_DD
1.189 + y933 = f933->LAT_DD
1.190 +
1.191 + id933@long_name = "SITE_ID"
1.192 +
1.193 +;change longitude from (-180,180) to (0,360)
1.194 +;for model data interpolation
1.195 +
1.196 + do i= 0,dimsizes(x933)-1
1.197 + if (x933(i) .lt. 0.) then
1.198 + x933(i) = x933(i)+ 360.
1.199 + end if
1.200 + end do
1.201 +;print (x933)
1.202 +;----------------------------------------
1.203 +;(3) global data, interpolated into model grid
1.204 + dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.205 + filglobe = "Npp_"+model_grid+"_mean.nc"
1.206 + fglobe = addfile (dirglobe+filglobe,"r")
1.207 +
1.208 + nppglobe0 = fglobe->NPP
1.209 + nppglobe = nppglobe0
1.210 +
1.211 +;***********************************************************************
1.212 +; interpolate model data into ob sites
1.213 +;***********************************************************************
1.214 +
1.215 + nppmod = nppmod0(0,:,:)
1.216 + rainmod = rainmod0(0,:,:)
1.217 + delete (nppmod0)
1.218 + delete (rainmod0)
1.219 +
1.220 + nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
1.221 +
1.222 + nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
1.223 +
1.224 + rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
1.225 +
1.226 + rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
1.227 +
1.228 +;**********************************************************
1.229 +; unified units
1.230 +;**********************************************************
1.231 +; Units for these variables are:
1.232 +;
1.233 +; rain81 : mm/year
1.234 +; rainmod : mm/s
1.235 +; npp81 : g C/m^2/year
1.236 +; nppmod81: g C/m^2/s
1.237 +; nppglobe: g C/m^2/year
1.238 +;
1.239 +; We want to convert these to "m/year" and "g C/m^2/year".
1.240 +
1.241 + nsec_per_year = 60*60*24*365
1.242 +
1.243 + rain81 = rain81 / 1000.
1.244 + rainmod81 = (rainmod81/ 1000.) * nsec_per_year
1.245 + nppmod81 = nppmod81 * nsec_per_year
1.246 +
1.247 + rain933 = rain933 / 1000.
1.248 + rainmod933 = (rainmod933/ 1000.) * nsec_per_year
1.249 + nppmod933 = nppmod933 * nsec_per_year
1.250 +
1.251 + nppmod = nppmod * nsec_per_year
1.252 +
1.253 + npp81@units = "gC/m^2/yr"
1.254 + nppmod81@units = "gC/m^2/yr"
1.255 + npp933@units = "gC/m^2/yr"
1.256 + nppmod933@units = "gC/m^2/yr"
1.257 + nppmod@units = "gC/m^2/yr"
1.258 + nppglobe@units = "gC/m^2/yr"
1.259 + rain81@units = "m/yr"
1.260 + rainmod81@units = "m/yr"
1.261 + rain933@units = "m/yr"
1.262 + rainmod933@units = "m/yr"
1.263 +
1.264 + npp81@long_name = "NPP (gC/m2/year)"
1.265 + npp933@long_name = "NPP (gC/m2/year)"
1.266 + nppmod81@long_name = "NPP (gC/m2/year)"
1.267 + nppmod933@long_name = "NPP (gC/m2/year)"
1.268 + nppmod@long_name = "NPP (gC/m2/year)"
1.269 + nppglobe@long_name = "NPP (gC/m2/year)"
1.270 + rain81@long_name = "PREC (m/year)"
1.271 + rain933@long_name = "PREC (m/year)"
1.272 + rainmod81@long_name = "PREC (m/year)"
1.273 + rainmod933@long_name = "PREC (m/year)"
1.274 +
1.275 +;*******************************************************************
1.276 +;(A)-1 html table of site81 -- observed
1.277 +;*******************************************************************
1.278 + output_html = "table_site81_ob.html"
1.279 +
1.280 + header = (/"<HTML>" \
1.281 + ,"<HEAD>" \
1.282 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.283 + ,"</HEAD>" \
1.284 + ,"<H1>Observation at 81 sites</H1>" \
1.285 + /)
1.286 + footer = "</HTML>"
1.287 +
1.288 + table_header = (/ \
1.289 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.290 + ,"<tr>" \
1.291 + ," <th bgcolor=DDDDDD >Site ID</th>" \
1.292 + ," <th bgcolor=DDDDDD >Latitude</th>" \
1.293 + ," <th bgcolor=DDDDDD >Longitude</th>" \
1.294 + ," <th bgcolor=DDDDDD >NPP(gC/m2.year)</th>" \
1.295 + ," <th bgcolor=DDDDDD >RAIN(m/year)</th>" \
1.296 + ,"</tr>" \
1.297 + /)
1.298 + table_footer = "</table>"
1.299 + row_header = "<tr>"
1.300 + row_footer = "</tr>"
1.301 +
1.302 + lines = new(50000,string)
1.303 + nline = 0
1.304 +
1.305 + set_line(lines,nline,header)
1.306 + set_line(lines,nline,table_header)
1.307 +;-----------------------------------------------
1.308 +;row of table
1.309 +
1.310 + nrow = dimsizes(id81)
1.311 + do n = 0,nrow-1
1.312 + set_line(lines,nline,row_header)
1.313 +
1.314 + txt1 = sprintf("%5.0f", id81(n))
1.315 + txt2 = sprintf("%5.2f", y81(n))
1.316 + txt3 = sprintf("%5.2f", x81(n))
1.317 + txt4 = sprintf("%5.2f", npp81(n))
1.318 + txt5 = sprintf("%5.2f", rain81(n))
1.319 +
1.320 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.321 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.322 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.323 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.324 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.325 +
1.326 + set_line(lines,nline,row_footer)
1.327 + end do
1.328 +;-----------------------------------------------
1.329 + set_line(lines,nline,table_footer)
1.330 + set_line(lines,nline,footer)
1.331 +
1.332 +; Now write to an HTML file.
1.333 + idx = ind(.not.ismissing(lines))
1.334 + if(.not.any(ismissing(idx))) then
1.335 + asciiwrite(output_html,lines(idx))
1.336 + else
1.337 + print ("error?")
1.338 + end if
1.339 +;*******************************************************************
1.340 +;(A)-2 html table of site933 -- observed
1.341 +;*******************************************************************
1.342 + output_html = "table_site933_ob.html"
1.343 +
1.344 + header = (/"<HTML>" \
1.345 + ,"<HEAD>" \
1.346 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.347 + ,"</HEAD>" \
1.348 + ,"<H1>Observation at 933 sites</H1>" \
1.349 + /)
1.350 + footer = "</HTML>"
1.351 +
1.352 + table_header = (/ \
1.353 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.354 + ,"<tr>" \
1.355 + ," <th bgcolor=DDDDDD >Site ID</th>" \
1.356 + ," <th bgcolor=DDDDDD >Latitude</th>" \
1.357 + ," <th bgcolor=DDDDDD >Longitude</th>" \
1.358 + ," <th bgcolor=DDDDDD >NPP(gC/m2.year)</th>" \
1.359 + ," <th bgcolor=DDDDDD >RAIN(m/year)</th>" \
1.360 + ,"</tr>" \
1.361 + /)
1.362 + table_footer = "</table>"
1.363 + row_header = "<tr>"
1.364 + row_footer = "</tr>"
1.365 +
1.366 + delete (lines)
1.367 + lines = new(50000,string)
1.368 + nline = 0
1.369 +
1.370 + set_line(lines,nline,header)
1.371 + set_line(lines,nline,table_header)
1.372 +;-----------------------------------------------
1.373 +;row of table
1.374 +
1.375 + nrow = dimsizes(id933)
1.376 + do n = 0,nrow-1
1.377 + set_line(lines,nline,row_header)
1.378 +
1.379 + txt1 = sprintf("%5.0f", id933(n))
1.380 + txt2 = sprintf("%5.2f", y933(n))
1.381 + txt3 = sprintf("%5.2f", x933(n))
1.382 + txt4 = sprintf("%5.2f", npp933(n))
1.383 + txt5 = sprintf("%5.2f", rain933(n))
1.384 +
1.385 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.386 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.387 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.388 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.389 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.390 +
1.391 + set_line(lines,nline,row_footer)
1.392 + end do
1.393 +;-----------------------------------------------
1.394 + set_line(lines,nline,table_footer)
1.395 + set_line(lines,nline,footer)
1.396 +
1.397 +; Now write to an HTML file.
1.398 + delete (idx)
1.399 + idx = ind(.not.ismissing(lines))
1.400 + if(.not.any(ismissing(idx))) then
1.401 + asciiwrite(output_html,lines(idx))
1.402 + else
1.403 + print ("error?")
1.404 + end if
1.405 +;*******************************************************************
1.406 +;(A)-3 html table of site81 -- model vs ob
1.407 +;*******************************************************************
1.408 + output_html = "table_site81_model_vs_ob.html"
1.409 +
1.410 + header_text = "<H1>Model "+model_name+" vs Observed at 81 sites</H1>"
1.411 +
1.412 + header = (/"<HTML>" \
1.413 + ,"<HEAD>" \
1.414 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.415 + ,"</HEAD>" \
1.416 + ,header_text \
1.417 + /)
1.418 + footer = "</HTML>"
1.419 +
1.420 + delete (table_header)
1.421 + table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>"
1.422 + table_header = (/ \
1.423 + "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
1.424 + ,"<tr>" \
1.425 + ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
1.426 + ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
1.427 + ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
1.428 + ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \
1.429 + ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \
1.430 + ,"</tr>" \
1.431 + ,"<tr>" \
1.432 + ," <th bgcolor=DDDDDD >observed</th>" \
1.433 + , table_header_text \
1.434 + ," <th bgcolor=DDDDDD >observed</th>" \
1.435 + , table_header_text \
1.436 + ,"</tr>" \
1.437 + /)
1.438 + table_footer = "</table>"
1.439 + row_header = "<tr>"
1.440 + row_footer = "</tr>"
1.441 +
1.442 + delete (lines)
1.443 + lines = new(50000,string)
1.444 + nline = 0
1.445 +
1.446 + set_line(lines,nline,header)
1.447 + set_line(lines,nline,table_header)
1.448 +;-----------------------------------------------
1.449 +;row of table
1.450 +
1.451 + nrow = dimsizes(id81)
1.452 + do n = 0,nrow-1
1.453 + set_line(lines,nline,row_header)
1.454 +
1.455 + txt1 = sprintf("%5.0f", id81(n))
1.456 + txt2 = sprintf("%5.2f", y81(n))
1.457 + txt3 = sprintf("%5.2f", x81(n))
1.458 + txt4 = sprintf("%5.2f", npp81(n))
1.459 + txt5 = sprintf("%5.2f", nppmod81(n))
1.460 + txt6 = sprintf("%5.2f", rain81(n))
1.461 + txt7 = sprintf("%5.2f", rainmod81(n))
1.462 +
1.463 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.464 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.465 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.466 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.467 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.468 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.469 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.470 +
1.471 + set_line(lines,nline,row_footer)
1.472 + end do
1.473 +;-----------------------------------------------
1.474 + set_line(lines,nline,table_footer)
1.475 + set_line(lines,nline,footer)
1.476 +
1.477 +; Now write to an HTML file.
1.478 + delete (idx)
1.479 + idx = ind(.not.ismissing(lines))
1.480 + if(.not.any(ismissing(idx))) then
1.481 + asciiwrite(output_html,lines(idx))
1.482 + else
1.483 + print ("error?")
1.484 + end if
1.485 +
1.486 +;*******************************************************************
1.487 +;(A)-4 html table of site933 -- model vs ob
1.488 +;*******************************************************************
1.489 + output_html = "table_site933_model_vs_ob.html"
1.490 +
1.491 + header_text = "<H1>Model "+model_name+" vs Observed at 933 sites</H1>"
1.492 +
1.493 + header = (/"<HTML>" \
1.494 + ,"<HEAD>" \
1.495 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.496 + ,"</HEAD>" \
1.497 + ,header_text \
1.498 + /)
1.499 + footer = "</HTML>"
1.500 +
1.501 + delete (table_header)
1.502 + table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>"
1.503 + table_header = (/ \
1.504 + "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
1.505 + ,"<tr>" \
1.506 + ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
1.507 + ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
1.508 + ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
1.509 + ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \
1.510 + ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \
1.511 + ,"</tr>" \
1.512 + ,"<tr>" \
1.513 + ," <th bgcolor=DDDDDD >observed</th>" \
1.514 + , table_header_text \
1.515 + ," <th bgcolor=DDDDDD >observed</th>" \
1.516 + , table_header_text \
1.517 + ,"</tr>" \
1.518 + /)
1.519 + table_footer = "</table>"
1.520 + row_header = "<tr>"
1.521 + row_footer = "</tr>"
1.522 +
1.523 + delete (lines)
1.524 + lines = new(50000,string)
1.525 + nline = 0
1.526 +
1.527 + set_line(lines,nline,header)
1.528 + set_line(lines,nline,table_header)
1.529 +;-----------------------------------------------
1.530 +;row of table
1.531 +
1.532 + nrow = dimsizes(id933)
1.533 + do n = 0,nrow-1
1.534 + set_line(lines,nline,row_header)
1.535 +
1.536 + txt1 = sprintf("%5.0f", id933(n))
1.537 + txt2 = sprintf("%5.2f", y933(n))
1.538 + txt3 = sprintf("%5.2f", x933(n))
1.539 + txt4 = sprintf("%5.2f", npp933(n))
1.540 + txt5 = sprintf("%5.2f", nppmod933(n))
1.541 + txt6 = sprintf("%5.2f", rain933(n))
1.542 + txt7 = sprintf("%5.2f", rainmod933(n))
1.543 +
1.544 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.545 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.546 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.547 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.548 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.549 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.550 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.551 +
1.552 + set_line(lines,nline,row_footer)
1.553 + end do
1.554 +;-----------------------------------------------
1.555 + set_line(lines,nline,table_footer)
1.556 + set_line(lines,nline,footer)
1.557 +
1.558 +; Now write to an HTML file.
1.559 + delete (idx)
1.560 + idx = ind(.not.ismissing(lines))
1.561 + if(.not.any(ismissing(idx))) then
1.562 + asciiwrite(output_html,lines(idx))
1.563 + else
1.564 + print ("error?")
1.565 + end if
1.566 +;***************************************************************************
1.567 +;(A)-5 scatter plot, model vs ob 81
1.568 +;***************************************************************************
1.569 +
1.570 + plot_name = "scatter_model_vs_ob_81"
1.571 + title = model_name +" vs ob 81 sites"
1.572 +
1.573 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.574 + res = True ; plot mods desired
1.575 + res@tiMainString = title ; add title
1.576 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.577 + res@xyMarkers = 16 ; choose type of marker
1.578 + res@xyMarkerColor = "red" ; Marker color
1.579 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.580 + res@tmLabelAutoStride = True ; nice tick mark labels
1.581 +
1.582 + res@gsnDraw = False
1.583 + res@gsnFrame = False ; don't advance frame yet
1.584 +;-------------------------------
1.585 +;compute correlation coef. and M
1.586 + ccr81 = esccr(nppmod81,npp81,0)
1.587 +;print (ccr81)
1.588 +
1.589 +;new eq
1.590 + bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
1.591 + M81 = (1. - (bias/dimsizes(y81)))*5.
1.592 +
1.593 + M_npp_S81 = sprintf("%.2f", M81)
1.594 + system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \
1.595 + "mv -f "+html_new+" "+html_name)
1.596 + print (M_npp_S81)
1.597 +
1.598 + tRes = True
1.599 + tRes@txFontHeightF = 0.025
1.600 +
1.601 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.602 +
1.603 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.604 +;--------------------------------
1.605 + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
1.606 +;-------------------------------
1.607 +; add polyline
1.608 +
1.609 + dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
1.610 +;-------------------------------
1.611 + draw (plot)
1.612 + frame(wks)
1.613 + clear (wks)
1.614 + delete (dum)
1.615 +
1.616 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.617 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.618 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.619 + "rm "+plot_name+"."+plot_type)
1.620 +;***************************************************************************
1.621 +;(A)-6 scatter plot, model vs ob 933
1.622 +;***************************************************************************
1.623 +
1.624 + plot_name = "scatter_model_vs_ob_933"
1.625 + title = model_name +" vs ob 933 sites"
1.626 +
1.627 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.628 + res = True ; plot mods desired
1.629 + res@tiMainString = title ; add title
1.630 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.631 + res@xyMarkers = 16 ; choose type of marker
1.632 + res@xyMarkerColor = "red" ; Marker color
1.633 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.634 + res@tmLabelAutoStride = True ; nice tick mark labels
1.635 +
1.636 + res@gsnDraw = False
1.637 + res@gsnFrame = False ; don't advance frame yet
1.638 +;-------------------------------
1.639 +;compute correlation coef. and M
1.640 + ccr933 = esccr(nppmod933,npp933,0)
1.641 +;print (ccr933)
1.642 +
1.643 +;new eq
1.644 + bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
1.645 + M933 = (1. - (bias/dimsizes(y933)))*5.
1.646 +
1.647 + M_npp_S933 = sprintf("%.2f", M933)
1.648 + system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \
1.649 + "mv -f "+html_new+" "+html_name)
1.650 + print (M_npp_S933)
1.651 +
1.652 + tRes = True
1.653 + tRes@txFontHeightF = 0.025
1.654 +
1.655 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.656 +
1.657 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.658 +;--------------------------------
1.659 + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
1.660 +;-------------------------------
1.661 +; add polyline
1.662 +
1.663 + dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
1.664 +;-------------------------------
1.665 + draw (plot)
1.666 + frame(wks)
1.667 + clear (wks)
1.668 + delete (dum)
1.669 +
1.670 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.671 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.672 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.673 + "rm "+plot_name+"."+plot_type)
1.674 +;**************************************************************************
1.675 +;(B) histogram 81
1.676 +;**************************************************************************
1.677 +;get data
1.678 + RAIN1_1D = ndtooned(rain81)
1.679 + RAIN2_1D = ndtooned(rainmod81)
1.680 + NPP1_1D = ndtooned(npp81)
1.681 + NPP2_1D = ndtooned(nppmod81)
1.682 +
1.683 +; number of bin
1.684 + nx = 8
1.685 +
1.686 + xvalues = new((/2,nx/),float)
1.687 + yvalues = new((/2,nx/),float)
1.688 + mn_yvalues = new((/2,nx/),float)
1.689 + mx_yvalues = new((/2,nx/),float)
1.690 + dx4 = new((/1/),float)
1.691 +
1.692 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.693 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.694 +
1.695 +;----------------------------------------
1.696 +;compute correlation coeff and M score
1.697 +
1.698 + u = yvalues(0,:)
1.699 + v = yvalues(1,:)
1.700 +
1.701 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.702 + uu = u(good)
1.703 + vv = v(good)
1.704 +
1.705 + ccr81h = esccr(uu,vv,0)
1.706 +;print (ccr81h)
1.707 +
1.708 +;new eq
1.709 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.710 + M81h = (1.- (bias/dimsizes(uu)))*5.
1.711 +
1.712 + M_npp_H81 = sprintf("%.2f", M81h)
1.713 + system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \
1.714 + "mv -f "+html_new+" "+html_name)
1.715 + print (M_npp_H81)
1.716 +
1.717 + delete (u)
1.718 + delete (v)
1.719 + delete (uu)
1.720 + delete (vv)
1.721 +;----------------------------------------------------------------------
1.722 +; histogram res
1.723 + resh = True
1.724 + resh@gsnMaximize = True
1.725 + resh@gsnDraw = False
1.726 + resh@gsnFrame = False
1.727 + resh@xyMarkLineMode = "Markers"
1.728 + resh@xyMarkerSizeF = 0.014
1.729 + resh@xyMarker = 16
1.730 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.731 + resh@trYMinF = min(mn_yvalues) - 10.
1.732 + resh@trYMaxF = max(mx_yvalues) + 10.
1.733 +
1.734 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.735 + resh@tiXAxisString = "Precipitation (m/year)"
1.736 +
1.737 + max_bar = new((/2,nx/),graphic)
1.738 + min_bar = new((/2,nx/),graphic)
1.739 + max_cap = new((/2,nx/),graphic)
1.740 + min_cap = new((/2,nx/),graphic)
1.741 +
1.742 + lnres = True
1.743 + line_colors = (/"brown","blue"/)
1.744 +
1.745 +;**************************************************************************
1.746 +;(B)-1 histogram plot, ob 81 site
1.747 +;**************************************************************************
1.748 +
1.749 + plot_name = "histogram_ob_81"
1.750 + title = "Observed 81 site"
1.751 + resh@tiMainString = title
1.752 +
1.753 + wks = gsn_open_wks (plot_type,plot_name)
1.754 +
1.755 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.756 +
1.757 +;-------------------------------
1.758 +;Attach the vertical bar and the horizontal cap line
1.759 +
1.760 + do nd=0,0
1.761 + lnres@gsLineColor = line_colors(nd)
1.762 + do i=0,nx-1
1.763 +
1.764 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.765 + .not.ismissing(mx_yvalues(nd,i))) then
1.766 +;
1.767 +; Attach the vertical bar, both above and below the marker.
1.768 +;
1.769 + x1 = xvalues(nd,i)
1.770 + y1 = yvalues(nd,i)
1.771 + y2 = mn_yvalues(nd,i)
1.772 + plx = (/x1,x1/)
1.773 + ply = (/y1,y2/)
1.774 + min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.775 +
1.776 + y2 = mx_yvalues(nd,i)
1.777 + plx = (/x1,x1/)
1.778 + ply = (/y1,y2/)
1.779 + max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.780 +;
1.781 +; Attach the horizontal cap line, both above and below the marker.
1.782 +;
1.783 + x1 = xvalues(nd,i) - dx4
1.784 + x2 = xvalues(nd,i) + dx4
1.785 + y1 = mn_yvalues(nd,i)
1.786 + plx = (/x1,x2/)
1.787 + ply = (/y1,y1/)
1.788 + min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.789 +
1.790 + y1 = mx_yvalues(nd,i)
1.791 + plx = (/x1,x2/)
1.792 + ply = (/y1,y1/)
1.793 + max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.794 + end if
1.795 + end do
1.796 + end do
1.797 +
1.798 + draw(xy)
1.799 + frame(wks)
1.800 + clear (wks)
1.801 +
1.802 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.803 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.804 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.805 + "rm "+plot_name+"."+plot_type)
1.806 +;****************************************************************************
1.807 +;(B)-2 histogram plot, model vs ob 81 site
1.808 +;****************************************************************************
1.809 +
1.810 + plot_name = "histogram_model_vs_ob_81"
1.811 + title = model_name+ " vs Observed 81 site"
1.812 + resh@tiMainString = title
1.813 +
1.814 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.815 +
1.816 +;-----------------------------
1.817 +; Add a boxed legend using the more simple method, which won't have
1.818 +; vertical lines going through the markers.
1.819 +
1.820 + resh@pmLegendDisplayMode = "Always"
1.821 +; resh@pmLegendWidthF = 0.1
1.822 + resh@pmLegendWidthF = 0.08
1.823 + resh@pmLegendHeightF = 0.05
1.824 + resh@pmLegendOrthogonalPosF = -1.17
1.825 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.826 +; resh@pmLegendParallelPosF = 0.18
1.827 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.828 +
1.829 +; resh@lgPerimOn = False
1.830 + resh@lgLabelFontHeightF = 0.015
1.831 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.832 +;-----------------------------
1.833 + tRes = True
1.834 + tRes@txFontHeightF = 0.025
1.835 +
1.836 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1.837 +
1.838 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.839 +
1.840 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.841 +;-------------------------------
1.842 +;Attach the vertical bar and the horizontal cap line
1.843 +
1.844 + do nd=0,1
1.845 + lnres@gsLineColor = line_colors(nd)
1.846 + do i=0,nx-1
1.847 +
1.848 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.849 + .not.ismissing(mx_yvalues(nd,i))) then
1.850 +;
1.851 +; Attach the vertical bar, both above and below the marker.
1.852 +;
1.853 + x1 = xvalues(nd,i)
1.854 + y1 = yvalues(nd,i)
1.855 + y2 = mn_yvalues(nd,i)
1.856 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.857 +
1.858 + y2 = mx_yvalues(nd,i)
1.859 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.860 +;
1.861 +; Attach the horizontal cap line, both above and below the marker.
1.862 +;
1.863 + x1 = xvalues(nd,i) - dx4
1.864 + x2 = xvalues(nd,i) + dx4
1.865 + y1 = mn_yvalues(nd,i)
1.866 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.867 +
1.868 + y1 = mx_yvalues(nd,i)
1.869 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.870 + end if
1.871 + end do
1.872 + end do
1.873 +
1.874 + draw(xy)
1.875 + frame(wks)
1.876 + clear(wks)
1.877 +
1.878 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.879 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.880 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.881 + "rm "+plot_name+"."+plot_type)
1.882 +
1.883 + delete (RAIN1_1D)
1.884 + delete (RAIN2_1D)
1.885 + delete (NPP1_1D)
1.886 + delete (NPP2_1D)
1.887 +;delete (range)
1.888 + delete (xvalues)
1.889 + delete (yvalues)
1.890 + delete (mn_yvalues)
1.891 + delete (mx_yvalues)
1.892 + delete (good)
1.893 + delete (max_bar)
1.894 + delete (min_bar)
1.895 + delete (max_cap)
1.896 + delete (min_cap)
1.897 +
1.898 +;**************************************************************************
1.899 +;(B) histogram 933
1.900 +;**************************************************************************
1.901 +
1.902 +; get data
1.903 + RAIN1_1D = ndtooned(rain933)
1.904 + RAIN2_1D = ndtooned(rainmod933)
1.905 + NPP1_1D = ndtooned(npp933)
1.906 + NPP2_1D = ndtooned(nppmod933)
1.907 +
1.908 +; number of bin
1.909 + nx = 10
1.910 +
1.911 + xvalues = new((/2,nx/),float)
1.912 + yvalues = new((/2,nx/),float)
1.913 + mn_yvalues = new((/2,nx/),float)
1.914 + mx_yvalues = new((/2,nx/),float)
1.915 + dx4 = new((/1/),float)
1.916 +
1.917 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.918 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.919 +
1.920 +;----------------------------------------
1.921 +;compute correlation coeff and M score
1.922 +
1.923 + u = yvalues(0,:)
1.924 + v = yvalues(1,:)
1.925 +
1.926 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.927 + uu = u(good)
1.928 + vv = v(good)
1.929 +
1.930 + ccr933h = esccr(uu,vv,0)
1.931 +;print (ccr933h)
1.932 +
1.933 +;new eq
1.934 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.935 + M933h = (1.- (bias/dimsizes(uu)))*5.
1.936 +
1.937 + M_npp_H933 = sprintf("%.2f", M933h)
1.938 + system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \
1.939 + "mv -f "+html_new+" "+html_name)
1.940 + print (M_npp_H933)
1.941 +
1.942 + delete (u)
1.943 + delete (v)
1.944 + delete (uu)
1.945 + delete (vv)
1.946 +;----------------------------------------------------------------------
1.947 +; histogram res
1.948 + delete (resh)
1.949 + resh = True
1.950 + resh@gsnMaximize = True
1.951 + resh@gsnDraw = False
1.952 + resh@gsnFrame = False
1.953 + resh@xyMarkLineMode = "Markers"
1.954 + resh@xyMarkerSizeF = 0.014
1.955 + resh@xyMarker = 16
1.956 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.957 + resh@trYMinF = min(mn_yvalues) - 10.
1.958 + resh@trYMaxF = max(mx_yvalues) + 10.
1.959 +
1.960 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.961 + resh@tiXAxisString = "Precipitation (m/year)"
1.962 +
1.963 + max_bar = new((/2,nx/),graphic)
1.964 + min_bar = new((/2,nx/),graphic)
1.965 + max_cap = new((/2,nx/),graphic)
1.966 + min_cap = new((/2,nx/),graphic)
1.967 +
1.968 + lnres = True
1.969 + line_colors = (/"brown","blue"/)
1.970 +
1.971 +;**************************************************************************
1.972 +;(B)-3 histogram plot, ob 933 site
1.973 +;**************************************************************************
1.974 +
1.975 + plot_name = "histogram_ob_933"
1.976 + title = "Observed 933 site"
1.977 + resh@tiMainString = title
1.978 +
1.979 + wks = gsn_open_wks (plot_type,plot_name)
1.980 +
1.981 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.982 +
1.983 +;-------------------------------
1.984 +;Attach the vertical bar and the horizontal cap line
1.985 +
1.986 + do nd=0,0
1.987 + lnres@gsLineColor = line_colors(nd)
1.988 + do i=0,nx-1
1.989 +
1.990 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.991 + .not.ismissing(mx_yvalues(nd,i))) then
1.992 +
1.993 +; Attach the vertical bar, both above and below the marker.
1.994 +
1.995 + x1 = xvalues(nd,i)
1.996 + y1 = yvalues(nd,i)
1.997 + y2 = mn_yvalues(nd,i)
1.998 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.999 +
1.1000 + y2 = mx_yvalues(nd,i)
1.1001 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1002 +
1.1003 +; Attach the horizontal cap line, both above and below the marker.
1.1004 +
1.1005 + x1 = xvalues(nd,i) - dx4
1.1006 + x2 = xvalues(nd,i) + dx4
1.1007 + y1 = mn_yvalues(nd,i)
1.1008 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1009 +
1.1010 + y1 = mx_yvalues(nd,i)
1.1011 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1012 + end if
1.1013 + end do
1.1014 + end do
1.1015 +
1.1016 + draw(xy)
1.1017 + frame(wks)
1.1018 + delete (xy)
1.1019 + clear (wks)
1.1020 +
1.1021 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1022 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1023 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1024 + "rm "+plot_name+"."+plot_type)
1.1025 +;****************************************************************************
1.1026 +;(B)-4 histogram plot, model vs ob 933 site
1.1027 +;****************************************************************************
1.1028 +
1.1029 + plot_name = "histogram_model_vs_ob_933"
1.1030 + title = model_name+ " vs Observed 933 site"
1.1031 + resh@tiMainString = title
1.1032 +
1.1033 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1034 +
1.1035 +;-----------------------------
1.1036 +; Add a boxed legend using the more simple method, which won't have
1.1037 +; vertical lines going through the markers.
1.1038 +
1.1039 + resh@pmLegendDisplayMode = "Always"
1.1040 +; resh@pmLegendWidthF = 0.1
1.1041 + resh@pmLegendWidthF = 0.08
1.1042 + resh@pmLegendHeightF = 0.05
1.1043 + resh@pmLegendOrthogonalPosF = -1.17
1.1044 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1045 +; resh@pmLegendParallelPosF = 0.18
1.1046 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.1047 +
1.1048 +; resh@lgPerimOn = False
1.1049 + resh@lgLabelFontHeightF = 0.015
1.1050 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.1051 +;-----------------------------
1.1052 + tRes = True
1.1053 + tRes@txFontHeightF = 0.025
1.1054 +
1.1055 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1.1056 +
1.1057 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1058 +
1.1059 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.1060 +;-------------------------------
1.1061 +;Attach the vertical bar and the horizontal cap line
1.1062 +
1.1063 + do nd=0,1
1.1064 + lnres@gsLineColor = line_colors(nd)
1.1065 + do i=0,nx-1
1.1066 +
1.1067 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1068 + .not.ismissing(mx_yvalues(nd,i))) then
1.1069 +;
1.1070 +; Attach the vertical bar, both above and below the marker.
1.1071 +;
1.1072 + x1 = xvalues(nd,i)
1.1073 + y1 = yvalues(nd,i)
1.1074 + y2 = mn_yvalues(nd,i)
1.1075 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1076 +
1.1077 + y2 = mx_yvalues(nd,i)
1.1078 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1079 +;
1.1080 +; Attach the horizontal cap line, both above and below the marker.
1.1081 +;
1.1082 + x1 = xvalues(nd,i) - dx4
1.1083 + x2 = xvalues(nd,i) + dx4
1.1084 + y1 = mn_yvalues(nd,i)
1.1085 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1086 +
1.1087 + y1 = mx_yvalues(nd,i)
1.1088 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1089 + end if
1.1090 + end do
1.1091 + end do
1.1092 +
1.1093 + draw(xy)
1.1094 + frame(wks)
1.1095 + delete(xy)
1.1096 + clear(wks)
1.1097 +
1.1098 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1099 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1100 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1101 + "rm "+plot_name+"."+plot_type)
1.1102 +;***************************************************************************
1.1103 +;(C) global contour
1.1104 +;***************************************************************************
1.1105 +
1.1106 +;res
1.1107 + resg = True ; Use plot options
1.1108 + resg@cnFillOn = True ; Turn on color fill
1.1109 + resg@gsnSpreadColors = True ; use full colormap
1.1110 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1111 +; resg@lbLabelAutoStride = True
1.1112 + resg@cnLinesOn = False ; Turn off contourn lines
1.1113 + resg@mpFillOn = False ; Turn off map fill
1.1114 +
1.1115 + resg@gsnSpreadColors = True ; use full colormap
1.1116 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1117 + resg@cnMinLevelValF = 0. ; Min level
1.1118 + resg@cnMaxLevelValF = 2200. ; Max level
1.1119 + resg@cnLevelSpacingF = 200. ; interval
1.1120 +
1.1121 +;****************************************************************************
1.1122 +;(C)-1 global contour plot, ob
1.1123 +;****************************************************************************
1.1124 +
1.1125 + delta = 0.00000000001
1.1126 + nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1.1127 +
1.1128 + plot_name = "global_ob"
1.1129 + title = "Observed MODIS MOD 17"
1.1130 + resg@tiMainString = title
1.1131 +
1.1132 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1133 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1134 +
1.1135 + plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1136 +
1.1137 + frame(wks)
1.1138 + clear (wks)
1.1139 +
1.1140 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1141 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1142 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1143 + "rm "+plot_name+"."+plot_type)
1.1144 +;****************************************************************************
1.1145 +;(C)-2 global contour plot, model
1.1146 +;****************************************************************************
1.1147 +
1.1148 + plot_name = "global_model"
1.1149 + title = "Model "+ model_name
1.1150 + resg@tiMainString = title
1.1151 +
1.1152 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1153 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1154 +
1.1155 + plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1156 +
1.1157 + frame(wks)
1.1158 + clear (wks)
1.1159 +
1.1160 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1161 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1162 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1163 + "rm "+plot_name+"."+plot_type)
1.1164 +;****************************************************************************
1.1165 +;(C)-3 global contour plot, model vs ob
1.1166 +;****************************************************************************
1.1167 +
1.1168 + plot_name = "global_model_vs_ob"
1.1169 +
1.1170 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1171 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1172 +
1.1173 + delete (plot)
1.1174 + plot=new(3,graphic) ; create graphic array
1.1175 +
1.1176 + resg@gsnFrame = False ; Do not draw plot
1.1177 + resg@gsnDraw = False ; Do not advance frame
1.1178 +
1.1179 +; compute correlation coef and M score
1.1180 +
1.1181 + uu1 = ndtooned(nppmod)
1.1182 + vv1 = ndtooned(nppglobe)
1.1183 +
1.1184 + delete (good)
1.1185 + good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1.1186 +
1.1187 + ug = uu1(good)
1.1188 + vg = vv1(good)
1.1189 +
1.1190 + ccrG = esccr(ug,vg,0)
1.1191 +; print (ccrG)
1.1192 +
1.1193 + MG = (ccrG*ccrG)* 5.0
1.1194 +
1.1195 + M_npp_G = sprintf("%.2f", MG)
1.1196 + system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \
1.1197 + "mv -f "+html_new+" "+html_name)
1.1198 + print (M_npp_G)
1.1199 +
1.1200 +; plot correlation coef
1.1201 +
1.1202 + gRes = True
1.1203 + gRes@txFontHeightF = 0.02
1.1204 + gRes@txAngleF = 90
1.1205 +
1.1206 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.1207 +
1.1208 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1209 +;--------------------------------------------------------------------
1.1210 +;(a) ob
1.1211 +
1.1212 + title = "Observed MODIS MOD 17"
1.1213 + resg@tiMainString = title
1.1214 +
1.1215 + plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1216 +
1.1217 +;(b) model
1.1218 +
1.1219 + title = "Model "+ model_name
1.1220 + resg@tiMainString = title
1.1221 +
1.1222 + plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1223 +
1.1224 +;(c) model-ob
1.1225 +
1.1226 + zz = nppmod
1.1227 + zz = nppmod - nppglobe
1.1228 + title = "Model_"+model_name+" - Observed"
1.1229 +
1.1230 + resg@cnMinLevelValF = -500 ; Min level
1.1231 + resg@cnMaxLevelValF = 500. ; Max level
1.1232 + resg@cnLevelSpacingF = 50. ; interval
1.1233 + resg@tiMainString = title
1.1234 +
1.1235 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1236 +
1.1237 + pres = True ; panel plot mods desired
1.1238 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1239 + ; indiv. plots in panel
1.1240 + pres@gsnMaximize = True ; fill the page
1.1241 +
1.1242 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1243 +
1.1244 + frame (wks)
1.1245 + clear (wks)
1.1246 + delete (plot)
1.1247 +
1.1248 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1249 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1250 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1251 + "rm "+plot_name+"."+plot_type)
1.1252 +;***************************************************************************
1.1253 +;(D)-1 zonal line plot, ob
1.1254 +;***************************************************************************
1.1255 +
1.1256 + vv = zonalAve(nppglobe)
1.1257 + vv@long_name = nppglobe@long_name
1.1258 +
1.1259 + plot_name = "zonal_ob"
1.1260 + title = "Observed MODIS MOD 17"
1.1261 +
1.1262 + resz = True
1.1263 + resz@tiMainString = title
1.1264 +
1.1265 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1266 +
1.1267 + plot = gsn_csm_xy (wks,ym,vv,resz)
1.1268 +
1.1269 + frame(wks)
1.1270 + clear (wks)
1.1271 +
1.1272 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1273 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1274 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1275 + "rm "+plot_name+"."+plot_type)
1.1276 +;****************************************************************************
1.1277 +;(D)-2 zonal line plot, model vs ob
1.1278 +;****************************************************************************
1.1279 +
1.1280 + s = new ((/2,dimsizes(ym)/), float)
1.1281 +
1.1282 + s(0,:) = zonalAve(nppglobe)
1.1283 + s(1,:) = zonalAve(nppmod)
1.1284 +
1.1285 + s@long_name = nppglobe@long_name
1.1286 +;-------------------------------------------
1.1287 +; compute correlation coef and M score
1.1288 +
1.1289 + ccrZ = esccr(s(1,:), s(0,:),0)
1.1290 +; print (ccrZ)
1.1291 +
1.1292 + MZ = (ccrZ*ccrZ)* 5.0
1.1293 +
1.1294 + M_npp_Z = sprintf("%.2f", MZ)
1.1295 + system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \
1.1296 + "mv -f "+html_new+" "+html_name)
1.1297 + print (M_npp_Z)
1.1298 +;-------------------------------------------
1.1299 + plot_name = "zonal_model_vs_ob"
1.1300 + title = "Zonal Average"
1.1301 + resz@tiMainString = title
1.1302 +
1.1303 + wks = gsn_open_wks (plot_type,plot_name)
1.1304 +
1.1305 +; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1.1306 +; resz@vpWidthF = 0.7
1.1307 +
1.1308 + resz@xyMonoLineColor = "False" ; want colored lines
1.1309 + resz@xyLineColors = (/"black","red"/) ; colors chosen
1.1310 +; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1.1311 + resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1.1312 + resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1.1313 +
1.1314 + resz@tiYAxisString = s@long_name ; add a axis title
1.1315 + resz@txFontHeightF = 0.0195 ; change title font heights
1.1316 +
1.1317 +; Legent
1.1318 + resz@pmLegendDisplayMode = "Always" ; turn on legend
1.1319 + resz@pmLegendSide = "Top" ; Change location of
1.1320 +; resz@pmLegendParallelPosF = .45 ; move units right
1.1321 + resz@pmLegendParallelPosF = .82 ; move units right
1.1322 + resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1.1323 +
1.1324 + resz@pmLegendWidthF = 0.10 ; Change width and
1.1325 + resz@pmLegendHeightF = 0.10 ; height of legend.
1.1326 + resz@lgLabelFontHeightF = .02 ; change font height
1.1327 +; resz@lgTitleOn = True ; turn on legend title
1.1328 +; resz@lgTitleString = "Example" ; create legend title
1.1329 +; resz@lgTitleFontHeightF = .025 ; font of legend title
1.1330 + resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
1.1331 +;--------------------------------------------------------------------
1.1332 + zRes = True
1.1333 + zRes@txFontHeightF = 0.025
1.1334 +
1.1335 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1.1336 +
1.1337 + gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1.1338 +;--------------------------------------------------------------------
1.1339 +
1.1340 + plot = gsn_csm_xy (wks,ym,s,resz)
1.1341 +
1.1342 + frame(wks)
1.1343 + clear (wks)
1.1344 +
1.1345 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1346 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.1347 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.1348 + "rm "+plot_name+"."+plot_type)
1.1349 +;***************************************************************************
1.1350 +; tar up output plots
1.1351 +;***************************************************************************
1.1352 +
1.1353 + temp_name = "npp." + model_name
1.1354 + system("mkdir -p " + temp_name+";"+ \
1.1355 + "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
1.1356 + "mv table_*.html " + temp_name +";"+ \
1.1357 + "mv *.png " + temp_name +";"+ \
1.1358 + "tar cf "+ temp_name +".tar " + temp_name)
1.1359 +;***************************************************************************
1.1360 +end
1.1361 +