<DIV> Hi:</DIV>
<DIV> 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.</DIV>
<DIV> But after I compared the output data with the contour plot that I ploted ,them are not <SPAN closure_uid_407902488="7">corresponding.</SPAN></DIV>
<DIV><SPAN closure_uid_407902488="7"> I thought maybe I confused the lat and lon ,but after I switch it ,the output data are not <SPAN closure_uid_407902488="7">correspond to the plots too.</SPAN></SPAN></DIV>
<DIV><SPAN closure_uid_407902488="7"><SPAN closure_uid_407902488="7"> I can not figure out where the problem is? I think maybe it is because of the zoom ?</SPAN></SPAN></DIV>
<DIV><SPAN closure_uid_407902488="7"><SPAN closure_uid_407902488="7"> But I do not know how to fix it,can you guys fix this?</SPAN></SPAN></DIV>
<DIV> </DIV>
<DIV><BR></DIV>
<DIV>This is my script:</DIV>
<DIV>
<DIV>; Example script to produce plots for a WRF real-data run,<BR>; with the ARW coordinate dynamics option.</DIV>
<DIV>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" <BR>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" <BR>load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"</DIV>
<DIV>begin<BR>;<BR>; The WRF ARW input file. <BR>; This needs to have a ".nc" appended, so just do it.<BR> a = addfile("/backup/hl/z_20min/wrfout_d03_2010-08-11_00:00:00"+".nc","r") </DIV>
<DIV><BR>; We generate plots, but what kind do we prefer?<BR>; type = "eps"<BR> type = "pdf"<BR>; type = "ps"<BR>; type = "ncgm"<BR> wks = gsn_open_wks(type,"plt_Cloudgraup_20min_20100811")<BR> ;colors = (/"white","black","white","royal blue","light sky blue",\<BR> ; "powder blue","light sea green","pale green","wheat","brown",\<BR> ; "pink"/)<BR> gsn_define_colormap(wks,"BlAqGrYeOrRevi200") ; overwrite the .hluresfile color map<BR><BR>; Set some basic resources<BR> res = True<BR> <A href="mailto:res@MainTitle">res@MainTitle</A> = "REAL-TIME WRF"</DIV>
<DIV> mpres = True ; Map resources<BR> <A href="mailto:mpres@mpOutlineOn">mpres@mpOutlineOn</A> = False ; Turn off map outlines<BR> <A href="mailto:mpres@mpFillOn">mpres@mpFillOn</A> = False ; Turn off map fill<BR> <A href="mailto:mpres@mpGridAndLimbOn">mpres@mpGridAndLimbOn</A> = True<BR> ;res@mpProjection = "Lambert"<BR> pltres = True ; Plot resources<BR> <A href="mailto:pltres@PanelPlot">pltres@PanelPlot</A> = True ; Tells wrf_map_overlays not to remove overlays</DIV>
<DIV><BR>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<BR>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</DIV>
<DIV>; What times and how many time steps are in the data set?<BR> times = wrf_user_getvar(a,"times",-1) ; get all times in the file<BR> ntimes = dimsizes(times) ; number of times in the file<BR> lat = new(ntimes,float) ;creat new variable to store the lat of the max qi<BR> lon = new(ntimes,float)<BR>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</DIV>
<DIV> do it =0, ntimes-1,1 ; TIME LOOP</DIV>
<DIV> print("Working on time: " + times(it) )<BR> <A href="mailto:res@TimeLabel">res@TimeLabel</A> = times(it) ; Set Valid time to use on plots<BR>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<BR>; First get the variables we will need <BR> if(isfilevar(a,"QGRAUP"))<BR> qi = wrf_user_getvar(a,"QGRAUP",it )<BR> qi = qi*1000.<BR> <A href="mailto:qi@units">qi@units</A> = "g/kg" <BR> ; qiav = dim_avg_n_Wrap(qi, 0)<BR> end if<BR> tes=True<BR> <A href="mailto:tes@returnInt">tes@returnInt</A> = False<BR> loc1=wrf_user_ll_to_ij(a,105.0,32.0,tes)<BR> ;print("X/Y location is: "+ loc1)<BR> <BR> loc2=wrf_user_ll_to_ij(a,112.0,38.0,tes)<BR> ; print("X/Y location is: "+ loc2)<BR> x_start = 15<BR> x_end = 227<BR> y_start = 71<BR> y_end = 274</DIV>
<DIV> level = 13</DIV>
<DIV> display_level = level + 2<BR> opts = res<BR> <A href="mailto:opts@cnFillOn">opts@cnFillOn</A> = True<BR> <A href="mailto:opts@gsnSpreadColors">opts@gsnSpreadColors</A> = False<BR> <A href="mailto:opts@ContourParameters">opts@ContourParameters</A> = (/ 1, 5, 0.5 /)<BR> <A href="mailto:opts@PlotLevelID">opts@PlotLevelID</A> = " Level -10 - -20 ~S~o~N~C " <BR> <A href="mailto:opts@gsnDraw">opts@gsnDraw</A> = False <BR> <A href="mailto:opts@gsnFrame">opts@gsnFrame</A> = False</DIV>
<DIV><BR> if (isvar("qi"))<BR> qis = qi(level + (/0,1,2,3,4,5,6/),:,:)<BR> qisum = dim_sum_n_Wrap(qis, 0)<BR> mpres1 = True<BR> <A href="mailto:mpres1@ZoomIn">mpres1@ZoomIn</A> = True<BR> <A href="mailto:mpres1@Xstart">mpres1@Xstart</A> = x_start<BR> <A href="mailto:mpres1@Ystart">mpres1@Ystart</A> = y_start<BR> <A href="mailto:mpres1@Xend">mpres1@Xend</A> = x_end<BR> <A href="mailto:mpres1@Yend">mpres1@Yend</A> = y_end<BR> ; printVarSummary(qisum)<BR> qi_zoom = qisum(y_start:y_end,x_start:x_end)<BR> ; printVarSummary(qi_zoom)<BR> ;;;;;;;convert to 1D;;;;;;;;<BR> qi_zoom1D = ndtooned(qi_zoom)<BR> dsizes_a = dimsizes(qi_zoom)<BR> ;;;;;;resolve the 1D indices back to their original 2D arry.;;;<BR> indices = ind_resolve(maxind(qi_zoom1D),dsizes_a)<BR> ; print(indices)<BR> 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 <BR> lon(it)= ilatlon(0)<BR> lat(it)= ilatlon(1)<BR> <BR> contour = wrf_contour(a,wks,qi_zoom,opts)<BR> plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres1)</DIV>
<DIV> delete(contour)<BR> end if<BR>;>============================================================<<BR>; add map<BR>;>------------------------------------------------------------<</DIV>
<DIV> <BR> shp_name2 = "/public/home/huanglei/map/shaanxi_city_l.shp"</DIV>
<DIV> prres=True<BR> <A href="mailto:prres@gsLineThicknessF">prres@gsLineThicknessF</A> = 1.0 <BR> <A href="mailto:prres@gsLineColor">prres@gsLineColor</A> = "black"<BR> plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)</DIV>
<DIV> txres2 = True<BR> <A href="mailto:txres2@txFont">txres2@txFont</A> = 10<BR> <A href="mailto:txres2@txFontHeightF">txres2@txFontHeightF</A> =0.01<BR> <A href="mailto:txres2@txFontColor">txres2@txFontColor</A> = "Blue"<BR> txdum1 =gsn_add_text(wks, plot, "Xian", 108.93,34.27, txres2)</DIV>
<DIV> draw(plot) ; This will draw the map and the shapefile outlines.<BR> frame(wks)<BR> delete(opts)</DIV>
<DIV> </DIV>
<DIV>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;</DIV>
<DIV> end do ; END OF TIME LOOP <BR> <BR> output_file = "/public/home/huanglei/z_data/20100811graup.txt"<BR> write_table(output_file,"w",[/times,lon,lat/],"%s%9.4f%9.4f")<BR>end<BR></DIV></DIV>
<DIV><BR></DIV>
<DIV><SPAN style="LINE-HEIGHT: 23px; FONT-FAMILY: 'lucida Grande', Verdana, 'Microsoft YaHei'"> How can I slove the problem?</SPAN></DIV>
<DIV></DIV>