[ncl-talk] error when mask data with shapefile according to the example on the website.

grace 313695096 at qq.com
Thu Aug 4 21:02:08 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 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/4daf31ba/attachment.html 


More information about the ncl-talk mailing list