[ncl-talk] data did not match with original plot after mask data with shapefile

Kim, Changhwan changhwan.kim at metoffice.gov.uk
Fri Aug 5 03:36:29 MDT 2016


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.
[cid:image001.jpg at 01D1EF04.A60EF200]
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/c822f2f0/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.jpg
Type: image/jpeg
Size: 109260 bytes
Desc: image001.jpg
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160805/c822f2f0/attachment-0001.jpg 


More information about the ncl-talk mailing list