<div dir="ltr"><div dir="ltr"><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:18.6667px">Hello</span></pre><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:14pt;font-family:""">I use ncl 6.5.0 working on EOF analysis, I got some issues when my masking my area of interest. My problem is how to mask out the analysis output according to my shapefile since the output displayed on places out of interest. I have attached the picture and script.</span></pre><pre><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:14pt;font-family:""">Please I need assistance,Thanks in Advance.</span></pre><pre><font color="#000000"><span style="font-size:14px">;==============================================================
;Calculate EOFs of the JFM Precipitation
; ==============================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
; ==============================================================
;Defined pancleters
; ==============================================================
yrStrt = 1961
yrLast = 2011
latS = -12.0
latN = 6.8
lonL = 30.0
lonR = 41.0
f= addfile("/mnt/c/Data/<a href="http://cru_ts3.23.1901.2014.pre.dat.nc">cru_ts3.23.1901.2014.pre.dat.nc</a>","r")
shp_n= addfile("Export3_Output.shp","r")
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
season="JFM"; choose June-July-Aug seasonal mean
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
neof = 2 ; number of EOFs
optEOF = True
optEOF@jopt = 0
optETS = False
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
lat2d = f->lat({latS:latN})
lon2d = f->lon({lonL:lonR})
;print(lon2d)
;exit
; Open the file: Read only the user specified period
;**************************************************************
TIME = f->time
opt=True
YYYY = cd_calendar(TIME,-1)/100 ; entire file
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
yrfrac = fspan(yrStrt, yrLast, yrLast-yrStrt+1)
precip1 = f->pre(iYYYY,{latS:latN},{lonL:lonR})
;precip2=calculate_monthly_values(precip1, "sum", 0, opt)
printVarSummary(precip1)
;exit
;******************************************
;Create the anomalies
;*****************************************
mean = clmMonTLL(precip1)
PRE_anom = calcMonAnomTLL(precip1,mean)
precip_JJAS_med=new((/3,51,24,22/),double)
do i=0,2
do j=0,50
precip_JJAS_med(i,j,:,:)=precip1(0+i+12*j,:,:)
end do
end do
precip_jjas = dim_avg_n_Wrap(precip_JJAS_med,0) ; jjas precipitation
precip_jjas!0 = "time"
precip_jjas!1 = "lat"
precip_jjas!2 = "lon"
precip_jjas@lat= lat2d
precip_jjas@lon= lon2d
;*****************************************************************
; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
; *****************************************************************
rad = 4.*atan(1.)/180.
clat = f->lat({latS:latN})
clat = sqrt( cos(rad*clat) )
;printVarSummary(clat)
;exit ; gw for gaussian grid
;*****************************************************************
; weight all observations
;*****************************************************************
WPRECIP = precip_jjas*conform(precip_jjas, clat, 1)
copy_VarMeta( precip_jjas, WPRECIP)
WPRECIP@long_name = "Wgt: "+WPRECIP@long_name
;printVarSummary(WPRECIP)
; exit
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
xw = WPRECIP({lat|latS:latN},{lon|lonL:lonR},time|:)
eof = eofunc_Wrap(xw, neof, optEOF)
eof = -1*eof
eof_ts = eofunc_ts_Wrap (xw, eof, optETS)
eof_ts = dim_standardize_n( eof_ts, 0, 1)
asciiwrite("GPCCeof2.txt",eof_ts)
;printMinMax(eof,0)
;printVarSummary(eof_ts)
;exit
;============================================================
; PLOTS
;============================================================
wks = gsn_open_wks("png","xxx_jjas_CR") ; send graphics to PNG file
gsn_define_colormap(wks,"BlueDarkred18")
;gsn_draw_colormap(wks)
;gsn_reverse_colormap(wks)
;gsn_draw_colormap(wks)
plot = new(neof,graphic) ; create graphic array
; only needed if paneling ; only needed if paneling ; create graphic array
; only needed if paneling
; EOF pattern
res = True
res@gsnMaximize = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnRightString = ""
res@gsnLeftString = ""
res@gsnAddCyclic = False
res@cnInfoLabelOn = False
res@cnFillOn = True
res@cnLinesOn = False
;res@cnLevelSelectionMode = "ExplicitLevels"
;res@cnLevels = fspan(-.15,.15,15)
res@cnFillDrawOrder = "PreDraw"
res@cnLineDrawOrder = "PreDraw"
res@cnLabelDrawOrder = "PreDraw"
res@cnLineLabelsOn = False
res@tmLabelAutoStride = 2
; res@tiMainString = "Correlation between Rainfall anomaly with DIMI"
;res@tiMainFontHeightF = 0.02
res@lbAutoManage = False
res@lbLabelStride = 1
res@lbLabelBarOn = False
res@lbBoxSizing = "UniformSizing" ; No single label bar
res@lbBoxSeparatorLinesOn = True
res@lbOrientation = "Vertical"
;res@pmLabelBarDisplayMode = "Always"
res@mpDataSetName = "Earth..4"
res@mpDataBaseVersion = "MediumRes"
res@mpOutlineOn = True
res@mpMinLatF = -12
res@mpMaxLatF = 6.8
res@mpMinLonF = 30
res@mpMaxLonF = 41
res@mpGridSpacingF = 0.55
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"tanzania"/)
res@mpLandFillColor = "White"
res@mpInlandWaterFillColor = "white"
res@mpOceanFillColor = "white"
res@mpOutlineBoundarySets = "National"
res@mpLandFillColor = "white"
res@pmTickMarkDisplayMode = "Always"
;res@lbLabelBarOn = False ; turn off individual lb's
res@mpOutlineOn = True ; Use outlines
;res@cnFillDrawOrder = "PreDraw"
;res@cnLineDrawOrder = "PreDraw"
;res@cnLabelDrawOrder = "PreDraw"
; rs@tmXBTickSpacingF =0.5
; exit()
;symMinMaxPlt(eof, 30, False, res)
res@tmXBLabelFontThicknessF = 1
res@tmXBLabelFontHeightF = 0.02
res@tmBorderThicknessF = 2
res@tmXBLabelStride = 2
res@tmYLLabelStride = 2
res@tmXBMinorOn = False
res@tmXTOn=False
res@tmYROn=False
;res@lbLabelStrings = sprintf ("%.1f",res@cnLevels)
res@mpGeophysicalLineThicknessF= 3
res@tmXBLabelFont = "times-bold"
res@tmYLLabelFont = "times-bold"
res@tiMainFont = "times-bold"
res@gsnStringFont = "times-bold"
res@cnLevelSelectionMode = "ManualLevels"
;res@cnLevels = fspan(-5,5,20)
res@cnMinLevelValF = -0.15 ; set min contour level
res@cnMaxLevelValF = 0.15 ; set max contour level
res@cnLevelSpacingF = 0.01
res@lbLabelStride= 2
res@gsnSpreadColorEnd = 2
res@gsnSpreadColorStart= -2
;res@lbAutoManage=False
;*************************************************************
;pannel plot resource
;*************************************************************
resP = True ; modify the panel plot
resP@gsnPanelLabelBar = True ; Add common label bar
;resP@gsnPanelFigureStrings =(/"a), b"
resP@gsnPanelFigureStringsFontHeightF=0.025
resP@gsnPanelFigureStringsBackgroundFillColor=-1
resP@gsnPanelFigureStringsPerimOn= False
resP@gsnPaperOrientation = "auto" ; automaticaly arange the plots in vertical or horizontaly
resP@gsnMaximize = True
resP@gsnStringFont = "times-bold"
;resP@lbLabelStride = 1
;resP@lbLabelStrings = sprintf ("%.1f",res@cnLevels)
resP@lbLabelFontHeightF = 0.02
resP@lbBoxSeparatorLinesOn = False
resP@pmLabelBarOrthogonalPosF = 0.01
resP@pmLabelBarWidthF = 0.8 ; default is shorter
resP@pmLabelBarHeightF = 0.1
resP@lbAutoManage = False
;resP@lbOrientation = "Vertical"
resP@lbLabelFontThicknessF = 2
resP@lbLabelFont = "times-bold"
;For Creating Panel
resp = True
resp@gsnFrame = False
resp@gsnPanelLabelBar =True
resp@gsnPanelBottom =0.05
resp@gsnPanelFigureStringsFontHeightF=0.015
;resp@amJust ="Topright"
resp@lbLabelStride = 2
resp@pmLabelBarOrthogonalPosF =0.01
resp@lbLabelFontHeightF = 0.015 ; set font height of Label Bar labels
resp@lbLabelFont = "times-bold"
resp@gsnPanelMainFont = "times-bold"
resp@txStringFont = "times-bold"
;*******************************************
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; first plot
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
do n=0,neof-1
res@gsnLeftString = "EOF "+(n+1)
res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n)=gsn_csm_contour_map(wks,eof(n,:,:),res)
end do
; ---Attach shapefile
lnres = True
lnres@gsLineColor = "black"
lnres@gsFillColor = "black"
lnres@gsLineThicknessF = 4
;---Setting lat/lon limits helps drawing go faster
lnres@minlat = -12
lnres@maxlat = -6.8
lnres@minlon = 30
lnres@maxlon = 41
;---Open shapefile and read lat/lon values.
dir_shp = "C:/cygwin/home/lenovo/"
shp_f= "Export3_Output.shp"
dum = gsn_add_shapefile_polylines(wks, plot, shp_f, lnres)
gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
; ;*******************************************
; second plot
;*******************************************
rts = True
rts@gsnDraw = False ; don't draw yet
rts@gsnFrame = False ; don't advance frame yet
rts@gsnScale = True ; force text scaling
rts@gsnYRefLine = 0 ; reference line
rts@gsnXYBarChart = True ; create bar chart
rts@gsnAboveYRefLineColor = "blue" ; above ref line fill red
rts@gsnBelowYRefLineColor = "red" ; below ref line fill blue
rts@tmYLMinorOn = False
rts@tmXBMinorOn = True
rts@tmXBLabelFontThicknessF = 2
rts@tmXBLabelFontHeightF = 0.02
rts@tmBorderThicknessF = 4
rts@tmYLLabelFontThicknessF = 2
rts@tmYLLabelFontHeightF = 0.02
rts@tmXTOn = False
rts@tmYROn = False
rts@vpHeightF = 0.40 ; Changes the aspect ratio
rts@vpWidthF = 0.85
rts@vpXF = 0.10 ; change start locations
rts@vpYF = 0.75 ; the plot
rts@tiYAxisString = "" ; y-axis label
rts@tmYLLabelFont = "times-bold"
rts@tiMainFont = "times-bold"
rts@gsnStringFont = "times-bold"
rts@tmXBLabelFont = "times-bold"
;*************************************************************
;pannel plot resource
;*************************************************************
rtsP = True ; modify the panel plot
rtsP@gsnMaximize = True ; large format
; rtsP@txString = "SLP: "+season+": "+yStrt+"-"+yLast
rtsP@gsnPanelFigureStrings = (/"a)","b)"/)
rtsP@gsnPanelFigureStringsFontHeightF = 0.025
rtsP@gsnPanelFigureStringsBackgroundFillColor= -1
rtsP@gsnPanelFigureStringsPerimOn = False
rtsP@gsnStringFont = "times-bold"
rtsP@lbLabelFont = "times-bold"
rtsP@gsnPanelYWhiteSpacePercent = 3 ;space between two plots
rtsP@gsnPanelXWhiteSpacePercent = 3
year = yrfrac
;print(year)
;exit
;create individual plots
do n=0,neof-1
rts@gsnLeftString = "EOF "+(n+1)
rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n) = gsn_csm_xy (wks,year,eof_ts(n,:),rts)
end do
gsn_panel(wks,plot,(/neof,1/),rtsP) ; now draw as one plot
poly=gsn_add_shapefile_polygons(wks,plot, shp_f, lnres)
;draw(plot)
frame(wks)
;exit
end<br></span></font></pre><font color="#000000"><span style="font-size:14px">
</span></font></pre><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:14pt;font-family:"""> </span></pre><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:14pt;font-family:"""> </span></pre><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:14pt;font-family:""">Best</span></pre><pre style="color:rgb(0,0,0);font-size:14px"><span style="font-size:14pt;font-family:""">Witty</span></pre></div></div>