clamp: npp/99.all.ncl.2@4be95183fbcd

     1 ;*****************************************************

     2 ; combine scatter, histogram, global and zonal plots

     3 ; compute all correlation coef and M score

     4 ; add 1-to-1 line in scatter plots

     5 ;*****************************************************

     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"

     7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"

     8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     9 ;*****************************************************

    10 procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \

    11                  ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \

    12                  ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \

    13                  ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \

    14                  ,dx4[1]:numeric \

    15                   )

    16 begin

    17 ; Calculaee "nice" bins for binning the data in equally spaced ranges.

    18 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D

    19 ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4

    21   nbins       = 15     ; Number of bins to use.

    23   nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)

    24   nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)

    25   range       = fspan(nicevals(0),nicevals(1),nvals)

    27 ; Use this range information to grab all the values in a

    28 ; particular range, and then take an average.

    30   nr           = dimsizes(range)

    31   nx           = nr-1

    32 ; print (nx)

    34 ; xvalues      = new((/2,nx/),typeof(RAIN1_1D))

    35   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

    36   dx           = xvalues(0,1) - xvalues(0,0)       ; range width

    37   dx4          = dx/4                              ; 1/4 of the range

    38   xvalues(1,:) = xvalues(0,:) - dx/5.

    39 ; yvalues      = new((/2,nx/),typeof(RAIN1_1D))

    40 ; mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

    41 ; mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

    43   do nd=0,1

    45 ; See if we are doing model or observational data.

    47     if(nd.eq.0) then

    48       data     = RAIN1_1D

    49       npp_data = NPP1_1D

    50     else

    51       data     = RAIN2_1D

    52       npp_data = NPP2_1D

    53     end if

    55 ; Loop through each range and check for values.

    57     do i=0,nr-2

    58       if (i.ne.(nr-2)) then

    59 ;        print("")

    60 ;        print("In range ["+range(i)+","+range(i+1)+")")

    61         idx = ind((range(i).le.data).and.(data.lt.range(i+1)))

    62       else

    63 ;        print("")

    64 ;        print("In range ["+range(i)+",)")

    65         idx = ind(range(i).le.data)

    66       end if

    68 ; Calculate average, and get min and max.

    70       if(.not.any(ismissing(idx))) then

    71         yvalues(nd,i)    = avg(npp_data(idx))

    72         mn_yvalues(nd,i) = min(npp_data(idx))

    73         mx_yvalues(nd,i) = max(npp_data(idx))

    74         count = dimsizes(idx)

    75       else

    76         count            = 0

    77         yvalues(nd,i)    = yvalues@_FillValue

    78         mn_yvalues(nd,i) = yvalues@_FillValue

    79         mx_yvalues(nd,i) = yvalues@_FillValue

    80       end if

    82 ; Print out information.

    83 ;      print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))

    84 ;      print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

    86 ; Clean up for next time in loop.

    88       delete(idx)

    89     end do

    90     delete(data)

    91     delete(npp_data)

    92   end do

    93 end

    94 ;****************************************************************************

    95 ; Main code.

    96 ;****************************************************************************

    97 begin

    99  plot_type     = "ps"

   100  plot_type_new = "png"

   102 ;************************************************

   103 ; read model data

   104 ;************************************************

   106 ;film = "i01.06cn_1798-2004_ANN_climo.nc"

   107 ;model_name = "06cn"

   108 ;model_grid = "T42"

   110 ;film = "i01.06casa_1798-2004_ANN_climo.nc"

   111 ;model_name = "06casa"

   112 ;model_grid = "T42"

   114 ;film = "i01.10casa_1948-2004_ANN_climo.nc"

   115 ;model_name = "10casa"

   116 ;model_grid = "T42"

   118  film = "i01.10cn_1948-2004_ANN_climo.nc"

   119  model_name = "10cn"

   120  model_grid = "T42"

   122  html_name = "table.html." + model_name

   123  html_new  = html_name +".new"

   124  system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new)

   125  system("mv -f "+html_new+" "+html_name)

   126 ;--------------------------------------------------

   127  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

   128  fm   = addfile (dirm+film,"r")

   130  nppmod0  = fm->NPP

   131  rainmod0 = fm->RAIN

   132  xm       = fm->lon

   133  ym       = fm->lat

   135 ;************************************************

   136 ; read ob data

   137 ;************************************************

   139 ;(1) data at 81 sites

   140  dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

   141  fil81 = "data.81.nc"

   142  f81   = addfile (dir81+fil81,"r")

   144  id81   = f81->SITE_ID

   145  npp81  = f81->TNPP_C

   146  rain81 = tofloat(f81->PREC_ANN)

   147  x81    = f81->LONG_DD

   148  y81    = f81->LAT_DD

   150  id81@long_name  = "SITE_ID"

   152 ;change longitude from (-180,180) to (0,360)

   153 ;for model data interpolation

   155  do i= 0,dimsizes(x81)-1

   156     if (x81(i) .lt. 0.) then

   157         x81(i) = x81(i)+ 360.

   158     end if

   159  end do

   160 ;print (x81)

   161 ;-------------------------------------

   162 ;(2) data at 933 sites

   163  dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

   164  fil933 = "data.933.nc"

   165  f933   = addfile (dir933+fil933,"r")

   167  id933   = f933->SITE_ID

   168  npp933  = f933->TNPP_C

   169  rain933 = f933->PREC

   170  x933    = f933->LONG_DD

   171  y933    = f933->LAT_DD

   173  id933@long_name  = "SITE_ID"

   175 ;change longitude from (-180,180) to (0,360)

   176 ;for model data interpolation

   178  do i= 0,dimsizes(x933)-1

   179     if (x933(i) .lt. 0.) then

   180         x933(i) = x933(i)+ 360.

   181     end if

   182  end do

   183 ;print (x933)

   184 ;----------------------------------------

   185 ;(3) global data, interpolated into model grid

   186  dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

   187  filglobe = "Npp_"+model_grid+"_mean.nc"

   188  fglobe   = addfile (dirglobe+filglobe,"r")

   190  nppglobe0 = fglobe->NPP

   191  nppglobe  = nppglobe0

   193 ;***********************************************************************

   194 ; interpolate model data into ob sites

   195 ;***********************************************************************

   197  nppmod   = nppmod0(0,:,:)

   198  rainmod  = rainmod0(0,:,:)

   199  delete (nppmod0)

   200  delete (rainmod0)

   202  nppmod81  =linint2_points(xm,ym,nppmod,True,x81,y81,0)

   204  nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)

   206  rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)

   208  rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)

   210 ;**********************************************************

   211 ; unified units

   212 ;**********************************************************

   213 ; Units for these variables are:

   214 ;

   215 ; rain81  : mm/year

   216 ; rainmod : mm/s

   217 ; npp81   : g C/m^2/year

   218 ; nppmod81: g C/m^2/s

   219 ; nppglobe: g C/m^2/year

   220 ;

   221 ; We want to convert these to "m/year" and "g C/m^2/year".

   223   nsec_per_year = 60*60*24*365

   225   rain81    = rain81 / 1000.

   226   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   227   nppmod81  = nppmod81 * nsec_per_year

   229   rain933    = rain933 / 1000.

   230   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   231   nppmod933  = nppmod933 * nsec_per_year

   233   nppmod  = nppmod * nsec_per_year

   235   npp81@units      = "gC/m^2/yr"

   236   nppmod81@units   = "gC/m^2/yr"

   237   npp933@units     = "gC/m^2/yr"

   238   nppmod933@units  = "gC/m^2/yr"

   239   nppmod@units     = "gC/m^2/yr"

   240   nppglobe@units   = "gC/m^2/yr"

   241   rain81@units     = "m/yr"

   242   rainmod81@units  = "m/yr"

   243   rain933@units    = "m/yr"

   244   rainmod933@units = "m/yr"

   246   npp81@long_name      = "NPP (gC/m2/year)"

   247   npp933@long_name     = "NPP (gC/m2/year)"

   248   nppmod81@long_name   = "NPP (gC/m2/year)"

   249   nppmod933@long_name  = "NPP (gC/m2/year)"

   250   nppmod@long_name     = "NPP (gC/m2/year)"

   251   nppglobe@long_name   = "NPP (gC/m2/year)"

   252   rain81@long_name     = "PREC (m/year)"

   253   rain933@long_name    = "PREC (m/year)"

   254   rainmod81@long_name  = "PREC (m/year)"

   255   rainmod933@long_name = "PREC (m/year)"

   257 ;*******************************************************************

   258 ;(A)-1 table of site81 (1) (the first 41 sites)

   259 ;*******************************************************************

   261   table_length = 0.95

   263 ; table header name

   264   table_header_name = "Site ID"

   266 ; column (not including header column)

   267   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)

   268   ncol = dimsizes(col_header)

   270 ; row (not including header row)

   271   nrow       = 41

   272   row_header = new ((/nrow/),string )

   273   row_header(0:nrow-1) = id81(0:nrow-1)

   275 ; Table header

   276   ncr1  = (/1,1/)               ; 1 row, 1 column

   277   x1    = (/0.005,0.15/)        ; Start and end X

   278   y1    = (/0.900,0.95/)       ; Start and end Y

   279   text1 = table_header_name

   280   res1               = True

   281   res1@txFontHeightF = 0.015

   282   res1@gsFillColor   = "CornFlowerBlue"

   284 ; Column header (equally space in x2)

   285   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   286   x2    = (/x1(1),0.995/)       ; start from end of x1

   287   y2    = y1                    ; same as y1

   288   text2 = col_header

   289   res2               = True

   290   res2@txFontHeightF = 0.015

   291   res2@gsFillColor   = "Gray"

   293 ; Row header (equally space in y2)

   294   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   295   x3    = x1                         ; same as x1

   296   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   297   text3 = row_header

   298   res3               = True

   299   res3@txFontHeightF = 0.010

   300   res3@gsFillColor   = "Gray"

   302 ; Main table body

   303   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   304   x4    = x2                      ; Start and end x

   305   y4    = y3                      ; Start and end Y

   306   text4 = new((/nrow,ncol/),string)

   308   color_fill4           = new((/nrow,ncol/),string)

   309   color_fill4           = "white"

   310 ; color_fill4(:,ncol-1) = "grey"

   311 ; color_fill4(nrow-1,:) = "green"

   313   res4               = True       ; Set up resource list

   314 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   315   res4@txFontHeightF = 0.015

   316   res4@gsFillColor   = color_fill4

   318   delete (color_fill4)

   319 ;-------------------------------------------------------------------

   320 ; for table value

   322   do n = 0,nrow-1

   323      text4(n,0) = sprintf("%5.2f", y81(n))

   324      text4(n,1) = sprintf("%5.2f", x81(n))

   325      text4(n,2) = sprintf("%5.2f", npp81(n))

   326      text4(n,3) = sprintf("%5.2f", rain81(n))

   327   end do

   328 ;-------------------------------------------------------------------

   330   plot_name = "table_site81_ob1"

   332   wks = gsn_open_wks (plot_type,plot_name)

   333 ;------------------------------------------

   334 ; for table title

   336   gRes               = True

   337   gRes@txFontHeightF = 0.02

   338 ; gRes@txAngleF      = 90

   340   title_text = "Observation at 81 Sites (1)"

   342   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   343 ;------------------------------------------

   345   gsn_table(wks,ncr1,x1,y1,text1,res1)

   346   gsn_table(wks,ncr2,x2,y2,text2,res2)

   347   gsn_table(wks,ncr3,x3,y3,text3,res3)

   348   gsn_table(wks,ncr4,x4,y4,text4,res4)

   350   frame(wks)

   351   clear (wks)

   353   delete (row_header)

   354   delete (res3)

   355   delete (ncr3)

   356   delete (text3)

   357   delete (res4)

   358   delete (ncr4)

   359   delete (text4)

   361   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   362   system("rm "+plot_name+"."+plot_type)

   363   system("rm "+plot_name+"-1."+plot_type_new)

   364   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   366 ;*******************************************************************

   367 ;(A)-2 table of site81 (2) (the last 40 sites)

   368 ;*******************************************************************

   370   table_length = 0.95

   372 ; table header name

   373   table_header_name = "Site ID"

   375 ; column (not including header column)

   376   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)

   377   ncol = dimsizes(col_header)

   379 ; row (not including header row)

   380   nrow       = 40

   381   row_header = new ((/nrow/),string )

   382   row_header(0:nrow-1) = id81(41:41+nrow-1)

   384 ; Table header

   385   ncr1  = (/1,1/)               ; 1 row, 1 column

   386   x1    = (/0.005,0.15/)        ; Start and end X

   387   y1    = (/0.900,0.95/)       ; Start and end Y

   388   text1 = table_header_name

   389   res1               = True

   390   res1@txFontHeightF = 0.015

   391   res1@gsFillColor   = "CornFlowerBlue"

   393 ; Column header (equally space in x2)

   394   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   395   x2    = (/x1(1),0.995/)       ; start from end of x1

   396   y2    = y1                    ; same as y1

   397   text2 = col_header

   398   res2               = True

   399   res2@txFontHeightF = 0.015

   400   res2@gsFillColor   = "Gray"

   402 ; Row header (equally space in y2)

   403   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   404   x3    = x1                         ; same as x1

   405   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   406   text3 = row_header

   407   res3               = True

   408   res3@txFontHeightF = 0.010

   409   res3@gsFillColor   = "Gray"

   411 ; Main table body

   412   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   413   x4    = x2                      ; Start and end x

   414   y4    = y3                      ; Start and end Y

   415   text4 = new((/nrow,ncol/),string)

   417   color_fill4           = new((/nrow,ncol/),string)

   418   color_fill4           = "white"

   419 ; color_fill4(:,ncol-1) = "grey"

   420 ; color_fill4(nrow-1,:) = "green"

   422   res4               = True       ; Set up resource list

   423 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   424   res4@txFontHeightF = 0.015

   425   res4@gsFillColor   = color_fill4

   427   delete (color_fill4)

   428 ;-------------------------------------------------------------------

   429 ; for table value

   431   do n = 0,nrow-1

   432      text4(n,0) = sprintf("%5.2f", y81(n+41))

   433      text4(n,1) = sprintf("%5.2f", x81(n+41))

   434      text4(n,2) = sprintf("%5.2f", npp81(n+41))

   435      text4(n,3) = sprintf("%5.2f", rain81(n+41))

   436   end do

   437 ;-------------------------------------------------------------------

   439   plot_name = "table_site81_ob2"

   441   wks = gsn_open_wks (plot_type,plot_name)

   442 ;------------------------------------------

   443 ; for table title

   445   gRes               = True

   446   gRes@txFontHeightF = 0.02

   447 ; gRes@txAngleF      = 90

   449   title_text = "Observation at 81 Sites (2)"

   451   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   452 ;------------------------------------------

   454   gsn_table(wks,ncr1,x1,y1,text1,res1)

   455   gsn_table(wks,ncr2,x2,y2,text2,res2)

   456   gsn_table(wks,ncr3,x3,y3,text3,res3)

   457   gsn_table(wks,ncr4,x4,y4,text4,res4)

   459   frame(wks)

   460   clear (wks)

   462   delete (row_header)

   463   delete (res3)

   464   delete (ncr3)

   465   delete (text3)

   466   delete (res4)

   467   delete (ncr4)

   468   delete (text4)

   470   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   471   system("rm "+plot_name+"."+plot_type)

   472   system("rm "+plot_name+"-1."+plot_type_new)

   473   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   475 ;**************************************************************************

   476 ;(A)-3 scatter plot, ob 933

   477 ;**************************************************************************

   479  plot_name = "scatter_ob_933"

   480  title     = "Observed 933 sites"

   482  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

   483  res                   = True                  ; plot mods desired

   484  res@tiMainString      = title                 ; add title

   485  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   486  res@xyMarkers         =  16                   ; choose type of marker

   487  res@xyMarkerColor     = "red"                 ; Marker color

   488  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   489  res@tmLabelAutoStride = True                  ; nice tick mark labels

   491  plot  = gsn_csm_xy (wks,id933,npp933,res)     ; create plot

   492  frame(wks)

   493  clear (wks)

   495  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   496  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   497  system("rm "+plot_name+"-*."+plot_type_new)

   498  system("rm "+plot_name+"."+plot_type)

   500 ;**************************************************************************

   501 ;(A)-4 scatter plot, model 933

   502 ;**************************************************************************

   504  plot_name = "scatter_model_933"

   505  title     = "Model "+ model_name +" 933 sites"

   507  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

   508  res                   = True                  ; plot mods desired

   509  res@tiMainString      = title                 ; add title

   510  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   511  res@xyMarkers         =  16                   ; choose type of marker

   512  res@xyMarkerColor     = "red"                 ; Marker color

   513  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   514  res@tmLabelAutoStride = True                  ; nice tick mark labels

   516  plot  = gsn_csm_xy (wks,id933,nppmod933,res)    ; create plot

   517  frame(wks)

   518  clear (wks)

   520  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   521  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   522  system("rm "+plot_name+"-*."+plot_type_new)

   523  system("rm "+plot_name+"."+plot_type)

   525 ;***************************************************************************

   526 ;(A)-5 scatter plot, model vs ob 81

   527 ;***************************************************************************

   529  plot_name = "scatter_model_vs_ob_81"

   530  title     = model_name +" vs ob 81 sites"

   532  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

   533  res                   = True                  ; plot mods desired

   534  res@tiMainString      = title                 ; add title

   535  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   536  res@xyMarkers         =  16                   ; choose type of marker

   537  res@xyMarkerColor     = "red"                 ; Marker color

   538  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   539  res@tmLabelAutoStride = True                  ; nice tick mark labels

   541  res@gsnDraw           = False

   542  res@gsnFrame          = False                 ; don't advance frame yet

   543 ;-------------------------------

   544 ;compute correlation coef. and M

   545  ccr81 = esccr(nppmod81,npp81,0)

   546 ;print (ccr81)

   548 ;new eq

   549  bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))

   550  M81  = (1. - (bias/dimsizes(y81)))*5.

   552  M_npp_S81  = sprintf("%.2f", M81)

   553  system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new)

   554  system("mv -f "+html_new+" "+html_name)

   555  print (M_npp_S81)

   557  tRes  = True

   558  tRes@txFontHeightF = 0.025

   560  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"

   562  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

   563 ;--------------------------------

   564  plot  = gsn_csm_xy (wks,npp81,nppmod81,res)       ; create plot

   565 ;-------------------------------

   566 ; add polyline

   568  dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)

   569 ;-------------------------------

   570  draw (plot)

   571  frame(wks)

   572  clear (wks)

   573  delete (dum)

   575  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   576  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   577  system("rm "+plot_name+"-*."+plot_type_new)

   578  system("rm "+plot_name+"."+plot_type)

   580 ;***************************************************************************

   581 ;(A)-6 scatter plot, model vs ob 933

   582 ;***************************************************************************

   584  plot_name = "scatter_model_vs_ob_933"

   585  title     = model_name +" vs ob 933 sites"

   587  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

   588  res                   = True                  ; plot mods desired

   589  res@tiMainString      = title                 ; add title

   590  res@xyMarkLineModes   = "Markers"             ; choose which have markers

   591  res@xyMarkers         =  16                   ; choose type of marker

   592  res@xyMarkerColor     = "red"                 ; Marker color

   593  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

   594  res@tmLabelAutoStride = True                  ; nice tick mark labels

   596  res@gsnDraw           = False

   597  res@gsnFrame          = False                 ; don't advance frame yet

   598 ;-------------------------------

   599 ;compute correlation coef. and M

   600  ccr933 = esccr(nppmod933,npp933,0)

   601 ;print (ccr933)

   603 ;new eq

   604  bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))

   605  M933 = (1. - (bias/dimsizes(y933)))*5.

   607  M_npp_S933  = sprintf("%.2f", M933)

   608  system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new)

   609  system("mv -f "+html_new+" "+html_name)

   610  print (M_npp_S933)

   612  tRes  = True

   613  tRes@txFontHeightF = 0.025

   615  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"

   617  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

   618 ;--------------------------------

   619  plot  = gsn_csm_xy (wks,npp933,nppmod933,res)    ; create plot

   620 ;-------------------------------

   621 ; add polyline

   623  dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)

   624 ;-------------------------------

   625  draw (plot)

   626  frame(wks)

   627  clear (wks)

   628  delete (dum)

   630  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   631  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   632  system("rm "+plot_name+"-*."+plot_type_new)

   633  system("rm "+plot_name+"."+plot_type)

   635 ;*******************************************************************

   636 ;(A)-7 table of site81 (1) (the first 41 sites), model vs ob

   637 ;*******************************************************************

   639   table_length = 0.95

   641 ; table header name

   642   table_header_name = "Site ID"

   644 ; row (not including header row)

   645   nrow       = 41

   646   row_header = new ((/nrow/),string )

   647   row_header(0:nrow-1) = id81(0:nrow-1)

   649 ; Table header

   650   ncr1  = (/1,1/)               ; 1 row, 1 column

   651   x1    = (/0.005,0.15/)        ; Start and end X

   652   y1    = (/0.85 ,0.95/)        ; Start and end Y

   653   text1 = table_header_name

   654   res1               = True

   655   res1@txFontHeightF = 0.015

   656   res1@gsFillColor   = "CornFlowerBlue"

   658 ; Column header1 (equally space in x2)

   660   delete (col_header)

   661   delete (ncr2)

   662   delete (text2)

   664   col_header = (/"Latitude","Longitude"/)

   665   ncol = dimsizes(col_header)

   667   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   668   x2    = (/x1(1),0.35/)        ; start from end of x1

   669   y2    = y1                    ; same as y1

   670   text2 = col_header

   671   res2               = True

   672   res2@txFontHeightF = 0.015

   673   res2@gsFillColor   = "Gray"

   675 ; Column header2 (equally space in x2)

   677   col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   678   col_header2 = (/"ob",model_name \

   679                  ,"ob",model_name \

   680                  /)

   682   ncol1 = dimsizes(col_header1)

   683   ncol2 = dimsizes(col_header2)

   685   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

   686   x21    = (/x2(1),0.995/)        ; start from end of x2

   687   yhalf  = (y1(0)+y1(1))*0.5

   688   y21    = (/yhalf,y1(1)/)        ; top half of y1

   689   text21 = col_header1

   690   res21               = True

   691   res21@txFontHeightF = 0.015

   692   res21@gsFillColor   = "Gray"

   694   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   695   x22    = x21                    ; start from end of x1

   696   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   697   text22 = col_header2

   698   res22               = True

   699   res22@txFontHeightF = 0.015

   700   res22@gsFillColor   = "Gray"

   702 ; Row header (equally space in y2)

   703   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   704   x3    = x1                         ; same as x1

   705   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   706   text3 = row_header

   707   res3               = True

   708   res3@txFontHeightF = 0.010

   709   res3@gsFillColor   = "Gray"

   711 ; Main table body-1

   712   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   713   x4    = x2                      ; Start and end x

   714   y4    = y3                      ; Start and end Y

   715   text4 = new((/nrow,ncol/),string)

   717   color_fill4           = new((/nrow,ncol/),string)

   718   color_fill4           = "white"

   719 ; color_fill4(:,ncol-1) = "grey"

   720 ; color_fill4(nrow-1,:) = "green"

   722   res4               = True       ; Set up resource list

   723 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   724   res4@txFontHeightF = 0.015

   725   res4@gsFillColor   = color_fill4

   727   delete (color_fill4)

   729 ; Main table body-2

   730   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   731   x42    = x21                      ; Start and end x

   732   y42    = y3                       ; Start and end Y

   733   text42 = new((/nrow,ncol2/),string)

   735   color_fill42           = new((/nrow,ncol2/),string)

   736   color_fill42           = "white"

   737 ; color_fill42(:,ncol-1) = "grey"

   738 ; color_fill42(nrow-1,:) = "green"

   740   res42               = True       ; Set up resource list

   741 ; res42@gsnDebug      = True       ; Useful to print NDC row,col values used.

   742   res42@txFontHeightF = 0.015

   743   res42@gsFillColor   = color_fill42

   745   delete (color_fill42)

   746 ;-------------------------------------------------------------------

   747 ; for table value

   749   do n = 0,nrow-1

   750      text4(n,0) = sprintf("%5.2f", y81(n))

   751      text4(n,1) = sprintf("%5.2f", x81(n))

   752   end do

   754   do n = 0,nrow-1

   755      text42(n,0) = sprintf("%5.2f", npp81(n))

   756      text42(n,1) = sprintf("%5.2f", nppmod81(n))

   757      text42(n,2) = sprintf("%5.2f", rain81(n))

   758      text42(n,3) = sprintf("%5.2f", rainmod81(n))

   759   end do

   760 ;---------------------------------------------------------------------------

   762   plot_name = "table_site81_model_vs_ob1"

   764   wks = gsn_open_wks (plot_type,plot_name)

   765 ;------------------------------------------

   766 ; for table title

   768   gRes               = True

   769   gRes@txFontHeightF = 0.02

   770 ; gRes@txAngleF      = 90

   772   title_text = "Model vs Observation at 81 Sites (1)"

   774   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   775 ;------------------------------------------

   777   gsn_table(wks,ncr1,x1,y1,text1,res1)

   778   gsn_table(wks,ncr2,x2,y2,text2,res2)

   779   gsn_table(wks,ncr21,x21,y21,text21,res21)

   780   gsn_table(wks,ncr22,x22,y22,text22,res22)

   781   gsn_table(wks,ncr3,x3,y3,text3,res3)

   782   gsn_table(wks,ncr4,x4,y4,text4,res4)

   783   gsn_table(wks,ncr42,x42,y42,text42,res42)

   785   frame(wks)

   786   clear (wks)

   788   delete (col_header)

   789   delete (row_header)

   790   delete (res1)

   791   delete (ncr1)

   792   delete (text1)

   793   delete (res2)

   794   delete (ncr2)

   795   delete (text2)

   796   delete (res21)

   797   delete (ncr21)

   798   delete (text21)

   799   delete (res22)

   800   delete (ncr22)

   801   delete (text22)

   802   delete (res3)

   803   delete (ncr3)

   804   delete (text3)

   805   delete (res4)

   806   delete (ncr4)

   807   delete (text4)

   808   delete (res42)

   809   delete (ncr42)

   810   delete (text42)

   812   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   813   system("rm "+plot_name+"."+plot_type)

   814   system("rm "+plot_name+"-1."+plot_type_new)

   815   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   817 ;*******************************************************************

   818 ;(A)-8 table of site81 (2) (the last 40 sites), model vs ob

   819 ;*******************************************************************

   821   table_length = 0.95

   823 ; table header name

   824   table_header_name = "Site ID"

   826 ; row (not including header row)

   827   nrow       = 40

   828   row_header = new ((/nrow/),string )

   829   row_header(0:nrow-1) = id81(41:41+nrow-1)

   831 ; Table header

   832   ncr1  = (/1,1/)               ; 1 row, 1 column

   833   x1    = (/0.005,0.15/)        ; Start and end X

   834   y1    = (/0.85 ,0.95/)        ; Start and end Y

   835   text1 = table_header_name

   836   res1               = True

   837   res1@txFontHeightF = 0.015

   838   res1@gsFillColor   = "CornFlowerBlue"

   840 ; Column header1 (equally space in x2)

   842   col_header = (/"Latitude","Longitude"/)

   843   ncol = dimsizes(col_header)

   845   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

   846   x2    = (/x1(1),0.35/)        ; start from end of x1

   847   y2    = y1                    ; same as y1

   848   text2 = col_header

   849   res2               = True

   850   res2@txFontHeightF = 0.015

   851   res2@gsFillColor   = "Gray"

   853 ; Column header2 (equally space in x2)

   855 ; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   856   col_header2 = (/"ob",model_name \

   857                  ,"ob",model_name \

   858                  /)

   860   ncol1 = dimsizes(col_header1)

   861   ncol2 = dimsizes(col_header2)

   863   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

   864   x21    = (/x2(1),0.995/)        ; start from end of x2

   865   yhalf  = (y1(0)+y1(1))*0.5

   866   y21    = (/yhalf,y1(1)/)        ; top half of y1

   867   text21 = col_header1

   868   res21               = True

   869   res21@txFontHeightF = 0.015

   870   res21@gsFillColor   = "Gray"

   872   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   873   x22    = x21                    ; start from end of x1

   874   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   875   text22 = col_header2

   876   res22               = True

   877   res22@txFontHeightF = 0.015

   878   res22@gsFillColor   = "Gray"

   880 ; Row header (equally space in y2)

   881   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   882   x3    = x1                         ; same as x1

   883   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   884   text3 = row_header

   885   res3               = True

   886   res3@txFontHeightF = 0.010

   887   res3@gsFillColor   = "Gray"

   889 ; Main table body-1

   890   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   891   x4    = x2                      ; Start and end x

   892   y4    = y3                      ; Start and end Y

   893   text4 = new((/nrow,ncol/),string)

   895   color_fill4           = new((/nrow,ncol/),string)

   896   color_fill4           = "white"

   897 ; color_fill4(:,ncol-1) = "grey"

   898 ; color_fill4(nrow-1,:) = "green"

   900   res4               = True       ; Set up resource list

   901 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   902   res4@txFontHeightF = 0.015

   903   res4@gsFillColor   = color_fill4

   905   delete (color_fill4)

   907 ; Main table body-2

   908   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   909   x42    = x21                      ; Start and end x

   910   y42    = y3                       ; Start and end Y

   911   text42 = new((/nrow,ncol2/),string)

   913   color_fill42           = new((/nrow,ncol2/),string)

   914   color_fill42           = "white"

   915 ; color_fill42(:,ncol-1) = "grey"

   916 ; color_fill42(nrow-1,:) = "green"

   918   res42               = True       ; Set up resource list

   919 ; res42@gsnDebug      = True       ; Useful to print NDC row,col values used.

   920   res42@txFontHeightF = 0.015

   921   res42@gsFillColor   = color_fill42

   923   delete (color_fill42)

   924 ;-------------------------------------------------------------------

   925 ; for table value

   927   do n = 0,nrow-1

   928      text4(n,0) = sprintf("%5.2f", y81(n+41))

   929      text4(n,1) = sprintf("%5.2f", x81(n+41))

   930   end do

   932   do n = 0,nrow-1

   933      text42(n,0) = sprintf("%5.2f", npp81(n+41))

   934      text42(n,1) = sprintf("%5.2f", nppmod81(n+41))

   935      text42(n,2) = sprintf("%5.2f", rain81(n+41))

   936      text42(n,3) = sprintf("%5.2f", rainmod81(n+41))

   937   end do

   938 ;-------------------------------------------------------------------

   940   plot_name = "table_site81_model_vs_ob2"

   942   wks = gsn_open_wks (plot_type,plot_name)

   943 ;------------------------------------------

   944 ; for table title

   946   gRes               = True

   947   gRes@txFontHeightF = 0.02

   948 ; gRes@txAngleF      = 90

   950   title_text = "Model vs Observation at 81 Sites (2)"

   952   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

   953 ;------------------------------------------

   955   gsn_table(wks,ncr1,x1,y1,text1,res1)

   956   gsn_table(wks,ncr2,x2,y2,text2,res2)

   957   gsn_table(wks,ncr21,x21,y21,text21,res21)

   958   gsn_table(wks,ncr22,x22,y22,text22,res22)

   959   gsn_table(wks,ncr3,x3,y3,text3,res3)

   960   gsn_table(wks,ncr4,x4,y4,text4,res4)

   961   gsn_table(wks,ncr42,x42,y42,text42,res42)

   963   frame(wks)

   964   clear (wks)

   966   delete (col_header)

   967   delete (row_header)

   968   delete (res1)

   969   delete (ncr1)

   970   delete (text1)

   971   delete (res2)

   972   delete (ncr2)

   973   delete (text2)

   974   delete (res21)

   975   delete (ncr21)

   976   delete (text21)

   977   delete (res22)

   978   delete (ncr22)

   979   delete (text22)

   980   delete (res3)

   981   delete (ncr3)

   982   delete (text3)

   983   delete (res4)

   984   delete (ncr4)

   985   delete (text4)

   986   delete (res42)

   987   delete (ncr42)

   988   delete (text42)

   990   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   991   system("rm "+plot_name+"."+plot_type)

   992   system("rm "+plot_name+"-1."+plot_type_new)

   993   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   995 ;**************************************************************************

   996 ;(B) histogram 81

   997 ;**************************************************************************

   998 ;get data

   999   RAIN1_1D = ndtooned(rain81)

  1000   RAIN2_1D = ndtooned(rainmod81)

  1001   NPP1_1D  = ndtooned(npp81)

  1002   NPP2_1D  = ndtooned(nppmod81)

  1004 ; number of bin

  1005   nx = 8

  1007   xvalues      = new((/2,nx/),float)

  1008   yvalues      = new((/2,nx/),float)

  1009   mn_yvalues   = new((/2,nx/),float)

  1010   mx_yvalues   = new((/2,nx/),float)

  1011   dx4          = new((/1/),float)

  1013   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

  1014          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

  1016 ;----------------------------------------

  1017 ;compute correlation coeff and M score

  1019  u = yvalues(0,:)

  1020  v = yvalues(1,:)

  1022  good = ind(.not.ismissing(u) .and. .not.ismissing(v))

  1023  uu = u(good)

  1024  vv = v(good)

  1026  ccr81h = esccr(uu,vv,0)

  1027 ;print (ccr81h)

  1029 ;new eq

  1030  bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))

  1031  M81h = (1.- (bias/dimsizes(uu)))*5.

  1033  M_npp_H81  = sprintf("%.2f", M81h)

  1034  system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new)

  1035  system("mv -f "+html_new+" "+html_name)

  1036  print (M_npp_H81)

  1038  delete (u)

  1039  delete (v)

  1040  delete (uu)

  1041  delete (vv)

  1042 ;----------------------------------------------------------------------

  1043 ; histogram res

  1044   resh                = True

  1045   resh@gsnMaximize    = True

  1046   resh@gsnDraw        = False

  1047   resh@gsnFrame       = False

  1048   resh@xyMarkLineMode = "Markers"

  1049   resh@xyMarkerSizeF  = 0.014

  1050   resh@xyMarker       = 16

  1051   resh@xyMarkerColors = (/"Brown","Blue"/)

  1052   resh@trYMinF        = min(mn_yvalues) - 10.

  1053   resh@trYMaxF        = max(mx_yvalues) + 10.

  1055   resh@tiYAxisString  = "NPP (g C/m2/year)"

  1056   resh@tiXAxisString  = "Precipitation (m/year)"

  1058   max_bar = new((/2,nx/),graphic)

  1059   min_bar = new((/2,nx/),graphic)

  1060   max_cap = new((/2,nx/),graphic)

  1061   min_cap = new((/2,nx/),graphic)

  1063   lnres = True

  1064   line_colors = (/"brown","blue"/)

  1066 ;**************************************************************************

  1067 ;(B)-1 histogram plot, ob 81 site

  1068 ;**************************************************************************

  1070   plot_name = "histogram_ob_81"

  1071   title     = "Observed 81 site"

  1072   resh@tiMainString  = title

  1074   wks   = gsn_open_wks (plot_type,plot_name)

  1076   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

  1078 ;-------------------------------

  1079 ;Attach the vertical bar and the horizontal cap line

  1081   delete (x1)

  1082   delete (x2)

  1083   delete (y1)

  1084   delete (y2)

  1086   do nd=0,0

  1087     lnres@gsLineColor = line_colors(nd)

  1088     do i=0,nx-1

  1090       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1091          .not.ismissing(mx_yvalues(nd,i))) then

  1092 ;

  1093 ; Attach the vertical bar, both above and below the marker.

  1094 ;

  1095         x1 = xvalues(nd,i)

  1096         y1 = yvalues(nd,i)

  1097         y2 = mn_yvalues(nd,i)

  1098         plx = (/x1,x1/)

  1099         ply = (/y1,y2/)

  1100         min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1102         y2 = mx_yvalues(nd,i)

  1103         plx = (/x1,x1/)

  1104         ply = (/y1,y2/)

  1105         max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1106 ;

  1107 ; Attach the horizontal cap line, both above and below the marker.

  1108 ;

  1109         x1 = xvalues(nd,i) - dx4

  1110         x2 = xvalues(nd,i) + dx4

  1111         y1 = mn_yvalues(nd,i)

  1112         plx = (/x1,x2/)

  1113         ply = (/y1,y1/)

  1114         min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1116         y1 = mx_yvalues(nd,i)

  1117         plx = (/x1,x2/)

  1118         ply = (/y1,y1/)

  1119         max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

  1120       end if

  1121     end do

  1122   end do

  1124   draw(xy)

  1125   frame(wks)

  1126   clear (wks)

  1128  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1129  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1130  system("rm "+plot_name+"-*."+plot_type_new)

  1131  system("rm "+plot_name+"."+plot_type)

  1133 ;****************************************************************************

  1134 ;(B)-2 histogram plot, model vs ob 81 site

  1135 ;****************************************************************************

  1137   plot_name = "histogram_model_vs_ob_81"

  1138   title     = model_name+ " vs Observed 81 site"

  1139   resh@tiMainString  = title

  1141   wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

  1143 ;-----------------------------

  1144 ; Add a boxed legend using the more simple method, which won't have

  1145 ; vertical lines going through the markers.

  1147   resh@pmLegendDisplayMode    = "Always"

  1148 ; resh@pmLegendWidthF         = 0.1

  1149   resh@pmLegendWidthF         = 0.08

  1150   resh@pmLegendHeightF        = 0.05

  1151   resh@pmLegendOrthogonalPosF = -1.17

  1152 ; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

  1153 ; resh@pmLegendParallelPosF   =  0.18

  1154   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1156 ; resh@lgPerimOn              = False

  1157   resh@lgLabelFontHeightF     = 0.015

  1158   resh@xyExplicitLegendLabels = (/"observed",model_name/)

  1159 ;-----------------------------

  1160   tRes  = True

  1161   tRes@txFontHeightF = 0.025

  1163   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"

  1165   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

  1167   xy = gsn_csm_xy(wks,xvalues,yvalues,resh)

  1168 ;-------------------------------

  1169 ;Attach the vertical bar and the horizontal cap line

  1171   do nd=0,1

  1172     lnres@gsLineColor = line_colors(nd)

  1173     do i=0,nx-1

  1175       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1176          .not.ismissing(mx_yvalues(nd,i))) then

  1177 ;

  1178 ; Attach the vertical bar, both above and below the marker.

  1179 ;

  1180         x1 = xvalues(nd,i)

  1181         y1 = yvalues(nd,i)

  1182         y2 = mn_yvalues(nd,i)

  1183         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1185         y2 = mx_yvalues(nd,i)

  1186         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1187 ;

  1188 ; Attach the horizontal cap line, both above and below the marker.

  1189 ;

  1190         x1 = xvalues(nd,i) - dx4

  1191         x2 = xvalues(nd,i) + dx4

  1192         y1 = mn_yvalues(nd,i)

  1193         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1195         y1 = mx_yvalues(nd,i)

  1196         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1197       end if

  1198     end do

  1199   end do

  1201   draw(xy)

  1202   frame(wks)

  1203   clear(wks)

  1205  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1206  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1207  system("rm "+plot_name+"-*."+plot_type_new)

  1208  system("rm "+plot_name+"."+plot_type)

  1210  delete (RAIN1_1D)

  1211  delete (RAIN2_1D)

  1212  delete (NPP1_1D)

  1213  delete (NPP2_1D)

  1214 ;delete (range)

  1215  delete (xvalues)

  1216  delete (yvalues)

  1217  delete (mn_yvalues)

  1218  delete (mx_yvalues)

  1219  delete (good)

  1220  delete (max_bar)

  1221  delete (min_bar)

  1222  delete (max_cap)

  1223  delete (min_cap)

  1225 ;**************************************************************************

  1226 ;(B) histogram 933

  1227 ;**************************************************************************

  1229 ; get data

  1230   RAIN1_1D = ndtooned(rain933)

  1231   RAIN2_1D = ndtooned(rainmod933)

  1232   NPP1_1D  = ndtooned(npp933)

  1233   NPP2_1D  = ndtooned(nppmod933)

  1235 ; number of bin

  1236   nx = 10

  1238   xvalues      = new((/2,nx/),float)

  1239   yvalues      = new((/2,nx/),float)

  1240   mn_yvalues   = new((/2,nx/),float)

  1241   mx_yvalues   = new((/2,nx/),float)

  1242   dx4          = new((/1/),float)

  1244   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

  1245          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

  1247 ;----------------------------------------

  1248 ;compute correlation coeff and M score

  1250  u = yvalues(0,:)

  1251  v = yvalues(1,:)

  1253  good = ind(.not.ismissing(u) .and. .not.ismissing(v))

  1254  uu = u(good)

  1255  vv = v(good)

  1257  ccr933h = esccr(uu,vv,0)

  1258 ;print (ccr933h)

  1260 ;new eq

  1261  bias  = sum(abs(vv-uu)/(abs(vv)+abs(uu)))

  1262  M933h = (1.- (bias/dimsizes(uu)))*5.

  1264  M_npp_H933 = sprintf("%.2f", M933h)

  1265  system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new)

  1266  system("mv -f "+html_new+" "+html_name)

  1267  print (M_npp_H933)

  1269  delete (u)

  1270  delete (v)

  1271  delete (uu)

  1272  delete (vv)

  1273 ;----------------------------------------------------------------------

  1274 ; histogram res

  1275   delete (resh)

  1276   resh                = True

  1277   resh@gsnMaximize    = True

  1278   resh@gsnDraw        = False

  1279   resh@gsnFrame       = False

  1280   resh@xyMarkLineMode = "Markers"

  1281   resh@xyMarkerSizeF  = 0.014

  1282   resh@xyMarker       = 16

  1283   resh@xyMarkerColors = (/"Brown","Blue"/)

  1284   resh@trYMinF        = min(mn_yvalues) - 10.

  1285   resh@trYMaxF        = max(mx_yvalues) + 10.

  1287   resh@tiYAxisString  = "NPP (g C/m2/year)"

  1288   resh@tiXAxisString  = "Precipitation (m/year)"

  1290   max_bar = new((/2,nx/),graphic)

  1291   min_bar = new((/2,nx/),graphic)

  1292   max_cap = new((/2,nx/),graphic)

  1293   min_cap = new((/2,nx/),graphic)

  1295   lnres = True

  1296   line_colors = (/"brown","blue"/)

  1298 ;**************************************************************************

  1299 ;(B)-3 histogram plot, ob 933 site

  1300 ;**************************************************************************

  1302   plot_name = "histogram_ob_933"

  1303   title     = "Observed 933 site"

  1304   resh@tiMainString  = title

  1306   wks   = gsn_open_wks (plot_type,plot_name)

  1308   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

  1310 ;-------------------------------

  1311 ;Attach the vertical bar and the horizontal cap line

  1313   do nd=0,0

  1314     lnres@gsLineColor = line_colors(nd)

  1315     do i=0,nx-1

  1317       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1318          .not.ismissing(mx_yvalues(nd,i))) then

  1320 ; Attach the vertical bar, both above and below the marker.

  1322         x1 = xvalues(nd,i)

  1323         y1 = yvalues(nd,i)

  1324         y2 = mn_yvalues(nd,i)

  1325         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1327         y2 = mx_yvalues(nd,i)

  1328         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1330 ; Attach the horizontal cap line, both above and below the marker.

  1332         x1 = xvalues(nd,i) - dx4

  1333         x2 = xvalues(nd,i) + dx4

  1334         y1 = mn_yvalues(nd,i)

  1335         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1337         y1 = mx_yvalues(nd,i)

  1338         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1339       end if

  1340     end do

  1341   end do

  1343   draw(xy)

  1344   frame(wks)

  1345   delete (xy)

  1346   clear (wks)

  1348  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1349  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1350  system("rm "+plot_name+"-*."+plot_type_new)

  1351  system("rm "+plot_name+"."+plot_type)

  1353 ;****************************************************************************

  1354 ;(B)-4 histogram plot, model vs ob 933 site

  1355 ;****************************************************************************

  1357   plot_name = "histogram_model_vs_ob_933"

  1358   title     = model_name+ " vs Observed 933 site"

  1359   resh@tiMainString  = title

  1361   wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

  1363 ;-----------------------------

  1364 ; Add a boxed legend using the more simple method, which won't have

  1365 ; vertical lines going through the markers.

  1367   resh@pmLegendDisplayMode    = "Always"

  1368 ; resh@pmLegendWidthF         = 0.1

  1369   resh@pmLegendWidthF         = 0.08

  1370   resh@pmLegendHeightF        = 0.05

  1371   resh@pmLegendOrthogonalPosF = -1.17

  1372 ; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

  1373 ; resh@pmLegendParallelPosF   =  0.18

  1374   resh@pmLegendParallelPosF   =  0.88  ;(rightward)

  1376 ; resh@lgPerimOn              = False

  1377   resh@lgLabelFontHeightF     = 0.015

  1378   resh@xyExplicitLegendLabels = (/"observed",model_name/)

  1379 ;-----------------------------

  1380   tRes  = True

  1381   tRes@txFontHeightF = 0.025

  1383   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"

  1385   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

  1387   xy = gsn_csm_xy(wks,xvalues,yvalues,resh)

  1388 ;-------------------------------

  1389 ;Attach the vertical bar and the horizontal cap line

  1391   do nd=0,1

  1392     lnres@gsLineColor = line_colors(nd)

  1393     do i=0,nx-1

  1395       if(.not.ismissing(mn_yvalues(nd,i)).and. \

  1396          .not.ismissing(mx_yvalues(nd,i))) then

  1397 ;

  1398 ; Attach the vertical bar, both above and below the marker.

  1399 ;

  1400         x1 = xvalues(nd,i)

  1401         y1 = yvalues(nd,i)

  1402         y2 = mn_yvalues(nd,i)

  1403         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1405         y2 = mx_yvalues(nd,i)

  1406         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

  1407 ;

  1408 ; Attach the horizontal cap line, both above and below the marker.

  1409 ;

  1410         x1 = xvalues(nd,i) - dx4

  1411         x2 = xvalues(nd,i) + dx4

  1412         y1 = mn_yvalues(nd,i)

  1413         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1415         y1 = mx_yvalues(nd,i)

  1416         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

  1417       end if

  1418     end do

  1419   end do

  1421   draw(xy)

  1422   frame(wks)

  1423   delete(xy)

  1424   clear(wks)

  1426  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1427  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1428  system("rm "+plot_name+"-*."+plot_type_new)

  1429  system("rm "+plot_name+"."+plot_type)

  1431 ;***************************************************************************

  1432 ;(C) global contour

  1433 ;***************************************************************************

  1435 ;res

  1436   resg                     = True             ; Use plot options

  1437   resg@cnFillOn            = True             ; Turn on color fill

  1438   resg@gsnSpreadColors     = True             ; use full colormap

  1439 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

  1440 ; resg@lbLabelAutoStride   = True

  1441   resg@cnLinesOn           = False            ; Turn off contourn lines

  1442   resg@mpFillOn            = False            ; Turn off map fill

  1444   resg@gsnSpreadColors      = True            ; use full colormap

  1445   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  1446   resg@cnMinLevelValF       = 0.              ; Min level

  1447   resg@cnMaxLevelValF       = 2200.           ; Max level

  1448   resg@cnLevelSpacingF      = 200.            ; interval

  1450 ;****************************************************************************

  1451 ;(C)-1 global contour plot, ob

  1452 ;****************************************************************************

  1454   delta = 0.00000000001

  1455   nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)

  1457   plot_name = "global_ob"

  1458   title     = "Observed MODIS MOD 17"

  1459   resg@tiMainString  = title

  1461   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1462   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1464   plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1466   frame(wks)

  1467   clear (wks)

  1469  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1470  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1471  system("rm "+plot_name+"-*."+plot_type_new)

  1472  system("rm "+plot_name+"."+plot_type)

  1474 ;****************************************************************************

  1475 ;(C)-2 global contour plot, model

  1476 ;****************************************************************************

  1478   plot_name = "global_model"

  1479   title     = "Model "+ model_name

  1480   resg@tiMainString  = title

  1482   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1483   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1485   plot = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1487   frame(wks)

  1488   clear (wks)

  1490  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1491  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1492  system("rm "+plot_name+"-*."+plot_type_new)

  1493  system("rm "+plot_name+"."+plot_type)

  1495 ;****************************************************************************

  1496 ;(C)-3 global contour plot, model vs ob

  1497 ;****************************************************************************

  1499   plot_name = "global_model_vs_ob"

  1501   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1502   gsn_define_colormap(wks,"gui_default")     ; choose colormap

  1504   delete (plot)

  1505   plot=new(3,graphic)                        ; create graphic array

  1507   resg@gsnFrame             = False          ; Do not draw plot

  1508   resg@gsnDraw              = False          ; Do not advance frame

  1510 ; compute correlation coef and M score

  1512   uu1 = ndtooned(nppmod)

  1513   vv1 = ndtooned(nppglobe)

  1515   delete (good)

  1516   good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))

  1518   ug = uu1(good)

  1519   vg = vv1(good)

  1521   ccrG = esccr(ug,vg,0)

  1522 ; print (ccrG)

  1524   MG   = (ccrG*ccrG)* 5.0

  1526   M_npp_G = sprintf("%.2f", MG)

  1527   system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new)

  1528   system("mv -f "+html_new+" "+html_name)

  1529   print (M_npp_G)

  1531 ; plot correlation coef

  1533   gRes  = True

  1534   gRes@txFontHeightF = 0.02

  1535   gRes@txAngleF = 90

  1537   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"

  1539   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

  1540 ;--------------------------------------------------------------------

  1541 ;(a) ob

  1543   title     = "Observed MODIS MOD 17"

  1544   resg@tiMainString  = title

  1546   plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)

  1548 ;(b) model

  1550   title     = "Model "+ model_name

  1551   resg@tiMainString  = title

  1553   plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)

  1555 ;(c) model-ob

  1557   zz = nppmod

  1558   zz = nppmod - nppglobe

  1559   title = "Model_"+model_name+" - Observed"

  1561   resg@cnMinLevelValF  = -500           ; Min level

  1562   resg@cnMaxLevelValF  =  500.          ; Max level

  1563   resg@cnLevelSpacingF =  50.           ; interval

  1564   resg@tiMainString    = title

  1566   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

  1568   pres                            = True        ; panel plot mods desired

  1569 ; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

  1570                                                 ; indiv. plots in panel

  1571   pres@gsnMaximize                = True        ; fill the page

  1573   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

  1575   frame (wks)

  1576   clear (wks)

  1577   delete (plot)

  1579  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1580  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1581  system("rm "+plot_name+"-*."+plot_type_new)

  1582  system("rm "+plot_name+"."+plot_type)

  1584 ;***************************************************************************

  1585 ;(D)-1 zonal line plot, ob

  1586 ;***************************************************************************

  1588   vv     = zonalAve(nppglobe)

  1589   vv@long_name = nppglobe@long_name

  1591   plot_name = "zonal_ob"

  1592   title     = "Observed MODIS MOD 17"

  1594   resz = True

  1595   resz@tiMainString  = title

  1597   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  1599   plot  = gsn_csm_xy (wks,ym,vv,resz)

  1601   frame(wks)

  1602   clear (wks)

  1604  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1605  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1606  system("rm "+plot_name+"-*."+plot_type_new)

  1607  system("rm "+plot_name+"."+plot_type)

  1609 ;****************************************************************************

  1610 ;(D)-2 zonal line plot, model vs ob

  1611 ;****************************************************************************

  1613   s  = new ((/2,dimsizes(ym)/), float)

  1615   s(0,:) = zonalAve(nppglobe)

  1616   s(1,:) = zonalAve(nppmod)

  1618   s@long_name = nppglobe@long_name

  1619 ;-------------------------------------------

  1620 ; compute correlation coef and M score

  1622   ccrZ = esccr(s(1,:), s(0,:),0)

  1623 ; print (ccrZ)

  1625   MZ   = (ccrZ*ccrZ)* 5.0

  1627   M_npp_Z = sprintf("%.2f", MZ)

  1628   system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new)

  1629   system("mv -f "+html_new+" "+html_name)

  1630   print (M_npp_Z)

  1631 ;-------------------------------------------

  1632   plot_name = "zonal_model_vs_ob"

  1633   title     = "Zonal Average"

  1634   resz@tiMainString  = title

  1636   wks = gsn_open_wks (plot_type,plot_name)

  1638 ; resz@vpHeightF          = 0.4               ; change aspect ratio of plot

  1639 ; resz@vpWidthF           = 0.7

  1641   resz@xyMonoLineColor    = "False"           ; want colored lines

  1642   resz@xyLineColors       = (/"black","red"/) ; colors chosen

  1643 ; resz@xyLineThicknesses  = (/3.,3./)         ; line thicknesses

  1644   resz@xyLineThicknesses  = (/2.,2./)         ; line thicknesses

  1645   resz@xyDashPatterns     = (/0.,0./)         ; make all lines solid

  1647   resz@tiYAxisString      = s@long_name       ; add a axis title

  1648   resz@txFontHeightF      = 0.0195            ; change title font heights

  1650 ; Legent

  1651   resz@pmLegendDisplayMode    = "Always"      ; turn on legend

  1652   resz@pmLegendSide           = "Top"         ; Change location of

  1653 ; resz@pmLegendParallelPosF   = .45           ; move units right

  1654   resz@pmLegendParallelPosF   = .82           ; move units right

  1655   resz@pmLegendOrthogonalPosF = -0.4          ; move units down

  1657   resz@pmLegendWidthF         = 0.10          ; Change width and

  1658   resz@pmLegendHeightF        = 0.10          ; height of legend.

  1659   resz@lgLabelFontHeightF     = .02           ; change font height

  1660 ; resz@lgTitleOn              = True          ; turn on legend title

  1661 ; resz@lgTitleString          = "Example"     ; create legend title

  1662 ; resz@lgTitleFontHeightF     = .025          ; font of legend title

  1663   resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels

  1664 ;--------------------------------------------------------------------

  1665   zRes  = True

  1666   zRes@txFontHeightF = 0.025

  1668   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"

  1670   gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)

  1671 ;--------------------------------------------------------------------

  1673   plot  = gsn_csm_xy (wks,ym,s,resz)

  1675   frame(wks)

  1676   clear (wks)

  1678  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  1679  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  1680  system("rm "+plot_name+"-*."+plot_type_new)

  1681  system("rm "+plot_name+"."+plot_type)

  1683 ;***************************************************************************

  1684 ; tar up output plots

  1685 ;***************************************************************************

  1687   temp_name = "npp." + model_name

  1688   system("mkdir -p " + temp_name)

  1689   system("cp "+ html_name + " " +temp_name+"/table.html")

  1690   system("mv *.png " + temp_name)

  1691   system("tar cf "+ temp_name +".tar " + temp_name)

  1693 ;***************************************************************************

  1694 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--