clamp: npp/24.histogram_ob.ncl@4be95183fbcd

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

     2 ; histogram normalized by rain and compute correleration

     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"

     8 procedure pminmax(data:numeric,name:string)

     9 begin

    10   print ("min/max " + name + " = " + min(data) + "/" + max(data))

    11   if(isatt(data,"units")) then

    12     print (name + " units = " + data@units)

    13   end if

    14 end

    16 ;

    17 ; Main code.

    18 ;

    19 begin

    20   data_types    = (/ "Obs",        "Model"       /)

    21   data_names    = (/ "data.81.nc" , "i01.03cn_1545-1569_ANN_climo.nc" /)

    22 ; data_names    = (/ "data.933.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)

    23 ; data_names    = (/ "data.81.nc" , "i01.04casa_1605-1629_ANN_climo.nc" /)

    24 ; data_names    = (/ "data.933.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)

    25   filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)

    26   ndata_types   = dimsizes(data_types)

    28   data_file_obs = addfile(data_names(0),"r")    ; Open obs file

    29   data_file_mod = addfile(data_names(1),"r")    ; Open model file

    31 ;************************************************

    32 ; read in data: observed

    33 ;************************************************

    34  RAIN1 = tofloat(data_file_obs->PREC_ANN) ; for data.81

    35 ;RAIN1 = data_file_obs->PREC              ; for data.933

    36  NPP1  = data_file_obs->TNPP_C

    37  xo    = data_file_obs->LONG_DD

    38  yo    = data_file_obs->LAT_DD

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

    41  print (xo)

    42  nx = dimsizes(xo)

    43  do i= 0,nx-1

    44     if (xo(i) .lt. 0.) then

    45         xo(i) = xo(i)+ 360.

    46     end if

    47  end do

    48  print (xo)

    50 ;************************************************

    51 ; read in data: model

    52 ;************************************************

    53  ai    = data_file_mod->RAIN

    54  bi    = data_file_mod->NPP

    55  xi    = data_file_mod->lon

    56  yi    = data_file_mod->lat

    58 ;************************************************

    59 ; interpolate from model grid to observed grid

    60 ;************************************************

    61  RAIN2 = linint2_points(xi,yi,ai,True,xo,yo,0)

    62  NPP2  = linint2_points(xi,yi,bi,True,xo,yo,0)

    64 ;************************************************

    65 ; convert unit

    66 ;************************************************

    67 ; Units for these four variables are:

    68 ;

    69 ; RAIN1 : mm/year

    70 ; RAIN2 : mm/s

    71 ; NPP1  : g C/m^2/year

    72 ; NPP2  : g C/m^2/s

    73 ;

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

    75 ;

    76   nsec_per_year = 60*60*24*365                  ; # seconds per year

    78 ; Do the necessary conversions.

    79   RAIN1 = RAIN1 / 1000.

    80   RAIN2 = (RAIN2/ 1000.) * nsec_per_year

    81   NPP2     = NPP2 * nsec_per_year

    83 ; Redo the units.

    84   RAIN1@units = "m/yr"

    85   RAIN2@units = "m/yr"

    86   NPP1@units  = "gC/m^2/yr"

    87   NPP2@units  = "gC/m^2/yr"

    89 ;************************************************

    90 ; print min/max and unit

    91 ;************************************************

    92   pminmax(RAIN1,"RAIN1")

    93   pminmax(RAIN2,"RAIN2")

    94   pminmax(NPP1,"NPP1")

    95   pminmax(NPP2,"NPP2")

    97   RAIN1_1D = ndtooned(RAIN1)

    98   RAIN2_1D = ndtooned(RAIN2)

    99   NPP1_1D  = ndtooned(NPP1)

   100   NPP2_1D  = ndtooned(NPP2)

   101 ;

   102 ; Calculate some "nice" bins for binning the data in equally spaced

   103 ; ranges.

   104 ;

   105   nbins       = 15     ; Number of bins to use.

   107   nicevals    = nice_mnmxintvl(min(RAIN2_1D),max(RAIN2_1D),nbins,True)

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

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

   110 ;

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

   112 ; particular range, and then take an average.

   113 ;

   114   nr      = dimsizes(range)

   115   nx      = nr-1

   116   xvalues     = new((/2,nx/),typeof(RAIN2_1D))

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

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

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

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

   121   yvalues      = new((/2,nx/),typeof(RAIN2_1D))

   122   mn_yvalues   = new((/2,nx/),typeof(RAIN2_1D))

   123   mx_yvalues   = new((/2,nx/),typeof(RAIN2_1D))

   125   do nd=0,1

   126 ;

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

   128 ;

   129     if(nd.eq.0) then

   130       data     = RAIN1_1D

   131       npp_data = NPP1_1D

   132     else

   133       data     = RAIN2_1D

   134       npp_data = NPP2_1D

   135     end if

   136 ;

   137 ; Loop through each range and check for values.

   138 ;

   139     do i=0,nr-2

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

   141 ;        print("")

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

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

   144       else

   145 ;        print("")

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

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

   148       end if

   149 ;

   150 ; Calculate average, and get min and max.

   151 ;

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

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

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

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

   156         count = dimsizes(idx)

   157       else

   158         count            = 0

   159         yvalues(nd,i)    = yvalues@_FillValue

   160         mn_yvalues(nd,i) = yvalues@_FillValue

   161         mx_yvalues(nd,i) = yvalues@_FillValue

   162       end if

   163 ;

   164 ; Print out information.

   165 ;

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

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

   169 ;

   170 ; Clean up for next time in loop.

   171 ;

   172       delete(idx)

   173     end do

   174     delete(data)

   175     delete(npp_data)

   176   end do

   178 ;

   179 ; Start the graphics.

   180 ;

   181   wks = gsn_open_wks("png","hist")

   183   res             = True

   184   res@gsnMaximize = True

   185   res@gsnDraw     = False

   186   res@gsnFrame    = False

   187   res@xyMarkLineMode = "Markers"

   188   res@xyMarkerSizeF   = 0.014

   189   res@xyMarker       = 16

   190   res@xyMarkerColors = (/"Brown","Blue"/)

   191   res@trYMinF        = min(mn_yvalues) - 10.

   192   res@trYMaxF        = max(mx_yvalues) + 10.

   194   res@tiMainString  = "Observed 81 sites"

   195 ; res@tiMainString  = "Observed 933 sites"

   196   res@tiYAxisString = "NPP (g C/m2/year)"

   197   res@tiXAxisString = "Precipitation (m/year)"

   199 ;

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

   201 ; vertical lines going through the markers.

   202 ;

   203 ; res@pmLegendDisplayMode    = "Always"

   204 ;;res@pmLegendWidthF         = 0.1

   205 ; res@pmLegendWidthF         = 0.08

   206 ; res@pmLegendHeightF        = 0.05

   207 ; res@pmLegendOrthogonalPosF = -1.17

   208 ;;res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   209 ;;res@pmLegendParallelPosF   =  0.18

   210 ; res@pmLegendParallelPosF   =  0.23  ;(rightward)

   212 ; res@lgPerimOn              = False

   213 ; res@lgLabelFontHeightF     = 0.015

   214 ; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)

   215 ;;res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)

   217   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)

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

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

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

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

   224   lnres = True

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

   227 ; do nd=0,1

   228   do nd=0,0

   229     lnres@gsLineColor = line_colors(nd)

   230     do i=0,nx-1

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

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

   234 ;

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

   236 ;

   237         x1 = xvalues(nd,i)

   238         y1 = yvalues(nd,i)

   239         y2 = mn_yvalues(nd,i)

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

   242         y2 = mx_yvalues(nd,i)

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

   244 ;

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

   246 ;

   247         x1 = xvalues(nd,i) - dx4

   248         x2 = xvalues(nd,i) + dx4

   249         y1 = mn_yvalues(nd,i)

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

   252         y1 = mx_yvalues(nd,i)

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

   254       end if

   255     end do

   256   end do

   258 ;

   259 ; Here's how to do the legend by hand.

   260 ;

   261 ; mkres               = True    ; Marker resources

   262 ; txres               = True    ; Text resources

   263 ; mkres@gsMarkerIndex = 16

   264 ; mkres@gsMarkerSizeF = 0.02

   265 ; txres@txFontHeightF = 0.02

   266 ; txres@txJust        = "CenterLeft"

   267 ;

   268 ; Change these values if you want to move the marker legend location.

   269 ; These values are in the same data space as the plot.

   270 ;

   271 ; xlg1_cen = 0.2

   272 ; ylg1_cen = 900.

   274 ; xlg2_cen = 0.2

   275 ; ylg2_cen = 760.

   277 ; mkres@gsMarkerColor = "brown"

   278 ; lnres@gsLineColor   = "brown"

   280 ; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)

   281 ; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \

   282 ;                                           (/ylg1_cen-60,ylg1_cen+60/),lnres)

   283 ; lg_cap_legend11  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \

   284 ;                                           (/ylg1_cen-60,ylg1_cen-60/),lnres)

   285 ; lg_cap_legend12  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \

   286 ;                                           (/ylg1_cen+60,ylg1_cen+60/),lnres)

   288 ; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)

   290 ; mkres@gsMarkerColor = "blue"

   291 ; lnres@gsLineColor   = "blue"

   293 ; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)

   294 ; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \

   295 ;                                           (/ylg2_cen-60,ylg2_cen+60/),lnres)

   296 ; lg_cap_legend21  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \

   297 ;                                           (/ylg2_cen-60,ylg2_cen-60/),lnres)

   298 ; lg_cap_legend22  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \

   299 ;                                           (/ylg2_cen+60,ylg2_cen+60/),lnres)

   300 ; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)

   302   draw(xy)

   303   frame(wks)

   305 end

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