[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
Bill Ladwig
ladwig at ucar.edu
Tue Sep 27 10:20:49 MDT 2016
Hi Grace,
One thing that is unusual (and confusing) about the wrf_user_ll_to_ij and
wrf_user_ij_to_ll routines is that they use Fortran 1-based indexing,
rather than NCL's 0-based indexing. So, when you get the result from
wrf_user_ll_to_ij, you need to subtract 1 from that result. When you use
wrf_user_ij_to_ll, you need to add 1 to the i,j arguments. At a first
glance at your script, it doesn't appear you are doing this, so I'd try
that first and see if that helps.
Bill
On Mon, Sep 26, 2016 at 9:24 PM, grace <313695096 at qq.com> wrote:
> 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?
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160927/f8e01033/attachment.html
More information about the ncl-talk
mailing list