Re: Transect

The WRF grid is curvilinear ... you must interpolate
to a rectilinear grid prior to using the transect code.

 
  f = addfile ("some_RCM_file.nc", "r")
   lat2d = f->XLAT(0,:,:) ; size = (nlat,nlon)
   lon2d = f->XLON(0.:,:) ; size = (nlat,nlon)
   x = f->X

   xgrd = rcm2rgrid(lat2d,lon2d,x,lat,lon,0)

then do the transect

Good luck

On Dec 11, 2008 1:54:46 PM, gs6r_at_virginia.edu wrote:
> Hey
>
> I have tried to plot a transect of temperature profile. It is very similar to the example you give http://www.ncl.ucar.edu/Applications/transect.shtml
>
> However, for some reasons, my plot does not look right.
> I think the reason is that I cannot specify the vertical heights correctly. I simply use the height of one grid at a specific time to represent the one dimension variable height for all, which I guess makes my plot incorrect.
>
> Attached is the code I wrote, could you please help me to correct it? Thanks a lot,
>
> Guan Song
> Graduate Student
> Reply email: gs6r_at_virginia.edu
> New Clark 272
> Dept. of Environmental Sciences
> University of Virginia
> ;************************************
> ; CSM_Graphics: transect_1.ncl
> ;************************************
> 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/WRF_contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>
> ;************************************
> begin
> wks = gsn_open_wks("ps","trans30") ; open a ps file
> gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; choose colormap
>
> in = addfile ("/home/gs6r/test1/WRFV3/test/em_real/RUC/wrfout_d03_2007-01-31_00:00:00.nc", "r")
> ; What times and how many time steps are in the data set?
> times = wrf_user_list_times(in) ; get times in the file
> ntimes = dimsizes(times) ; number of times in the file
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> do ti = 0,ntimes-1 ; TIME LOOP
> T = in->T(ti,:,:,:)
> P = in->P(ti,:,:,:)
> TK = wrf_tk(P,T)
> lat = in->XLAT(0,:,0)
> lon = in->XLONG(0,0,:)
> z = wrf_user_getvar(in,"z",5)
> ;************************************
> ; calculate great circle along transect
> ;************************************
> leftlat = 40.654
> rightlat = 40.654
> rightlon = -111.65
> leftlon = -112.1500
>
> npts = 10 ; number of points in resulting transect
> dist = gc_latlon(leftlat,leftlon,rightlat,rightlon,npts,-4)
> points = ispan(0,npts-1,1)*1.0
> ;********************************
> ; interpolate data to great circle
> ;*******************************
> TK!0 = "z"
> TK!1 = "south_north"
> TK!2 = "west_east"
> TK&south_north=lat
> TK&west_east=lon
> TK&z=z(:,5,5) ; Possible WRONG,BUT NOT SURE HOW TO CORRECT IT
> trans = linint2_points(TK&west_east,TK&south_north,TK,True,dist_at_gclon,dist_at_gclat,2)
> copy_VarAtts(TK,trans) ; copy attributes
> trans!0 = "z" ; create named dimension and assign
> trans&z = TK&z ; coordinate variable for 0th dimension only
> ;********************************
> ; create plot
> ;********************************
> res = True ; plot mods desired
> res_at_tmXBMode = "Explicit" ; explicitly label x-axis
> res_at_tmXBValues = (/points(0),points(npts-1)/) ; points to label
> ; label values
> res_at_tmXBLabels = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)
> res_at_gsnYAxisIrregular2Linear = True
> res_at_cnFillOn = True ; turn on color
> res_at_lbLabelAutoStride = True ; nice label bar label stride
> res_at_gsnSpreadColors = True ; use full range of colormap
> res_at_cnLinesOn = False ; turn off countour lines
> res_at_lbOrientation = "vertical" ; vertical label bar
> res_at_pmLabelBarOrthogonalPosF = -0.05 ; move label bar closer to plot
> res_at_cnLevelSelectionMode = "ManualLevels" ; Manually set contour levels
> res_at_cnMinLevelValF = -5.0
> res_at_cnMaxLevelValF = 5.0
> res_at_cnLevelSpacingF = 1.0
> res_at_lbLabelAutoStride = True ; Control labelbar labels.
> res_at_tiMainString = "Transect" ; add title
> res_at_tiXAxisString = "lat/lon along transect"
> res_at_trYReverse = False ; reverse y axis
> ; res_at_trXReverse = True ; reverse x axis (neg longitudes)
> res_at_cnLevelSpacingF = 0.2 ; set contour spacing
>
> plot = gsn_csm_contour(wks,trans,res) ; create plot
> ;********************************
> ; show transect on a map
> ;********************************
> i = NhlNewColor(wks,0.8,0.8,0.8) ;add gray to continents
> mres = True ; plot mods desired
> mres_at_gsnFrame = False ; don't turn page yet
> mres_at_gsnDraw = False ; don't draw yet
> mres_at_tiMainString= "Transect Location" ; title
> ; map = gsn_csm_map_ce(wks,mres) ; create map
> map=wrf_map(wks,in,mres)
> ; add polyline
> pres = True ; polyline mods desired
> pres_at_gsLineColor = "red" ; color of lines
> pres_at_gsLineThicknessF = 2.0 ; line thickness
> gsn_polyline(wks,map,(/leftlon,rightlon/),(/leftlat,rightlat/),pres)
> gsres = True
> gsres_at_gsMarkerColor="NavyBlue"
> gsres_at_gsMarkerIndex=12
> gsres_at_gsMarkerSizeF=10
> gsres_at_gsMarkerThicknessF=3.0
> dum=gsn_add_polymarker(wks,map,-111.888,40.654,gsres)
> gsres_HDP=True
> gsres_HDP_at_gsMarkerColor="Red"
> gsres_HDP_at_gsMarkerIndex=5
> gsres_HDP_at_gsMarkerSizeF=10
> gsres_HDP_at_gsMarkerThicknessF=3.0
> dum_HDP=gsn_add_polymarker(wks,map,-111.65,40.56,gsres_HDP)
> draw(map)
> frame(wks)
> end do
> end
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Dec 14 2008 - 14:01:15 MST