<div>Hi:</div><div>&nbsp; All,I am trying to mask data with shapefile according &nbsp;to the example shapefiles_14_mask.ncl on the NCL website.</div><div>B<span style="line-height: 1.5;">ut the masked plot did not match the original plot,Ican not figure out what's wrong with it.</span></div><div>Left plot is the masked plot,right is the original one.</div><div><img src="cid:D0881338@B8CA7277.E358A457.jpg" filesize="109260" modifysize="64%" diffpixels="12px" scalingmode="zoom" style="width: 487px; height: 191px;"></div><div>This is my script:</div><div><div>; &nbsp; Example script to produce plots for a WRF real-data run,</div><div>; &nbsp; with the ARW coordinate dynamics option.</div><div><br></div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" </div><div>load"./shapefile_mask_data.ncl"</div><div>begin</div><div>;</div><div>; The WRF ARW input file. &nbsp;</div><div>; This needs to have a ".nc" appended, so just do it.</div><div>&nbsp; a = addfile("/public/home/huanglei/data/20160724/oldwrf_cu5/wrfout_d03_2016-07-24_00:00:00"+".nc","r")</div><div><br></div><div>; We generate plots, but what kind do we prefer?</div><div>&nbsp;; type = "x11"</div><div>&nbsp; type = "png"</div><div>; type = "ps"</div><div>; type = "ncgm"</div><div>&nbsp; wks = gsn_open_wks(type,"mask_plt_Precip_oldwrf5km_zoom_cu5_from24")</div><div>&nbsp; plot = new(3,graphic)</div><div><br></div><div>; Set some basic resources</div><div>&nbsp; res = True</div><div>&nbsp; res@MainTitle = "REAL-TIME WRF"</div><div><br></div><div>&nbsp; mpres &nbsp;= True &nbsp;; Map resources</div><div>&nbsp; mpres@mpOutlineOn = False &nbsp;; Turn off map outlines</div><div>&nbsp; mpres@mpFillOn &nbsp; &nbsp;= False &nbsp;; Turn off map fill</div><div>&nbsp; mpres@mpGridAndLimbOn = True</div><div>&nbsp;;res@mpProjection &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;= "Lambert"</div><div>&nbsp; pltres = True ; Plot resources</div><div>&nbsp; pltres@PanelPlot &nbsp;= True &nbsp; ; Tells wrf_map_overlays not to remove overlays</div><div><br></div><div><br></div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div><br></div><div>; What times and how many time steps are in the data set?</div><div>&nbsp; FirstTime = True</div><div>&nbsp; times = wrf_user_getvar(a,"times",-1) &nbsp;; get all times in the file</div><div>&nbsp; ntimes = dimsizes(times) &nbsp; &nbsp; &nbsp; &nbsp; ; number of times in the file</div><div>&nbsp;; print(times)</div><div>&nbsp;; exit</div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div>&nbsp; &nbsp;it_start = 10</div><div>&nbsp; &nbsp;it_end = 17</div><div>&nbsp; &nbsp;; print("Working on time: " + times(it) )</div><div>&nbsp; &nbsp;; res@TimeLabel = times(it) &nbsp; ; Set Valid time to use on plots</div><div><br></div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div>; First get the variables we will need &nbsp; &nbsp; &nbsp; &nbsp;</div><div><br></div><div><br></div><div>&nbsp; ; Get non-convective, convective and total precipitation of 5km &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;</div><div>&nbsp; &nbsp; rain_exp = wrf_user_getvar(a,"RAINNC",it_end)</div><div>&nbsp; &nbsp; rain_con = wrf_user_getvar(a,"RAINC",it_end)</div><div>&nbsp; &nbsp; rain_tot = rain_exp + rain_con</div><div>&nbsp; &nbsp; rain_tot@description = "Total Precipitation"</div><div><br></div><div><br></div><div><span class="Apple-tab-span" style="white-space:pre">        </span>;calculate the precipitation </div><div>&nbsp; &nbsp; &nbsp; &nbsp; rain_exp_save = wrf_user_getvar(a,"RAINNC",it_start)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rain_con_save = wrf_user_getvar(a,"RAINC",it_start)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rain_tot_save = rain_exp_save + rain_con_save</div><div><br></div><div><br></div><div>&nbsp; &nbsp; &nbsp; &nbsp; times_sav = times(it_start)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rain_tot_tend = rain_tot - rain_tot_save</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rainc_tend = rain_con - rain_con_save &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;; CUMULUS PRECIPITATION</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rainnc_tend= rain_exp - rain_exp_save &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;; SCALE PRECIPITATION</div><div><br></div><div>&nbsp; &nbsp; &nbsp; &nbsp; rain_tot_tend@description = "Precipitation of 5km(old wrf)"</div><div><div>&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;rain_tot_tend@lat2d = a-&gt;XLAT(1,:,:)</div><div>&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;rain_tot_tend@lon2d = a-&gt;XLONG(1,:,:)</div></div><div><br></div><div>&nbsp; &nbsp; &nbsp; &nbsp; rainc_tend@description = "RAINC of 5km(old wrf)"</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rainnc_tend@description = "RAINNC of 5km(old wrf)"</div><div><br></div><div>;;;;;;;;;;;;;;;;mask data with shapefile;;;;;;;;;;;;;;;;;;;;;;;;</div><div>shp_filename = "/public/home/huanglei/map/xian.shp"</div><div>rain_tot_mask = shapefile_mask_data(rain_tot_tend,shp_filename,True)</div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div><br></div><div><br></div><div>&nbsp; &nbsp; &nbsp; ; Plotting options for Precipitation</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r = res &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@UnitLabel &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;= "mm"</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@cnLevelSelectionMode = "ExplicitLevels"</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@cnLevels &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 12.8, 25.6, 51.2/)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@cnFillColors &nbsp; &nbsp; &nbsp; &nbsp; = (/"White","DarkOliveGreen1", \</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; "DarkOliveGreen3","Chartreuse", \</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; "Chartreuse3","Green","ForestGreen", \</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; "Yellow","Orange","Red","Violet"/)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@cnInfoLabelOn &nbsp; &nbsp; &nbsp; &nbsp;= False</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@cnConstFLabelOn &nbsp; &nbsp; &nbsp;= False</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@cnFillOn &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; = True</div><div>&nbsp; &nbsp; &nbsp; ; &nbsp;opts_r@vpHeightF &nbsp; = 0.1</div><div>&nbsp; &nbsp; &nbsp; ; &nbsp;opts_r@vpWidthF &nbsp; &nbsp;= 0.9 </div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@gsnDraw &nbsp; &nbsp; &nbsp;= &nbsp;False &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@gsnFrame &nbsp; &nbsp; = &nbsp;False</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@lbLabelBarOn &nbsp;= False</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@Footer &nbsp; &nbsp; = False</div><div>&nbsp; &nbsp; &nbsp; &nbsp; opts_r@NoHeaderFooter =True </div><div>&nbsp; &nbsp; &nbsp; &nbsp; ;;;;;;set zoom ;;;;;</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;tes = True</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;tes@returnInt = False</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;loc1=wrf_user_ll_to_ij(a,107.50,33.5,tes)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;print("X/Y location is: "+ loc1)</div><div>&nbsp; &nbsp;</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;loc2=wrf_user_ll_to_ij(a,110.0,35.2,tes)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;print("X/Y location is: "+ loc2)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ; &nbsp;exit</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;x_start = 35</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;x_end &nbsp; = 81</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;y_start = 47</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;y_end &nbsp; = 84</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mpres1 = True</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mpres1@ZoomIn = True</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mpres1@Xstart = x_start</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mpres1@Ystart = y_start</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mpres1@Xend &nbsp; = x_end</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; mpres1@Yend &nbsp; = y_end</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rain_tot_zoom = rain_tot_tend(y_start:y_end,x_start:x_end)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rainc_zoom = rainc_tend(y_start:y_end,x_start:x_end)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; rainnc_zoom = rainnc_tend(y_start:y_end,x_start:x_end)</div><div>&nbsp; &nbsp; &nbsp; ; Precipitation Tendencies</div><div>&nbsp; &nbsp; &nbsp;; &nbsp; opts_r@SubFieldTitle = "from " + times(it_start) + " to " + times(it_end)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; contour_tend = wrf_contour(a,wks, rain_tot_mask,opts_r) ; total (color)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; contour_rainc_tend = wrf_contour(a,wks, rainc_zoom,opts_r) ; total cumulus precipitation (color)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; contour_rainnc_tend = wrf_contour(a,wks, rainnc_zoom,opts_r) ; total scale precipitation(color)</div><div><br></div><div>&nbsp; &nbsp; &nbsp; &nbsp; delete(opts_r)</div><div><br></div><div>&nbsp; &nbsp; &nbsp; ; MAKE PLOTS &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; </div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;plot(0) = wrf_map_overlays(a,wks,contour_tend,pltres,mpres1)</div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;plot(1) = wrf_map_overlays(a,wks,contour_rainc_tend,pltres,mpres1) </div><div>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;plot(2) = wrf_map_overlays(a,wks,contour_rainnc_tend,pltres,mpres1)</div><div><br></div><div><br></div><div>&nbsp; </div><div>;&gt;============================================================&lt;</div><div>; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;add China map</div><div>;&gt;------------------------------------------------------------&lt;</div><div><br></div><div>&nbsp; &nbsp; &nbsp;</div><div>&nbsp; shp_name1 &nbsp; &nbsp;= "/public/home/huanglei/map/xian.shp"</div><div><br></div><div>&nbsp; lnres &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;= True</div><div>&nbsp; lnres@gsLineColor &nbsp; &nbsp; &nbsp;= "gray25"</div><div>&nbsp; lnres@gsLineThicknessF = 0.5 &nbsp; </div><div>&nbsp;id = new(3,graphic)</div><div>&nbsp;id(0) = gsn_add_shapefile_polylines(wks,plot(0),shp_name1,lnres)</div><div>&nbsp;id(1) = gsn_add_shapefile_polylines(wks,plot(1),shp_name1,lnres)</div><div>&nbsp;id(2) = gsn_add_shapefile_polylines(wks,plot(2),shp_name1,lnres)</div><div>&nbsp;</div><div><br></div><div><br></div><div>&nbsp;shp_name2 &nbsp; &nbsp;= "/public/home/huanglei/map/shaanxi_city_l.shp"</div><div><br></div><div>&nbsp; prres=True</div><div>&nbsp; prres@gsLineThicknessF = 1.0</div><div>&nbsp; prres@gsLineColor = "black"</div><div>&nbsp;id2 = new(3,graphic)</div><div>&nbsp;id2(0) = gsn_add_shapefile_polylines(wks,plot(0),shp_name2,prres)</div><div>&nbsp;id2(1) = gsn_add_shapefile_polylines(wks,plot(1),shp_name2,prres)</div><div>&nbsp;id2(2) = gsn_add_shapefile_polylines(wks,plot(2),shp_name2,prres)</div><div>&nbsp;</div><div><br></div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div><br></div><div>;;;;;;;;;;;;;;;;;;creat panel;;;;;;;;;;;;;;;</div><div>&nbsp; resP &nbsp; &nbsp; = True</div><div>&nbsp; resP@gsnMaximize = True</div><div>&nbsp; resP@gsnFrame &nbsp;=False</div><div>&nbsp; resP@gsnPanelLabelBar = True</div><div>&nbsp; resP@gsnPanelBottom = 0.05</div><div>&nbsp; resP@gsnPanelTop &nbsp; &nbsp; &nbsp; = 0.95 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ; Make sure not too close to </div><div>&nbsp; resP@gsnPanelBottom &nbsp; &nbsp;= 0.001 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; ; edge, so it maximizes better.</div><div>&nbsp; resP@gsnPanelLabelBar &nbsp;= True</div><div>&nbsp; resP@lbLabelFontHeightF = &nbsp;0.006</div><div>&nbsp; resP@gsnPanelFigureStrings= (/"Total precipitation","RAINC","RAINNC"/) ; add strings to panel</div><div>&nbsp; resP@amJust &nbsp; = "TopRight" </div><div>&nbsp; resP@gsnPanelFigureStringsFontHeightF = 0.005</div><div>&nbsp; resP@lbLabelAutoStride = True</div><div>&nbsp; resP@gsnPaperOrientation = "landscape"</div><div>&nbsp; gsn_panel(wks,plot,(/1,3/),resP)</div><div><br></div><div><br></div><div>&nbsp;;psres = True &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; </div><div>; maximize_output(wks,psres) &nbsp;; calls draw and frame for you </div><div><br></div><div><br></div><div>; draw(plot) &nbsp; &nbsp; &nbsp; ; This will draw the map and the shapefile outlines.</div><div>&nbsp; frame(wks)</div><div><br></div><div>end</div></div><div><br></div><div><span style="font-family: 'lucida Grande', Verdana, 'Microsoft YaHei'; line-height: 23.8px;">&nbsp;How can I slove the problem?</span></div><div></div>