[ncl-talk] used the funtion " ind_resolve", "maxind" and "wrf_user_ij_to_ll" to find out the lat and lon location of the max data, but the output data are not correspond to the plots

[grace]

Hi:
        All,I am trying to find out the lat and lon location of the max graup of wrfout data,I have used the funtion " ind_resolve","maxind" and "wrf_user_ij_to_ll" according to the web examples.
        But after I compared the output data with the contour plot that I ploted ,them are not corresponding.
        I thought maybe I confused the lat and lon ,but after I switch it ,the output data are not correspond to the plots too.
        I can not figure out where the problem is? I think maybe it is because of the zoom ?
        But I do not know how to fix it,can you guys fix this?
  
 

 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/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 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?
  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 at 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 at units = "g/kg"   
         ; qiav = dim_avg_n_Wrap(qi, 0)
    end if
   tes=True
   tes at 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 at cnFillOn         = True
      opts at gsnSpreadColors  = False
      opts at ContourParameters       = (/ 1, 5, 0.5 /)
      opts at PlotLevelID      = " Level -10 - -20 ~S~o~N~C " 
          opts at gsnDraw      =  False                  
      opts at 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 at ZoomIn = True
        mpres1 at Xstart = x_start
        mpres1 at Ystart = y_start
        mpres1 at Xend   = x_end
        mpres1 at 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 at gsLineThicknessF = 1.0       
  prres at gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)
    txres2  = True
   txres2 at txFont  = 10
   txres2 at txFontHeightF =0.01
   txres2 at 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


 

  How can I slove the problem?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160927/e31358b2/attachment.html