load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl" ; useful for plotting calendar time as Y axis ; import data f=addfile("/mnt/geog/k/kg/kg312/Exercise/onset_x.nc","r") x=f->precip ; SJ onset ; Perform a 10-day running average x_sm_SJ=runave_n_Wrap(x,10,0,0) ; this implies that 4 points will be lost at the leftsmost side and 5 at the rightmost side ; zonal average x_zonal_SJ=dim_avg_n_Wrap(x_sm_SJ,2) ;Interpolating data to pass from a (0.5,1.5,2.5,...) grid to a (0,1,2,3,...) grid xi_SJ=x_zonal_SJ&lat ; get the latitude coordinates of xi_SJ(input coordinates) xo_SJ=ispan(0,21,1) ; define the new latitude coordinates (output coordinates) copy_VarMeta(xi_SJ,xo_SJ) xo_SJ=dim_pqsort( xo_SJ,-1); Sort in the same order as xi_SJ x_regrid_SJ=linint1_n_Wrap(xi_SJ,x_zonal_SJ,False,xo_SJ,0,1); False is to say not to process cyclically, 0 is set by default; 1 means the dimension of the latitude in x_zonal_SJ ; Change the type of lat cordinate of x_regrid_SJ to double Lat=int2dble(x_regrid_SJ&lat) delete(x_regrid_SJ&lat) x_regrid_SJ&lat=Lat xx=x_regrid_SJ(4:dimsizes(x_regrid_SJ&time)-6,1:dimsizes(x_regrid_SJ&lat)-1) ; remove the missing values from the time series (missing values come from the time running average and also from the linear interpolation, on coordinates that are extrapolated) delete(x_regrid_SJ) x_regrid_SJ=xx delete(xx) ; Reorder the time series from (time,lat) into (lat,time) x_reorder= x_regrid_SJ(lat|:,time|:) delete(x_regrid_SJ) x_regrid_SJ=x_reorder ; Changing time to yyyymmdd format time=f->time ;Extract the year TIM=cd_calendar(time,0) year=floattointeger(TIM(0,0)) print("----------------------------------") print(year) print("----------------------------------") ; get onset dates ;onset dates are expected to occur after 1 May may_1=cd_inv_calendar (year,05,01,00,00,00,"days since 1900-01-01 00:00:00",0 ) ;SJ_10 if (year .gt. 1996) then do i=ind(x_regrid_SJ&time .eq. may_1),dimsizes( x_regrid_SJ&time)-7,1 ; loop from 1 May to the end(-7 because in the following line j goes till i+6) j=ispan(i,i+6,1) cond10=x_regrid_SJ({10},j).gt. x_regrid_SJ({5},j) if (all(cond10)) then break end if end do SJ_10=x_regrid_SJ&time(i) do i=ind(x_regrid_SJ&time .eq. may_1),dimsizes( x_regrid_SJ&time)-7,1 j=ispan(i,i+6,1) cond15=x_regrid_SJ({15},j).gt. x_regrid_SJ({5},j) if (all(cond15)) then break end if end do SJ_15=x_regrid_SJ&time(i) else SJ_10="no data for "+ tostring(year) SJ_15="no data for "+ tostring(year) end if ; save the onset dates dat=asciiread("/mnt/geog/k/kg/kg312/Exercise/SJ_in.txt",(/20,3/),"integer") ;sj_10=doubletoint(SJ_10) ;sj_15=doubletoint(SJ_15) if (year .gt. 1996) then dat(year-1996,1)=cd_calendar(SJ_10,-2) dat(year-1996,2)=cd_calendar(SJ_15,-2) option=True option@fout="/mnt/geog/k/kg/kg312/Exercise/SJ_in.txt" write_matrix(dat,"2I10",option) end if ;Plotting plot_dir="/home/k/kg/kg312/Exercise/plot/" wks=gsn_open_wks("png", plot_dir+"contours"+tostring(year) ) ; workstation for contours wks2=gsn_open_wks("png", plot_dir+"shaded"+tostring(year) ) ; workstation for shades wks1=gsn_open_wks ("png",plot_dir+"lat-plot"+tostring(year)) ;workstation for plot by latitude res = True ; plot mods desired res@tiMainString = "monsoon-"+tostring(year) ; title ;res@tmYBLabelStride = 1 ; tick mark label stride res@tmYLOn=True ;turns on the left tick marks ;res@tiYAxisString = "time" ; y axis title res@lbOrientation="Vertical" res@tmXBLabelAngleF = 90 ; tilt the X-axis label res@lbLabelStrings = tostring(ispan(0,30,3)) ;contour level res@cnFillOn = False ; color on res@cnLinesOn = True ; turn on contour lines res@cnFillPalette =(/"white","grey","green","navyblue","blue","yellow","salmon","red","black"/) res@cnLevelSpacingF = 2 ; contour level spacing res@cnLevelSelectionMode = "ManualLevels" ; manual contour levels res@cnMinLevelValF = 1 ; min level res@cnMaxLevelValF = floattointeger(max(x_regrid_SJ))+1 ; max level ; Tick mark ;res@tmYBMode = "Manual" ;res@tmYBTickStartF=0 ; Y label minimum value ;res@tmYBTickStartF=0 ;res@tmYBTickEndF=20 ;res@tmYBTickSpacingF=5 ;time axis: put it in a calendar form restick = True restick@ttmFormat = "%d%C" ; ddMonYYYY restick@ttmAxis="XB" ; place X-axis label at the Bottom restick@ttmValues = (/(/year,1,15,0,0,0/),(/year,2,15,0,0,0/),(/year,3,15,0,0,0/),(/year,4,15,0,0,0/),(/year,5,15,0,0,0/),(/year,6,15,0,0,0/),(/year,7,15,0,0,0/),(/year,8,15,0,0,0/),(/year,9,15,0,0,0/),(/year,10,15,0,0,0/),(/year,11,15,0,0,0/),(/year,12,15,0,0,0/)/) restick@ttmMinorStride=10 time_axis_labels(x_regrid_SJ&time,res,restick) plot=gsn_csm_lat_time(wks, x_regrid_SJ, res) ;reproduce the same plot but with filled contours res2=res res2@cnFillOn = True ; color on time_axis_labels(x_regrid_SJ&time,res2,restick) plot2=gsn_csm_lat_time(wks2, x_regrid_SJ, res2) ; xy plot x_regrid_SJ_xy=x_regrid_SJ({5:15:5},:) res1 = True ; plot mods desired res1@tiMainString = "lat-plot "+tostring(year) res1@xyLineThicknesses = (/1.0,1.0,1.0/) res1@xyLineColors = (/"blue","green","red"/) ; change line color res1@tmXBLabelAngleF = 90 ; tilt the X-axis label ;res1@xyMonoMarkLineMode=True ;res1@xyMarkLineMode="Lines" ;res1@xyMonoMarker=True ;res1@xyMarker=0 res1@xyDashPattern = 0 ; only lines will be plotted res1@trYMinF = 0 res1@trYMaxF = 20 if (year .gt. 1996) then res1@gsnXRefLine=(/SJ_10,SJ_15/) ;vertical line for the onset dates res1@gsnXRefLineColor=(/"Black","grey"/) res1@gsnXRefLineThicknessF=2.0 time_axis_labels(x_regrid_SJ_xy&time,res1,restick) res1@trXMinF=min(x_regrid_SJ_xy&time) res1@trXMaxF=max(x_regrid_SJ_xy&time) plot1=gsn_csm_xy(wks1,x_regrid_SJ_xy&time,x_regrid_SJ_xy,res1) else ; no onset date time_axis_labels(x_regrid_SJ_xy&time,res1,restick) res1@trXMinF=min(x_regrid_SJ_xy&time) res1@trXMaxF=max(x_regrid_SJ_xy&time) plot1=gsn_csm_xy(wks1,x_regrid_SJ_xy&time,x_regrid_SJ_xy,res1)