Re: superimposing Gaussian PDF onto histogram

Hi Kerrie,

In the histogram, the X axs is nothing more than X values that go
from roughly 0 to 1.0. You are trying to draw polylines and markers on
top of it whose X values go from 0 to 30000. These two axes are not in
the same range, so that's why your gsn_add_polyxxxx calls are not
working.

To get the gsn_add_polyxxx calls to work, you need to map the 0-30000
values into the actual range of the histogram plot. The histogram X
values has three sets of "x" values: the beginning of each bar
(plot@BeginBarLocs), the mid point of each bar (plot1@MidBarLocs), and
then end of each bar (plot1@EndBarLocs).

So, if we decided to map the range of your markers to the midpoint of
the histogram, then you would use something like:

nbars = dimsizes(plot1@MidBarLocs)
nx = dimsizes(x)
xrange_hist = plot1@MidBarLocs(nbars-1) - plot1@MidBarLocs(0)
xrange_xy = x(nx-1) - x(0)
xnorm = ((x-x(0))/xrange_xy) * xrange_hist + plot1@MidBarLocs(0)

do i=0,dimsizes(gausDJFL)-2
   xvalues=(/xnorm(i),xnorm(i+1)/)
   yvalues=(/gausDJFL(i),gausDJFL(i+1)/)
   dum2(i)=gsn_add_polyline(wks,plot1,xvalues,yvalues,gres)
   dum3(i)=gsn_add_polymarker(wks,plot1,xnorm,gausDJFL,mkres)
end do

--Mary

On Sep 12, 2011, at 6:02 PM, Kerrie Geil wrote:

> Hi,
>
> I'm trying to superimpose gausian and gamma PDFs on top of a
> histogram to show which is a better fit for a course I'm taking.
> I'm plotting the histogram and then using gsn_add_polyline to add
> the PDFs, but can't get the PDFs to scale correctly over the
> histogram. If I plot the PDFs alone, they look fine. Is there
> something with gsn_histogram that will not allow this?
>
> My code is attached along with the small dataset I'm using.
>
> Thanks for the help!
>
> Kerrie
>
> ;###############################################
> load "$NCARG_NCARG/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_NCARG/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>
> begin
>
> ;##### READ DATA ####
> lees_ferry_data=asciiread("/Users/geil/Documents/stats_class/HW1/
> lees_ferry.txt",-1,"string") ;#### EDIT LOCATION HERE #####
>
> nrows_lees=dimsizes(lees_ferry_data)
> time=new((/nrows_lees/),string)
>
> ;##### CREATE A TIME COORDINATE #####
> do i=0,nrows_lees-1
> if(toint(str_get_field(lees_ferry_data(i),2," ")).lt.10)then
> time(i)=str_get_field(lees_ferry_data(i),1," ")
> +"0"+str_get_field(lees_ferry_data(i),2," ")
> else
> time(i)=str_get_field(lees_ferry_data(i),1," ")
> +str_get_field(lees_ferry_data(i),2," ")
> end if
> end do
>
> ;##### SEPARATE OUT THE STREAM FLOW DATA #####
> flow_lees=tofloat(str_get_field(lees_ferry_data,3," "))
> flow_lees!0="time"
> flow_lees&time=time
>
> ;##### SEPARATE WINTER DATA #####
> DJFindicesL=ind( (toint(str_get_cols(flow_lees&time,4,5)) .le.
> 2) .or. (toint(str_get_cols(flow_lees&time,4,5)).eq.12) )
> DJF_L=flow_lees(DJFindicesL)
>
> ;##### CALCULATE GAUSSIAN #####
> mean_DJFL=sum(DJF_L)/dimsizes(.not.(ismissing(DJF_L)))
> stddev_DJFL=(sum((DJF_L-mean_DJFL)^2)/num(.not.
> (ismissing(DJF_L))))^0.5
>
> pi=3.14159
> qsort(DJF_L)
>
> x1=ispan(0,30000,100)
> gausDJFL=(1./(stddev_DJFL*((2*pi)^0.5)))*exp(-(((x1-mean_DJFL)^2)/
> (2*stddev_DJFL^2)))
>
>
> ;#################################################################################
> ;##### PLOTTING SECTION
> ##########################################################
> ;#################################################################################
> res=True
> res@gsnMaximize=True
> res@gsnHistogramSelectNiceIntervals=False
> res@tmXBLabelAngleF=270.
> res@gsFillColor=(/"transparent"/)
> res@gsEdgeColor="black"
> res@tmXBLabelFontHeightF=0.012
> res@tmXTBorderOn=False
> res@tmYRBorderOn=False
> ;res@tmYRMinorOn=False
> res@tmYROn=False
> res@tmXTOn=False
> res@tmXBLabelJust="CenterCenter"
> res@tmXBMajorLengthF=0.01
> res@tmYLMinorLengthF=0.01
> res@tiMainFontHeightF=0.012
> res@tiXAxisFontHeightF=0.012
> res@tiYAxisFontHeightF=0.012
>
>
> wks=gsn_open_wks("x11","DJF_L")
>
> res@gsnDraw=False
> res@gsnFrame=False
>
> txres=True
> txres@txJust="BottomCenter"
> txres@txFontHeightF=0.013
>
> plotmin=min(DJF_L)
> plotmax=max(DJF_L)
> h=1633.5 ;##### BIN WIDTH #####
> res@gsnHistogramNumberOfBins=round((plotmax-plotmin)/h,3)
> plot1=gsn_histogram(wks,decimalPlaces(DJF_L,1,True),res)
> labels=""+plot1@NumInBins
> dum1=gsn_add_text(wks,plot1,labels,plot1@MidBarLocs,plot1@NumInBins+0.5
> ,txres)
>
> gres=True
> gres@gsLineColor="blue"
> mkres=True
> mkres@gsMarkerIndex=17
> mkres@gsMarkerSizeF=0.035
> mkres@gsMarkerColor="blue"
>
> x=fspan(0,30000,301)
> dum2=new(dimsizes(gausDJFL)-1,graphic)
> dum3=new(dimsizes(gausDJFL)-1,graphic)
> gausDJFL=gausDJFL*102910.5 ;##### 102910.5 is No.Data Values * Bin
> Width=scaling factor #####
>
>
> do i=0,dimsizes(gausDJFL)-2
> xvalues=(/x(i),x(i+1)/)
> yvalues=(/gausDJFL(i),gausDJFL(i+1)/)
> dum2(i)=gsn_add_polyline(wks,plot1,xvalues,yvalues,gres)
> dum3(i)=gsn_add_polymarker(wks,plot1,x,gausDJFL,mkres)
> end do
>
> draw(plot1)
> frame(wks)
>
> plot=gsn_xy(wks,x,gausDJFL,res)
> draw(plot)
> frame(wks)
>
> end
> ;##########################################################
>
> --
> Kerrie Geil
>
> Master's Student
> Department of Atmospheric Sciences
> University of Arizona
> PAS Building Rm 526
> 1118 E 4th Street
> Tucson, AZ 85721
>
>
> <
> lees_ferry
> .txt><ncl_talk.ncl>_______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Sep 12 21:58:57 2011