1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/beta/05.biome.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,487 @@
1.4 +;********************************************************
1.5 +; required command line input parameters:
1.6 +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
1.7 +;
1.8 +; histogram normalized by rain and compute correleration
1.9 +;**************************************************************
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.13 +;**************************************************************
1.14 +procedure set_line(lines:string,nline:integer,newlines:string)
1.15 +begin
1.16 +; add line to ascci/html file
1.17 +
1.18 + nnewlines = dimsizes(newlines)
1.19 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.20 + print("set_line: bad index, not setting anything.")
1.21 + return
1.22 + end if
1.23 + lines(nline:nline+nnewlines-1) = newlines
1.24 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.25 + nline = nline + nnewlines
1.26 + return
1.27 +end
1.28 +;**************************************************************
1.29 +; Main code.
1.30 +begin
1.31 +
1.32 + plot_type = "ps"
1.33 + plot_type_new = "png"
1.34 +
1.35 +;************************************************
1.36 +; read data: model
1.37 +;************************************************
1.38 + co2_i = 283.1878
1.39 + co2_f = 364.1252
1.40 +
1.41 + model_grid = "T42"
1.42 +
1.43 + model_name_i = "i01.07cn"
1.44 + model_name_f = "i01.10cn"
1.45 +
1.46 +;model_name_i = "i01.07casa"
1.47 +;model_name_f = "i01.10casa"
1.48 +
1.49 + model_name = model_name_f
1.50 +
1.51 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.52 + film_i = model_name_i + "_1990-2004_ANN_climo.nc"
1.53 + film_f = model_name_f + "_1990-2004_ANN_climo.nc"
1.54 +
1.55 + fm_i = addfile (dirm+film_i,"r")
1.56 + fm_f = addfile (dirm+film_f,"r")
1.57 +
1.58 + xm = fm_f->lon
1.59 + ym = fm_f->lat
1.60 +
1.61 + npp_i = fm_i->NPP
1.62 + npp_f = fm_f->NPP
1.63 +
1.64 +;Units for these variables are:
1.65 +;npp_i: g C/m^2/s
1.66 +
1.67 + nsec_per_year = 60*60*24*365
1.68 +
1.69 + npp_i = npp_i * nsec_per_year
1.70 + npp_f = npp_f * nsec_per_year
1.71 +
1.72 +;===================================================
1.73 +; read data: observed -station
1.74 +;===================================================
1.75 +
1.76 + station = (/"DukeFACE" \
1.77 + ,"AspenFACE" \
1.78 + ,"ORNL-FACE" \
1.79 + ,"POP-EUROFACE" \
1.80 + /)
1.81 +
1.82 + lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
1.83 + lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
1.84 + lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
1.85 +;print (lon_ob)
1.86 +
1.87 + n_sta = dimsizes(station)
1.88 + beta_4_ob = new((/n_sta/),float)
1.89 + beta_4_ob = 0.60
1.90 +
1.91 +;===================================================
1.92 +; get model data at station
1.93 +;===================================================
1.94 +
1.95 + npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
1.96 +
1.97 + npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
1.98 +
1.99 +;print (npp_i_4)
1.100 +;print (npp_f_4)
1.101 +;============================
1.102 +;compute beta_4
1.103 +;============================
1.104 + beta_4 = new((/n_sta/),float)
1.105 +
1.106 + beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
1.107 +
1.108 + beta_4_avg = avg(beta_4)
1.109 +
1.110 +;print (beta_4)
1.111 +;print (beta_4_avg)
1.112 +
1.113 +;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
1.114 +
1.115 + bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
1.116 + M_beta = (1. - (bias/n_sta))*3.
1.117 +
1.118 + print (M_beta)
1.119 +
1.120 +;=========================
1.121 +; for html table - station
1.122 +;=========================
1.123 +
1.124 + output_html = "table_station.html"
1.125 +
1.126 +; column (not including header column)
1.127 +
1.128 + col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
1.129 +
1.130 + ncol = dimsizes(col_head)
1.131 +
1.132 +; row (not including header row)
1.133 + row_head = (/"DukeFACE" \
1.134 + ,"AspenFACE" \
1.135 + ,"ORNL-FACE" \
1.136 + ,"POP-EUROFACE" \
1.137 + ,"All Station" \
1.138 + /)
1.139 + nrow = dimsizes(row_head)
1.140 +
1.141 +; arrays to be passed to table.
1.142 + text4 = new ((/nrow, ncol/),string )
1.143 +
1.144 + do i=0,nrow-2
1.145 + text4(i,0) = sprintf("%.1f",lat_ob(i))
1.146 + text4(i,1) = sprintf("%.1f",lon_ob(i))
1.147 + text4(i,2) = sprintf("%.1f",co2_i)
1.148 + text4(i,3) = sprintf("%.1f",co2_f)
1.149 + text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
1.150 + text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
1.151 + text4(i,6) = sprintf("%.1f",beta_4(i))
1.152 + text4(i,7) = "-"
1.153 + end do
1.154 + text4(nrow-1,0) = "-"
1.155 + text4(nrow-1,1) = "-"
1.156 + text4(nrow-1,2) = "-"
1.157 + text4(nrow-1,3) = "-"
1.158 + text4(nrow-1,4) = "-"
1.159 + text4(nrow-1,5) = "-"
1.160 + text4(nrow-1,6) = sprintf("%.1f",beta_4_avg)
1.161 + text4(nrow-1,7) = sprintf("%.1f",avg(beta_4_ob))
1.162 +
1.163 +;-----------
1.164 +; html table
1.165 +;-----------
1.166 +
1.167 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.168 +
1.169 + header = (/"<HTML>" \
1.170 + ,"<HEAD>" \
1.171 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.172 + ,"</HEAD>" \
1.173 + ,header_text \
1.174 + /)
1.175 + footer = "</HTML>"
1.176 +
1.177 + table_header = (/ \
1.178 + "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
1.179 + ,"<tr>" \
1.180 + ," <th bgcolor=DDDDDD >Station</th>" \
1.181 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.182 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.183 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.184 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.185 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.186 + ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
1.187 + ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
1.188 + ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
1.189 + ,"</tr>" \
1.190 + /)
1.191 + table_footer = "</table>"
1.192 + row_header = "<tr>"
1.193 + row_footer = "</tr>"
1.194 +
1.195 + lines = new(50000,string)
1.196 + nline = 0
1.197 +
1.198 + set_line(lines,nline,header)
1.199 + set_line(lines,nline,table_header)
1.200 +;-----------------------------------------------
1.201 +;row of table
1.202 +
1.203 + do n = 0,nrow-1
1.204 + set_line(lines,nline,row_header)
1.205 +
1.206 + txt1 = row_head(n)
1.207 + txt2 = text4(n,0)
1.208 + txt3 = text4(n,1)
1.209 + txt4 = text4(n,2)
1.210 + txt5 = text4(n,3)
1.211 + txt6 = text4(n,4)
1.212 + txt7 = text4(n,5)
1.213 + txt8 = text4(n,6)
1.214 + txt9 = text4(n,7)
1.215 +
1.216 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.217 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.218 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.219 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.220 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.221 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.222 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.223 + set_line(lines,nline,"<th>"+txt8+"</th>")
1.224 + set_line(lines,nline,"<th>"+txt9+"</th>")
1.225 +
1.226 + set_line(lines,nline,row_footer)
1.227 + end do
1.228 +;-----------------------------------------------
1.229 + set_line(lines,nline,table_footer)
1.230 + set_line(lines,nline,footer)
1.231 +
1.232 +; Now write to an HTML file.
1.233 + idx = ind(.not.ismissing(lines))
1.234 + if(.not.any(ismissing(idx))) then
1.235 + asciiwrite(output_html,lines(idx))
1.236 + else
1.237 + print ("error?")
1.238 + end if
1.239 +
1.240 + delete (col_head)
1.241 + delete (row_head)
1.242 + delete (text4)
1.243 + delete (table_header)
1.244 + delete (idx)
1.245 +
1.246 +;************************************************
1.247 +; read data: observed-2
1.248 +;************************************************
1.249 +
1.250 + ob_name = "MODIS MOD 15A2 2000-2005"
1.251 +
1.252 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.253 + filo = "land_class_"+model_grid+".nc"
1.254 +
1.255 + fo = addfile(diro+filo,"r")
1.256 +
1.257 + classob = tofloat(fo->LAND_CLASS)
1.258 +
1.259 + nclass = 20
1.260 +
1.261 +;*******************************************************************
1.262 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.263 +;********************************************************************
1.264 +
1.265 + nclassn = nclass + 1
1.266 + range = fspan(0,nclassn-1,nclassn)
1.267 +; print (range)
1.268 +
1.269 +; Use this range information to grab all the values in a
1.270 +; particular range, and then take an average.
1.271 +
1.272 + nr = dimsizes(range)
1.273 + nx = nr-1
1.274 + xvalues = new((/2,nx/),float)
1.275 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.276 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.277 + dx4 = dx/4 ; 1/4 of the range
1.278 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.279 +; get data
1.280 +
1.281 + base_1D = ndtooned(classob)
1.282 + data1_1D = ndtooned(npp_i)
1.283 + data2_1D = ndtooned(npp_f)
1.284 +
1.285 +; output
1.286 +
1.287 + yvalues = new((/2,nx/),float)
1.288 +
1.289 + do nd=0,1
1.290 +
1.291 +; See if we are doing data1 (nd=0) or data2 (nd=1).
1.292 +
1.293 + base = base_1D
1.294 +
1.295 + if(nd.eq.0) then
1.296 + data = data1_1D
1.297 + else
1.298 + data = data2_1D
1.299 + end if
1.300 +
1.301 +; Loop through each range, using base.
1.302 +
1.303 + do i=0,nr-2
1.304 + if (i.ne.(nr-2)) then
1.305 +; print("")
1.306 +; print("In range ["+range(i)+","+range(i+1)+")")
1.307 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.308 + else
1.309 +; print("")
1.310 +; print("In range ["+range(i)+",)")
1.311 + idx = ind(base.ge.range(i))
1.312 + end if
1.313 +
1.314 +; Calculate average
1.315 +
1.316 + if(.not.any(ismissing(idx))) then
1.317 + yvalues(nd,i) = avg(data(idx))
1.318 + count = dimsizes(idx)
1.319 + else
1.320 + yvalues(nd,i) = yvalues@_FillValue
1.321 + count = 0
1.322 + end if
1.323 +
1.324 +;#############################################################
1.325 +; set the following 4 classes to _FillValue:
1.326 +; Water Bodies(0), Urban and Build-Up(13),
1.327 +; Permenant Snow and Ice(15), Unclassified(17)
1.328 +
1.329 + if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
1.330 + yvalues(nd,i) = yvalues@_FillValue
1.331 + count = 0
1.332 + end if
1.333 +;#############################################################
1.334 +
1.335 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.336 +
1.337 +; Clean up for next time in loop.
1.338 +
1.339 + delete(idx)
1.340 + end do
1.341 +
1.342 + delete(data)
1.343 + end do
1.344 +
1.345 +
1.346 +;============================
1.347 +;compute beta
1.348 +;============================
1.349 +
1.350 + u = yvalues(0,:)
1.351 + v = yvalues(1,:)
1.352 +
1.353 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.354 +
1.355 + uu = u(good)
1.356 + vv = v(good)
1.357 +
1.358 + n_biome = dimsizes(uu)
1.359 +
1.360 + beta_biome = new((/n_biome/),float)
1.361 +
1.362 + beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
1.363 +
1.364 + beta_biome_avg = avg(beta_biome)
1.365 +
1.366 +;print (beta_biome_avg)
1.367 +
1.368 +;===========================
1.369 +; for html table - biome
1.370 +;===========================
1.371 +
1.372 + output_html = "table_biome.html"
1.373 +
1.374 +; column (not including header column)
1.375 +
1.376 + col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
1.377 +
1.378 + ncol = dimsizes(col_head)
1.379 +
1.380 +; row (not including header row)
1.381 +; 3 classes removed: Water Bodies, Urban and Build-Up, Unclassified
1.382 +; function "good" removed the last 2 classes
1.383 +; text4(i,2) = sprintf("%.2f",uu(i+1)) remove the first class
1.384 +
1.385 + row_head = (/"Evergreen Needleleaf Forests" \
1.386 + ,"Evergreen Broadleaf Forests" \
1.387 + ,"Deciduous Needleleaf Forest" \
1.388 + ,"Deciduous Broadleaf Forests" \
1.389 + ,"Mixed Forests" \
1.390 + ,"Closed Bushlands" \
1.391 + ,"Open Bushlands" \
1.392 + ,"Woody Savannas (S. Hem.)" \
1.393 + ,"Savannas (S. Hem.)" \
1.394 + ,"Grasslands" \
1.395 + ,"Permanent Wetlands" \
1.396 + ,"Croplands" \
1.397 + ,"Cropland/Natural Vegetation Mosaic" \
1.398 + ,"Barren or Sparsely Vegetated" \
1.399 + ,"Woody Savannas (N. Hem.)" \
1.400 + ,"Savannas (N. Hem.)" \
1.401 + ,"All Biome" \
1.402 + /)
1.403 + nrow = dimsizes(row_head)
1.404 +
1.405 +; arrays to be passed to table.
1.406 + text4 = new ((/nrow, ncol/),string )
1.407 +
1.408 + do i=0,nrow-2
1.409 + text4(i,0) = sprintf("%.1f",co2_i)
1.410 + text4(i,1) = sprintf("%.1f",co2_f)
1.411 + text4(i,2) = sprintf("%.1f",uu(i))
1.412 + text4(i,3) = sprintf("%.1f",vv(i))
1.413 + text4(i,4) = sprintf("%.1f",beta_biome(i))
1.414 + end do
1.415 + text4(nrow-1,0) = "-"
1.416 + text4(nrow-1,1) = "-"
1.417 + text4(nrow-1,2) = "-"
1.418 + text4(nrow-1,3) = "-"
1.419 + text4(nrow-1,4) = sprintf("%.1f",beta_biome_avg)
1.420 +
1.421 +;**************************************************
1.422 +; html table
1.423 +;**************************************************
1.424 +
1.425 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.426 +
1.427 + header = (/"<HTML>" \
1.428 + ,"<HEAD>" \
1.429 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.430 + ,"</HEAD>" \
1.431 + ,header_text \
1.432 + /)
1.433 + footer = "</HTML>"
1.434 +
1.435 + table_header = (/ \
1.436 + "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
1.437 + ,"<tr>" \
1.438 + ," <th bgcolor=DDDDDD >Biome Class</th>" \
1.439 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.440 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.441 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.442 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.443 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.444 + ,"</tr>" \
1.445 + /)
1.446 + table_footer = "</table>"
1.447 + row_header = "<tr>"
1.448 + row_footer = "</tr>"
1.449 +
1.450 + lines = new(50000,string)
1.451 + nline = 0
1.452 +
1.453 + set_line(lines,nline,header)
1.454 + set_line(lines,nline,table_header)
1.455 +;-----------------------------------------------
1.456 +;row of table
1.457 +
1.458 + do n = 0,nrow-1
1.459 + set_line(lines,nline,row_header)
1.460 +
1.461 + txt1 = row_head(n)
1.462 + txt2 = text4(n,0)
1.463 + txt3 = text4(n,1)
1.464 + txt4 = text4(n,2)
1.465 + txt5 = text4(n,3)
1.466 + txt6 = text4(n,4)
1.467 +
1.468 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.469 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.470 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.471 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.472 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.473 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.474 +
1.475 + set_line(lines,nline,row_footer)
1.476 + end do
1.477 +;-----------------------------------------------
1.478 + set_line(lines,nline,table_footer)
1.479 + set_line(lines,nline,footer)
1.480 +
1.481 +; Now write to an HTML file.
1.482 + idx = ind(.not.ismissing(lines))
1.483 + if(.not.any(ismissing(idx))) then
1.484 + asciiwrite(output_html,lines(idx))
1.485 + else
1.486 + print ("error?")
1.487 + end if
1.488 +
1.489 +end
1.490 +