用NCL输出最大值的经纬度，但是输出的数据和我画出来的平面图不匹配

黄小仙儿

PS：刚刚发错版了。。发到GRADS那去了，不知道怎么把那边的删掉，只有跑到这里再发一次了。sorry啊版主
用NCL输出最大值的经纬度，但是输出的数据和我画出来的平面图不匹配，差的很多。
每个时次都对比了，但是都不匹配，差的距离有点大，不只是一两个格点的误差问题。
应该是我的程序有错，但是看不出是哪里出错了，因为我画图之前把数据zoom了一下，不知道是不是因为这个原因。但是不知道怎么改，请各位帮忙看看

;   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/csm/gsn_csm.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("/backup/hl/z_20min/wrfout_d03_2010-08-11_00:00:00"+".nc","r")            

; We generate plots, but what kind do we prefer?
;  type = "eps"
type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_Cloudgraup_20min_20100811")
   ;colors = (/"white","black","white","royal blue","light sky blue",\
           ;  "powder blue","light sea green","pale green","wheat","brown",\
            ; "pink"/)
             gsn_define_colormap(wks,"BlAqGrYeOrRevi200")   ; overwrite the .hluresfile color map

; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"
  mpres  = True  ; Map resources
  mpres@mpOutlineOn = False  ; Turn off map outlines
  mpres@mpFillOn    = False  ; Turn off map fill
  mpres@mpGridAndLimbOn = True
;res@mpProjection          = "Lambert"
  pltres = True ; Plot resources
  pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove overlays

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What 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
lat = new(ntimes,float)  ;creat new variable to store the lat of the max qi
lon = new(ntimes,float)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  do it =0, ntimes-1,1        ; TIME LOOP
    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
    if(isfilevar(a,"QGRAUP"))
          qi = wrf_user_getvar(a,"QGRAUP",it )
      qi = qi*1000.
      qi@units = "g/kg"   
         ; qiav = dim_avg_n_Wrap(qi, 0)
    end if
   tes=True
   tes@returnInt = False
   loc1=wrf_user_ll_to_ij(a,105.0,32.0,tes)
   ;print("X/Y location is: "+ loc1)
   
   loc2=wrf_user_ll_to_ij(a,112.0,38.0,tes)
  ; print("X/Y location is: "+ loc2)
   x_start = 15
   x_end   = 227
   y_start = 71
   y_end   = 274
     level = 13
      display_level = level + 2
      opts = res
      opts@cnFillOn         = True
      opts@gsnSpreadColors  = False
      opts@ContourParameters       = (/ 1, 5, 0.5 /)
      opts@PlotLevelID      = " Level -10 - -20 ~S~o~N~C "
          opts@gsnDraw      =  False                  
      opts@gsnFrame     =  False

      if (isvar("qi"))
        qis  = qi(level + (/0,1,2,3,4,5,6/),:,:)
        qisum = dim_sum_n_Wrap(qis, 0)
        mpres1 = True
        mpres1@ZoomIn = True
        mpres1@Xstart = x_start
        mpres1@Ystart = y_start
        mpres1@Xend   = x_end
        mpres1@Yend   = y_end
      ;  printVarSummary(qisum)
        qi_zoom = qisum(y_start:y_end,x_start:x_end)
       ; printVarSummary(qi_zoom)
       ;;;;;;;convert to 1D;;;;;;;;
        qi_zoom1D  = ndtooned(qi_zoom)
        dsizes_a = dimsizes(qi_zoom)
      ;;;;;;resolve the 1D indices back to their original 2D arry.;;;
        indices = ind_resolve(maxind(qi_zoom1D),dsizes_a)
       ; print(indices)
        ilatlon = wrf_user_ij_to_ll(a, indices(0,1),indices(0,0),True)  ; select the latitude index where the X array is at its' maximum
        lon(it)= ilatlon(0)
        lat(it)= ilatlon(1)
                                      
        contour = wrf_contour(a,wks,qi_zoom,opts)
        plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres1)
        delete(contour)
      end if
;>============================================================<
;                      add  map
;>------------------------------------------------------------<
     
  shp_name2    = "/public/home/huanglei/map/shaanxi_city_l.shp"
  prres=True
  prres@gsLineThicknessF = 1.0      
  prres@gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)
   txres2  = True
   txres2@txFont  = 10
   txres2@txFontHeightF =0.01
   txres2@txFontColor = "Blue"
   txdum1 =gsn_add_text(wks, plot, "Xian", 108.93,34.27, txres2)
  draw(plot)       ; This will draw the map and the shapefile outlines.
  frame(wks)
   delete(opts)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  end do        ; END OF TIME LOOP     

output_file = "/public/home/huanglei/z_data/20100811graup.txt"
write_table(output_file,"w",[/times,lon,lat/],"%s%9.4f%9.4f")
end

台风三维可视化

你好，我有问题也要咨询你，给你的新浪微博发了私信，不知看到没呢？