[ncl-talk] error when mask data with shapefile according to the example on the website.
Kim, Changhwan
changhwan.kim at metoffice.gov.uk
Fri Aug 5 02:26:55 MDT 2016
Hi,
The first parameter of shapefile_mask_data() function must have coordinate informations.
You should add coordinate variables to rain_tot_tend.
For example,
rain_tot_tend!0 = "lat2d"
rain_tot_tend!1 = "lon2d"
rain_tot_tend at lat2d = grid_y
rain_tot_tend at lon2d = grid_x
Cheers.
From: ncl-talk-bounces at ucar.edu [mailto:ncl-talk-bounces at ucar.edu] On Behalf Of grace
Sent: 05 August 2016 04:02
To: ncl-talk
Subject: [ncl-talk] error when mask data with shapefile according to the example on the website.
Hi:
All,I am trying to mask data with shapefile according to the example shapefiles_14_mask.ncl on the NCL website.
But it appears error: shapefile_mask_data: Error: not a valid rectilinear, curvilinear, or unstructured grid.
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)"
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/fbd25b1c/attachment.html
More information about the ncl-talk
mailing list