NCL画WRF输出数据，限定经纬度范围后，等值线不对了

黄小仙儿

想把WRF输出的数据限定自己想要的经纬度范围画出来，看了NCL的例子，http://www.ncl.ucar.edu/Document ... _map_overlays.shtml设定了自己想要的区域，图是画出来了，但是等值线不对了。虽然底图区域缩小了，但是等值线并没有相应对应到地图上，求助~~我的脚本：
;   Example script to produce dbz plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   November 2008

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
   a = addfile("/home/Huanglei/data/d032"+".nc","r")

; We generate plots, but what kind do we prefer?
  type = "pdf"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_dbzmax1")
  gsn_define_colormap(wks,"WhViBlGrYeOrReWh")       ; Overwrite the standard color map

; Set some basic resources
  res = True
  res@MainTitle                   = "REAL-TIME WRF"

   pltres = True
pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove overlays
  mpres = True
mpres@mpOutlineOn = False  ; Turn off map outlines
mpres@mpFillOn    = False  ; Turn off map fill

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Which times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  do it = 1,ntimes-1               ; TIME LOOP
;it=20
     print("Working on time: " + times(it) )
     res@TimeLabel = times(it)   ; Set Valid time to use on plots


; First get the variables we will need        
; Both dbz and mdbz will be calculated using constant intercept parameters
; (as in early Reisner-2), with values of 8x10^6, 2x10^7, and 4x10^6 m^-4,
; for rain, snow, and graupel, respectively.

        ; mdbz = wrf_user_getvar(a,(/"mdbz","1","1"/),it)
     dbz = wrf_user_getvar(a,(/"dbz","1","1"/),it)
  lat2d = wrf_user_getvar(a,"XLAT",it)   ; latitude
  lon2d = wrf_user_getvar(a,"XLONG",it)  ; longitude
     mpres  = True
     opts = res   
          opts@sfXArray = lon2d
  opts@sfYArray = lat2d
     opts@cnFillOn = True  
     opts@ContourParameters = (/ 5., 65., 5./)
         opts@lbAutoManage= True
      opts@lbOrientation="Vertical"
  opts@pmLabelBarParallelPosF = 1.08 ;default = 0.5
  opts@pmLabelBarOrthogonalPosF =-0.9  ;lb的位置 default = 0.02
      opts@pmLabelBarWidthF      = 0.2               ; default is shorter
      opts@pmLabelBarHeightF     = 0.6               ; default is taller
      opts@lbTitleOn        = False                  ; turn on titl
          opts@lbLabelFontHeightF    = 0.000001             ; default is HUGE
          opts@lbTitlePosition  =  "Right"
          opts@lbTitleDirection = "Down"
      opts@mpLimitMode = "LatLon"    ; Limit the map view.
      mpres@mpLeftCornerLatF  =  29.32
      mpres@mpRightCornerLatF =  32.
      mpres@mpLeftCornerLonF  = 102.45
      mpres@mpRightCornerLonF =105.57

      ;opts@lbTitleOn="False"
     ;contour = wrf_contour(a,wks,dbz(1,:,:),opts)     ; plot only lowest level
    ; plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

     contour = wrf_contour(a,wks,dbz(1,:,:),opts)
     plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

;>============================================================<
;                      add China map
;>------------------------------------------------------------<


  shp_name1    = "/home/Huanglei/map/China/diquJie_polyline.shp"

  lnres                  = True
  lnres@gsLineColor      = "gray25"
  lnres@gsLineThicknessF = 0.5   

id = gsn_add_shapefile_polylines(wks,plot,shp_name1,lnres)
  shp_name2    = "/home/Huanglei/map/China/cnmap/cnhimap.shp"

  prres=True
  prres@gsLineThicknessF = 2.0      
  prres@gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)
  draw(plot)
  frame(wks)
delete(opts)
   end do        ; END OF TIME LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end

画出的图：

askira2011

貌似你改了区域后只有地图被放大了，等值线没变，你可以试试这个http://www2.mmm.ucar.edu/wrf/OnL ... PECIAL/wrf_Zoom.htm

黄小仙儿

askira2011 发表于 2014-12-3 11:28
貌似你改了区域后只有地图被放大了，等值线没变，你可以试试这个http://www2.mmm.ucar.edu/wrf/OnLineTutor ...

好的，谢谢，我看了看那个
  dims = dimsizes(ter)
    ; As an example, we are looking for the lower right 1/4 of the domain
    x_start = dims(1)/2
    x_end   = dims(1)-1
    y_start = 0
    y_end   = dims(0)/2
这个x_start ，x_end 。。。是格点数吗？

askira2011

黄小仙儿 发表于 2014-12-3 11:44
好的，谢谢，我看了看那个
  dims = dimsizes(ter)
    ; As an example, we are looking for the lowe ...

对的  这个要你自己算

黄小仙儿

askira2011 发表于 2014-12-3 12:09
对的  这个要你自己算

恩，好的 多谢~