[ncl-talk] data did not match with original plot after mask data with shapefile
grace
313695096 at qq.com
Fri Aug 5 03:14:11 MDT 2016
Hi:
All,I am trying to mask data with shapefile according to the example shapefiles_14_mask.ncl on the NCL website.
But the masked plot did not match the original plot,Ican not figure out what's wrong with it.
Left plot is the masked plot,right is the original one.
This is my script:
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load"./shapefile_mask_data.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
a = addfile("/public/home/huanglei/data/20160724/oldwrf_cu5/wrfout_d03_2016-07-24_00:00:00"+".nc","r")
; We generate plots, but what kind do we prefer?
; type = "x11"
type = "png"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"mask_plt_Precip_oldwrf5km_zoom_cu5_from24")
plot = new(3,graphic)
; Set some basic resources
res = True
res at MainTitle = "REAL-TIME WRF"
mpres = True ; Map resources
mpres at mpOutlineOn = False ; Turn off map outlines
mpres at mpFillOn = False ; Turn off map fill
mpres at mpGridAndLimbOn = True
;res at mpProjection = "Lambert"
pltres = True ; Plot resources
pltres at PanelPlot = True ; Tells wrf_map_overlays not to remove overlays
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
FirstTime = True
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
; print(times)
; exit
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
it_start = 10
it_end = 17
; print("Working on time: " + times(it) )
; res at TimeLabel = times(it) ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
; Get non-convective, convective and total precipitation of 5km
rain_exp = wrf_user_getvar(a,"RAINNC",it_end)
rain_con = wrf_user_getvar(a,"RAINC",it_end)
rain_tot = rain_exp + rain_con
rain_tot at description = "Total Precipitation"
;calculate the precipitation
rain_exp_save = wrf_user_getvar(a,"RAINNC",it_start)
rain_con_save = wrf_user_getvar(a,"RAINC",it_start)
rain_tot_save = rain_exp_save + rain_con_save
times_sav = times(it_start)
rain_tot_tend = rain_tot - rain_tot_save
rainc_tend = rain_con - rain_con_save ; CUMULUS PRECIPITATION
rainnc_tend= rain_exp - rain_exp_save ; SCALE PRECIPITATION
rain_tot_tend at description = "Precipitation of 5km(old wrf)"
rain_tot_tend at lat2d = a->XLAT(1,:,:)
rain_tot_tend at lon2d = a->XLONG(1,:,:)
rainc_tend at description = "RAINC of 5km(old wrf)"
rainnc_tend at description = "RAINNC of 5km(old wrf)"
;;;;;;;;;;;;;;;;mask data with shapefile;;;;;;;;;;;;;;;;;;;;;;;;
shp_filename = "/public/home/huanglei/map/xian.shp"
rain_tot_mask = shapefile_mask_data(rain_tot_tend,shp_filename,True)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Plotting options for Precipitation
opts_r = res
opts_r at UnitLabel = "mm"
opts_r at cnLevelSelectionMode = "ExplicitLevels"
opts_r at cnLevels = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \
12.8, 25.6, 51.2/)
opts_r at cnFillColors = (/"White","DarkOliveGreen1", \
"DarkOliveGreen3","Chartreuse", \
"Chartreuse3","Green","ForestGreen", \
"Yellow","Orange","Red","Violet"/)
opts_r at cnInfoLabelOn = False
opts_r at cnConstFLabelOn = False
opts_r at cnFillOn = True
; opts_r at vpHeightF = 0.1
; opts_r at vpWidthF = 0.9
opts_r at gsnDraw = False
opts_r at gsnFrame = False
opts_r at lbLabelBarOn = False
opts_r at Footer = False
opts_r at NoHeaderFooter =True
;;;;;;set zoom ;;;;;
tes = True
tes at returnInt = False
loc1=wrf_user_ll_to_ij(a,107.50,33.5,tes)
print("X/Y location is: "+ loc1)
loc2=wrf_user_ll_to_ij(a,110.0,35.2,tes)
print("X/Y location is: "+ loc2)
; exit
x_start = 35
x_end = 81
y_start = 47
y_end = 84
mpres1 = True
mpres1 at ZoomIn = True
mpres1 at Xstart = x_start
mpres1 at Ystart = y_start
mpres1 at Xend = x_end
mpres1 at Yend = y_end
rain_tot_zoom = rain_tot_tend(y_start:y_end,x_start:x_end)
rainc_zoom = rainc_tend(y_start:y_end,x_start:x_end)
rainnc_zoom = rainnc_tend(y_start:y_end,x_start:x_end)
; Precipitation Tendencies
; opts_r at SubFieldTitle = "from " + times(it_start) + " to " + times(it_end)
contour_tend = wrf_contour(a,wks, rain_tot_mask,opts_r) ; total (color)
contour_rainc_tend = wrf_contour(a,wks, rainc_zoom,opts_r) ; total cumulus precipitation (color)
contour_rainnc_tend = wrf_contour(a,wks, rainnc_zoom,opts_r) ; total scale precipitation(color)
delete(opts_r)
; MAKE PLOTS
plot(0) = wrf_map_overlays(a,wks,contour_tend,pltres,mpres1)
plot(1) = wrf_map_overlays(a,wks,contour_rainc_tend,pltres,mpres1)
plot(2) = wrf_map_overlays(a,wks,contour_rainnc_tend,pltres,mpres1)
;>============================================================<
; add China map
;>------------------------------------------------------------<
shp_name1 = "/public/home/huanglei/map/xian.shp"
lnres = True
lnres at gsLineColor = "gray25"
lnres at gsLineThicknessF = 0.5
id = new(3,graphic)
id(0) = gsn_add_shapefile_polylines(wks,plot(0),shp_name1,lnres)
id(1) = gsn_add_shapefile_polylines(wks,plot(1),shp_name1,lnres)
id(2) = gsn_add_shapefile_polylines(wks,plot(2),shp_name1,lnres)
shp_name2 = "/public/home/huanglei/map/shaanxi_city_l.shp"
prres=True
prres at gsLineThicknessF = 1.0
prres at gsLineColor = "black"
id2 = new(3,graphic)
id2(0) = gsn_add_shapefile_polylines(wks,plot(0),shp_name2,prres)
id2(1) = gsn_add_shapefile_polylines(wks,plot(1),shp_name2,prres)
id2(2) = gsn_add_shapefile_polylines(wks,plot(2),shp_name2,prres)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;creat panel;;;;;;;;;;;;;;;
resP = True
resP at gsnMaximize = True
resP at gsnFrame =False
resP at gsnPanelLabelBar = True
resP at gsnPanelBottom = 0.05
resP at gsnPanelTop = 0.95 ; Make sure not too close to
resP at gsnPanelBottom = 0.001 ; edge, so it maximizes better.
resP at gsnPanelLabelBar = True
resP at lbLabelFontHeightF = 0.006
resP at gsnPanelFigureStrings= (/"Total precipitation","RAINC","RAINNC"/) ; add strings to panel
resP at amJust = "TopRight"
resP at gsnPanelFigureStringsFontHeightF = 0.005
resP at lbLabelAutoStride = True
resP at gsnPaperOrientation = "landscape"
gsn_panel(wks,plot,(/1,3/),resP)
;psres = True
; maximize_output(wks,psres) ; calls draw and frame for you
; draw(plot) ; This will draw the map and the shapefile outlines.
frame(wks)
end
How can I slove the problem?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160805/d9466573/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/jpeg
Size: 109260 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160805/d9466573/attachment.jpe
More information about the ncl-talk
mailing list