[ncl-talk] (no subject)
witness massawe
witnessmassawe94 at gmail.com
Tue Mar 5 06:41:38 MST 2019
Hello
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.
Please I need assistance,Thanks in Advance.
;==============================================================
;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/cru_ts3.23.1901.2014.pre.dat.nc","r")
shp_n= addfile("Export3_Output.shp","r")
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
season="JFM"; choose June-July-Aug seasonal mean
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
neof = 2 ; number of EOFs
optEOF = True
optEOF at 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 at lat= lat2d
precip_jjas at 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 at long_name = "Wgt: "+WPRECIP at 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 at gsnMaximize = True
res at gsnDraw = False
res at gsnFrame = False
res at gsnRightString = ""
res at gsnLeftString = ""
res at gsnAddCyclic = False
res at cnInfoLabelOn = False
res at cnFillOn = True
res at cnLinesOn = False
;res at cnLevelSelectionMode = "ExplicitLevels"
;res at cnLevels = fspan(-.15,.15,15)
res at cnFillDrawOrder = "PreDraw"
res at cnLineDrawOrder = "PreDraw"
res at cnLabelDrawOrder = "PreDraw"
res at cnLineLabelsOn = False
res at tmLabelAutoStride = 2
; res at tiMainString = "Correlation between Rainfall anomaly
with DIMI"
;res at tiMainFontHeightF = 0.02
res at lbAutoManage = False
res at lbLabelStride = 1
res at lbLabelBarOn = False
res at lbBoxSizing = "UniformSizing" ; No
single label bar
res at lbBoxSeparatorLinesOn = True
res at lbOrientation = "Vertical"
;res at pmLabelBarDisplayMode = "Always"
res at mpDataSetName = "Earth..4"
res at mpDataBaseVersion = "MediumRes"
res at mpOutlineOn = True
res at mpMinLatF = -12
res at mpMaxLatF = 6.8
res at mpMinLonF = 30
res at mpMaxLonF = 41
res at mpGridSpacingF = 0.55
res at mpAreaMaskingOn = True
res at mpMaskAreaSpecifiers = (/"tanzania"/)
res at mpLandFillColor = "White"
res at mpInlandWaterFillColor = "white"
res at mpOceanFillColor = "white"
res at mpOutlineBoundarySets = "National"
res at mpLandFillColor = "white"
res at pmTickMarkDisplayMode = "Always"
;res at lbLabelBarOn = False ; turn off individual lb's
res at mpOutlineOn = True ; Use outlines
;res at cnFillDrawOrder = "PreDraw"
;res at cnLineDrawOrder = "PreDraw"
;res at cnLabelDrawOrder = "PreDraw"
; rs at tmXBTickSpacingF =0.5
; exit()
;symMinMaxPlt(eof, 30, False, res)
res at tmXBLabelFontThicknessF = 1
res at tmXBLabelFontHeightF = 0.02
res at tmBorderThicknessF = 2
res at tmXBLabelStride = 2
res at tmYLLabelStride = 2
res at tmXBMinorOn = False
res at tmXTOn=False
res at tmYROn=False
;res at lbLabelStrings = sprintf ("%.1f",res at cnLevels)
res at mpGeophysicalLineThicknessF= 3
res at tmXBLabelFont = "times-bold"
res at tmYLLabelFont = "times-bold"
res at tiMainFont = "times-bold"
res at gsnStringFont = "times-bold"
res at cnLevelSelectionMode = "ManualLevels"
;res at cnLevels = fspan(-5,5,20)
res at cnMinLevelValF = -0.15 ; set min contour level
res at cnMaxLevelValF = 0.15 ; set max contour level
res at cnLevelSpacingF = 0.01
res at lbLabelStride= 2
res at gsnSpreadColorEnd = 2
res at gsnSpreadColorStart= -2
;res at lbAutoManage=False
;*************************************************************
;pannel plot resource
;*************************************************************
resP = True ; modify the panel plot
resP at gsnPanelLabelBar = True ; Add common label bar
;resP at gsnPanelFigureStrings =(/"a), b"
resP at gsnPanelFigureStringsFontHeightF=0.025
resP at gsnPanelFigureStringsBackgroundFillColor=-1
resP at gsnPanelFigureStringsPerimOn= False
resP at gsnPaperOrientation = "auto" ; automaticaly arange the
plots in vertical or horizontaly
resP at gsnMaximize = True
resP at gsnStringFont = "times-bold"
;resP at lbLabelStride = 1
;resP at lbLabelStrings = sprintf ("%.1f",res at cnLevels)
resP at lbLabelFontHeightF = 0.02
resP at lbBoxSeparatorLinesOn = False
resP at pmLabelBarOrthogonalPosF = 0.01
resP at pmLabelBarWidthF = 0.8 ; default is shorter
resP at pmLabelBarHeightF = 0.1
resP at lbAutoManage = False
;resP at lbOrientation = "Vertical"
resP at lbLabelFontThicknessF = 2
resP at lbLabelFont = "times-bold"
;For Creating Panel
resp = True
resp at gsnFrame = False
resp at gsnPanelLabelBar =True
resp at gsnPanelBottom =0.05
resp at gsnPanelFigureStringsFontHeightF=0.015
;resp at amJust ="Topright"
resp at lbLabelStride = 2
resp at pmLabelBarOrthogonalPosF =0.01
resp at lbLabelFontHeightF = 0.015 ; set font height of Label Bar labels
resp at lbLabelFont = "times-bold"
resp at gsnPanelMainFont = "times-bold"
resp at txStringFont = "times-bold"
;*******************************************
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; first plot
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
do n=0,neof-1
res at gsnLeftString = "EOF "+(n+1)
res at gsnRightString = sprintf("%5.1f", eof at pcvar(n)) +"%"
plot(n)=gsn_csm_contour_map(wks,eof(n,:,:),res)
end do
; ---Attach shapefile
lnres = True
lnres at gsLineColor = "black"
lnres at gsFillColor = "black"
lnres at gsLineThicknessF = 4
;---Setting lat/lon limits helps drawing go faster
lnres at minlat = -12
lnres at maxlat = -6.8
lnres at minlon = 30
lnres at 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 at gsnDraw = False ; don't draw yet
rts at gsnFrame = False ; don't advance frame yet
rts at gsnScale = True ; force text scaling
rts at gsnYRefLine = 0 ; reference line
rts at gsnXYBarChart = True ; create bar chart
rts at gsnAboveYRefLineColor = "blue" ; above ref line fill red
rts at gsnBelowYRefLineColor = "red" ; below ref line fill blue
rts at tmYLMinorOn = False
rts at tmXBMinorOn = True
rts at tmXBLabelFontThicknessF = 2
rts at tmXBLabelFontHeightF = 0.02
rts at tmBorderThicknessF = 4
rts at tmYLLabelFontThicknessF = 2
rts at tmYLLabelFontHeightF = 0.02
rts at tmXTOn = False
rts at tmYROn = False
rts at vpHeightF = 0.40 ; Changes the aspect ratio
rts at vpWidthF = 0.85
rts at vpXF = 0.10 ; change start locations
rts at vpYF = 0.75 ; the plot
rts at tiYAxisString = "" ; y-axis label
rts at tmYLLabelFont = "times-bold"
rts at tiMainFont = "times-bold"
rts at gsnStringFont = "times-bold"
rts at tmXBLabelFont = "times-bold"
;*************************************************************
;pannel plot resource
;*************************************************************
rtsP = True ; modify the panel plot
rtsP at gsnMaximize = True ; large format
; rtsP at txString = "SLP: "+season+": "+yStrt+"-"+yLast
rtsP at gsnPanelFigureStrings = (/"a)","b)"/)
rtsP at gsnPanelFigureStringsFontHeightF = 0.025
rtsP at gsnPanelFigureStringsBackgroundFillColor= -1
rtsP at gsnPanelFigureStringsPerimOn = False
rtsP at gsnStringFont = "times-bold"
rtsP at lbLabelFont = "times-bold"
rtsP at gsnPanelYWhiteSpacePercent = 3 ;space between two plots
rtsP at gsnPanelXWhiteSpacePercent = 3
year = yrfrac
;print(year)
;exit
;create individual plots
do n=0,neof-1
rts at gsnLeftString = "EOF "+(n+1)
rts at gsnRightString = sprintf("%5.1f", eof at 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
Best
Witty
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190305/121b7981/attachment.html>
More information about the ncl-talk
mailing list