[ncl-talk] data did not match with original plot after mask data with shapefile
Mary Haley
haley at ucar.edu
Thu Aug 11 02:02:04 MDT 2016
Hi,
If you continue to struggle with the shapefile_mask_data function, then it
would really help if you could provide the latest version of your script,
along with all the required data files. I'm the person who wrote this
function, so I can help you debug the problem pretty quickly.
You can use our ftp site, if your data files aren't too large:
http://www.ncl.ucar.edu/report_bug.shtml#HowToFTP
--Mary
On Fri, Aug 5, 2016 at 11:36 AM, Kim, Changhwan <
changhwan.kim at metoffice.gov.uk> wrote:
> Hi,
>
>
>
> I don’t know exactly what reason is cause that problem.
>
> I recommend to check coordinate variables of rain_tot_mask(using
> printVarSummary() function).
>
> Maybe those are absent.
>
>
>
> Or use this command before plotting script.
>
>
>
> Copy_varCoords(rain_tot_tend, rain_tot_mask)
>
> or
>
> Copy_varMeta(rain_tot_tend, rain_tot_mask)
>
>
>
> Cheers!
>
>
>
> *From:* ncl-talk-bounces at ucar.edu [mailto:ncl-talk-bounces at ucar.edu] *On
> Behalf Of *grace
> *Sent:* 05 August 2016 10:14
> *To:* ncl-talk
> *Subject:* [ncl-talk] data did not match with original plot after mask
> data with shapefile
>
>
>
> 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?
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160811/357ade77/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.jpg
Type: image/jpeg
Size: 109260 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160811/357ade77/attachment-0001.jpg
More information about the ncl-talk
mailing list