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") 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