npp/15.scatter_bias.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:a6a37ba36f8d
       
     1 ; ***********************************************
       
     2 ; combine all scatter
       
     3 ; ***********************************************
       
     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
     7 ;************************************************
       
     8 begin
       
     9 
       
    10  plot_type     = "ps"
       
    11  plot_type_new = "png"
       
    12 
       
    13 ;************************************************
       
    14 ; read in ob data
       
    15 ;************************************************
       
    16 
       
    17 ;(A) plot ob scatter_81
       
    18 
       
    19  dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
       
    20  fil81 = "data.81.nc"
       
    21  f81   = addfile (dir81+fil81,"r")
       
    22 
       
    23  x81   = f81->SITE_ID  
       
    24  y81   = f81->TNPP_C
       
    25  xo81  = f81->LONG_DD  
       
    26  yo81  = f81->LAT_DD
       
    27 
       
    28  x81@long_name = "SITE_ID"
       
    29  y81@long_name = "NPP (gC/m2/year)"
       
    30    
       
    31  plot_name = "scatter_ob_81"
       
    32  title     = "Observed 81 sites"
       
    33 
       
    34  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
       
    35  res                   = True                  ; plot mods desired
       
    36  res@tiMainString      = title                 ; add title
       
    37  res@xyMarkLineModes   = "Markers"             ; choose which have markers
       
    38  res@xyMarkers         =  16                   ; choose type of marker  
       
    39  res@xyMarkerColor     = "red"                 ; Marker color
       
    40  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
       
    41  res@tmLabelAutoStride = True                  ; nice tick mark labels
       
    42 
       
    43  plot  = gsn_csm_xy (wks,x81,y81,res)              ; create plot
       
    44  frame(wks)
       
    45 
       
    46  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
       
    47  system("rm "+plot_name+"."+plot_type)
       
    48  system("rm "+plot_name+"-1."+plot_type_new)
       
    49  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
       
    50 
       
    51  clear (wks)
       
    52  
       
    53 ;(B) plot ob scatter_933
       
    54 
       
    55  dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
       
    56  fil933 = "data.933.nc"
       
    57  f933   = addfile (dir933+fil933,"r")
       
    58 
       
    59  x933  = f933->SITE_ID  
       
    60  y933  = f933->TNPP_C
       
    61  xo933 = f933->LONG_DD  
       
    62  yo933 = f933->LAT_DD 
       
    63 
       
    64  x933@long_name = "SITE_ID"
       
    65  y933@long_name = "NPP (gC/m2/year)"
       
    66 
       
    67  plot_name = "scatter_ob_933"
       
    68  title     = "Observed 933 sites"
       
    69 
       
    70  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
       
    71  res                   = True                  ; plot mods desired
       
    72  res@tiMainString      = title                 ; add title
       
    73  res@xyMarkLineModes   = "Markers"             ; choose which have markers
       
    74  res@xyMarkers         =  16                   ; choose type of marker  
       
    75  res@xyMarkerColor     = "red"                 ; Marker color
       
    76  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
       
    77  res@tmLabelAutoStride = True                  ; nice tick mark labels
       
    78 
       
    79  plot  = gsn_csm_xy (wks,x933,y933,res)              ; create plot
       
    80  frame(wks)
       
    81 
       
    82  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
       
    83  system("rm "+plot_name+"."+plot_type)
       
    84  system("rm "+plot_name+"-1."+plot_type_new)
       
    85  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
       
    86 
       
    87  clear (wks)
       
    88 ;************************************************
       
    89 ; read in model data
       
    90 ;************************************************
       
    91 
       
    92 ;(C) plot model scatter_81
       
    93 
       
    94  model_name = "b30.061n"
       
    95 
       
    96  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    97  film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
       
    98  fm   = addfile (dirm+film,"r")
       
    99   
       
   100  ymod = fm->NPP
       
   101  xm   = fm->lon  
       
   102  ym   = fm->lat
       
   103 
       
   104  nx = dimsizes(xo81)
       
   105  do i= 0,nx-1
       
   106     if (xo81(i) .lt. 0.) then
       
   107         xo81(i) = xo81(i)+ 360.
       
   108     end if
       
   109  end do
       
   110  print (xo81)
       
   111 
       
   112  sec_to_year = 86400.*365.
       
   113 
       
   114  ymod81 = linint2_points(xm,ym,ymod(0,:,:),True,xo81,yo81,0) * sec_to_year
       
   115  xmod81 = x81
       
   116 
       
   117  ymod81@long_name = "NPP (gC/m2/year)"
       
   118  xmod81@long_name = "SITE_ID"
       
   119 print (ymod81)
       
   120 print (xmod81)
       
   121   
       
   122  plot_name = "scatter_model_81"
       
   123  title     = "Model "+ model_name +" 81 sites"
       
   124 
       
   125  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
       
   126  res                   = True                  ; plot mods desired
       
   127  res@tiMainString      = title                 ; add title
       
   128  res@xyMarkLineModes   = "Markers"             ; choose which have markers
       
   129  res@xyMarkers         =  16                   ; choose type of marker  
       
   130  res@xyMarkerColor     = "red"                 ; Marker color
       
   131  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
       
   132  res@tmLabelAutoStride = True                  ; nice tick mark labels
       
   133 
       
   134  plot  = gsn_csm_xy (wks,xmod81,ymod81,res)    ; create plot
       
   135  frame(wks)
       
   136 
       
   137  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
       
   138  system("rm "+plot_name+"."+plot_type)
       
   139  system("rm "+plot_name+"-1."+plot_type_new)
       
   140  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
       
   141 
       
   142  clear (wks)
       
   143  
       
   144 ;(D) plot model scatter_933
       
   145 
       
   146  model_name = "b30.061n"
       
   147 
       
   148  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
   149  film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
       
   150  fm   = addfile (dirm+film,"r")
       
   151   
       
   152  ymod = fm->NPP
       
   153  xm   = fm->lon  
       
   154  ym   = fm->lat
       
   155 
       
   156  nx = dimsizes(xo933)
       
   157  do i= 0,nx-1
       
   158     if (xo933(i) .lt. 0.) then
       
   159         xo933(i) = xo933(i)+ 360.
       
   160     end if
       
   161  end do
       
   162  print (xo933)
       
   163 
       
   164  sec_to_year = 86400.*365.
       
   165 
       
   166  ymod933 = linint2_points(xm,ym,ymod(0,:,:),True,xo933,yo933,0) * sec_to_year
       
   167  xmod933 = x933
       
   168 
       
   169  ymod933@long_name = "NPP (gC/m2/year)"
       
   170  xmod933@long_name = "SITE_ID"
       
   171 ;print (ymod933)
       
   172 ;print (xmod933)
       
   173   
       
   174  plot_name = "scatter_model_933"
       
   175  title     = "Model "+ model_name +" 933 sites"
       
   176 
       
   177  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
       
   178  res                   = True                  ; plot mods desired
       
   179  res@tiMainString      = title                 ; add title
       
   180  res@xyMarkLineModes   = "Markers"             ; choose which have markers
       
   181  res@xyMarkers         =  16                   ; choose type of marker  
       
   182  res@xyMarkerColor     = "red"                 ; Marker color
       
   183  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
       
   184  res@tmLabelAutoStride = True                  ; nice tick mark labels
       
   185 
       
   186  plot  = gsn_csm_xy (wks,xmod933,ymod933,res)    ; create plot
       
   187  frame(wks)
       
   188 
       
   189  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
       
   190  system("rm "+plot_name+"."+plot_type)
       
   191  system("rm "+plot_name+"-1."+plot_type_new)
       
   192  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
       
   193 
       
   194  clear (wks)
       
   195 
       
   196 ;(E) scatter model vs ob 81
       
   197 
       
   198  ymod81@long_name = "NPP (gC/m2/year)"
       
   199  y81@long_name    = "NPP (gC/m2/year)"
       
   200 ;print (ymod81)
       
   201 ;print (y81)
       
   202   
       
   203  plot_name = "scatter_model_vs_ob_81"
       
   204  title     = model_name +" vs ob 81 sites"
       
   205 
       
   206  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
       
   207  res                   = True                  ; plot mods desired
       
   208  res@tiMainString      = title                 ; add title
       
   209  res@xyMarkLineModes   = "Markers"             ; choose which have markers
       
   210  res@xyMarkers         =  16                   ; choose type of marker  
       
   211  res@xyMarkerColor     = "red"                 ; Marker color
       
   212  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
       
   213  res@tmLabelAutoStride = True                  ; nice tick mark labels
       
   214 ;-------------------------------
       
   215 ;compute correlation coef. and M
       
   216  ccr81 = esccr(ymod81,y81,0)
       
   217  print (ccr81)
       
   218 
       
   219 ;new eq
       
   220  bias = sum(abs(ymod81-y81)/(ymod81+y81))
       
   221  M81  = (1. - (bias/dimsizes(y81)))*5.
       
   222  print (M81)
       
   223 
       
   224  tRes  = True
       
   225  tRes@txFontHeightF = 0.025
       
   226 
       
   227  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
       
   228 
       
   229  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
       
   230 ;--------------------------------
       
   231  plot  = gsn_csm_xy (wks,y81,ymod81,res)       ; create plot
       
   232  frame(wks)
       
   233 
       
   234  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
       
   235  system("rm "+plot_name+"."+plot_type)
       
   236  system("rm "+plot_name+"-1."+plot_type_new)
       
   237  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
       
   238 
       
   239  clear (wks)
       
   240  
       
   241 ;(F) scatter model vs ob 933
       
   242  
       
   243  ymod933@long_name = "NPP (gC/m2/year)"
       
   244  y933@long_name    = "NPP (gC/m2/year)"
       
   245 ;print (ymod33)
       
   246 ;print (y933)
       
   247   
       
   248  plot_name = "scatter_model_vs_ob_933"
       
   249  title     = model_name +" vs ob 933 sites"
       
   250 
       
   251  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation
       
   252  res                   = True                  ; plot mods desired
       
   253  res@tiMainString      = title                 ; add title
       
   254  res@xyMarkLineModes   = "Markers"             ; choose which have markers
       
   255  res@xyMarkers         =  16                   ; choose type of marker  
       
   256  res@xyMarkerColor     = "red"                 ; Marker color
       
   257  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)
       
   258  res@tmLabelAutoStride = True                  ; nice tick mark labels
       
   259 ;-------------------------------
       
   260 ;compute correlation coef. and M
       
   261  ccr933 = esccr(ymod933,y933,0)
       
   262  print (ccr933)
       
   263 
       
   264 ;new eq
       
   265  bias = sum(abs(ymod933-y933)/(ymod933+y933))
       
   266  M933 = (1. - (bias/dimsizes(y933)))*5.
       
   267  print (M933)
       
   268 
       
   269  tRes  = True
       
   270  tRes@txFontHeightF = 0.025
       
   271 
       
   272  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
       
   273 
       
   274  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
       
   275 ;--------------------------------
       
   276  plot  = gsn_csm_xy (wks,y933,ymod933,res)    ; create plot
       
   277  frame(wks)
       
   278 
       
   279  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
       
   280  system("rm "+plot_name+"."+plot_type)
       
   281  system("rm "+plot_name+"-1."+plot_type_new)
       
   282  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
       
   283 
       
   284  clear (wks)
       
   285 end