;************************************************* ; regress_1.ncl ;************************************************* load "/home/server/student/homes/naraligi/contributed/gsn_code.ncl" load "/home/server/student/homes/naraligi/contributed/gsn_csm.ncl" load "/home/server/student/homes/naraligi/contributed/shea_util.ncl" load "/home/server/student/homes/naraligi/contributed/contributed.ncl" begin nmb="0n155w" nmb1="0-155" cap=" SSUM" ;************************************************ ; create pointer to file and read in temperature ;************************************************ in = addfile("~/research/taodata/final/ncdf/intlt/lt"+nmb+".dat.new.nc.flt.nc","r") in1 = addfile("~/research/trmm/inttrmm/rain"+nmb1+".dat.nc.flt.nc","r") ;****************************************************** ;separate th data into seasons nsum=MJJASO ssum=NDJFMA ;************************************************************* time_lt=in->time time_rn=in1->time ymd_lt = ut_calendar(time_lt,0) ymd_rn = ut_calendar(time_rn,0) mon_lt = ymd_lt(:,1) mon_rn = ymd_rn(:,1) i = ind(mon_lt.eq.1 .or. mon_lt.eq.2 .or. \ mon_lt.eq.3 .or. mon_lt.eq.4 .or. mon_lt.eq.11 .or. mon_lt.eq.12) ; i = ind(mon_lt.eq.5 .or. mon_lt.eq.6 .or. \ ; mon_lt.eq.7 .or. mon_lt.eq.8 .or. mon_lt.eq.9 .or. mon_lt.eq.10) j = ind(mon_rn.eq.1 .or. mon_rn.eq.2 .or. \ mon_rn.eq.3 .or. mon_rn.eq.4 .or. mon_rn.eq.11 .or. mon_rn.eq.12) ; j = ind(mon_rn.eq.5 .or. mon_rn.eq.6 .or. \ ; mon_rn.eq.7 .or. mon_rn.eq.8 .or. mon_rn.eq.9 .or. mon_rn.eq.10) lt = in->lt(i) ; extract time series rn = in1->rain(j) rn = 28.935*rn rn@_FillValue = -9999.00 lt@_FillValue = -999.99 print(sizeof(lt)) print(sizeof(rn)) print(max(rn)) print(min(rn)) ;************************************************ ;note regline works on one dimensional arrays. ;************************************************ rc = regline(lt,rn) ; rcf = rc*(stddev(lt)/stddev(rn)) rcf = esccv(lt,rn,0)/(stddev(lt)*stddev(rn)) print(rc) ;************************************************ ; create an array to hold both the original data ; and the calculated regression line ;************************************************ data = new ( (/2,dimsizes(rn)/), typeof(rn)) data(0,:) = rn ; y = mx+b ; m is the slope: rc returned from regline ; b is the y intercept: rc@yave attribute of rc returned from regline data(1,:) = rc*(lt) + rc@yave ;************************************************ ; plotting parameters ;************************************************ wks = gsn_open_wks("ps",nmb+"_ssum") ; specifies a ps plot plot = new(1,graphic) res = True ; plot mods desired res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw res@xyMarkLineModes = (/"Markers","Lines"/) ; choose which have markers res@xyMarkers = 16 ; choose type of marker res@xyMarkerColor = "red" ; Marker color res@xyMarkerSizeF = 0.005 ; Marker size (default 0.01) res@xyDashPatterns = 1 ; solid line res@xyLineThicknesses = (/1,2/) ; set second line to 2 res@tiYAxisString = "Precipitation" plot = gsn_xy (wks,lt,data,res) ; create plot getvalues plot(0) "trXMaxF" : xmax "trYMaxF" : ymax end getvalues print(ymax) tres = True tres@txFontHeightF = 0.015 tres@txJust = "TopRight" gsn_text(wks,plot,rcf,xmax,ymax-50,tres) ;************************************************ resP = True ; modify the panel plot resP@txString = nmb+cap gsn_panel(wks,plot,(/1,1/),resP) end