[ncl-talk] Vertical cross section plot of WRF_Reflectivity (10 Cm) (REFL_10CM) through various longitude corresponding to a single lat
Bill Ladwig
ladwig at ucar.edu
Wed Sep 19 14:32:00 MDT 2018
Is there a reason why you can't use the cross section capability of
wrf_user_intrp3d? See Example 1 here:
https://www.ncl.ucar.edu/Document/Functions/WRF_arw/wrf_user_intrp3d.shtml
If you need specific height levels, rather than the 1% increments the
routine does by default, I can show you how to accomplish that (see code
below, or do a search, since I've guided others through this in the past).
Also, WRF uses PH + PHB to get the geopotential, not ZNW, so to get the
height at each grid point, you do (PH + PHB) / 9.81 (or just use
wrf_user_getvar with 'z' as the product type). (If you don't have both of
those variables in your WRF file, you're going to be in trouble). Also, the
right two dimensions in a WRF file aren't lat and lon, they're south_north,
and west_east, which is a bit of a misnomer since it's really used for
"bottom to top" and "left to right" for a grid. WRF uses curvilinear
coordinates, so the latitude and longitude values require two dimensions to
determine them (this is why XLAT and XLONG are two dimensions).
For reference, the code that computes the cross section in wrf_user_intrp3d
is here (removed the horizontal interpolation stuff):
undef("wrf_user_intrp3d")
> function wrf_user_intrp3d( var3d:numeric, z_in:numeric, \
> plot_type:string, \
> loc_param:numeric, angle:numeric, opts:logical
> )
> ; var3d - 3d field to interpolate (all input fields must be
> unstaggered)
> ; z_in - interpolate to this field (either p/z)
> ; plot_type - interpolate horizontally "h", or vertically "v"
> ; loc_param - level(s) for horizontal plots (eg. 500hPa ; 3000m -
> scalar),
> ; plane for vertical plots (2 values representing an xy point
>
> ; on the model domain through which the vertical plane will
> pass
> ; OR 4 values specifying start and end values
> ; angle - 0.0 for horizontal plots, and
> ; an angle for vertical plots - 90 represent a WE cross
> section
> ; opts Used IF opts is TRUE, else use loc_param and angle to
> determine crosssection
> begin
>
> if(plot_type .eq. "v" ) then ; vertical cross section needed
> dims = dimsizes(var3d)
> nd = dimsizes(dims)
> dimX = dims(nd-1)
> dimY = dims(nd-2)
> dimZ = dims(nd-3)
> if ( nd .eq. 4 ) then
> if ( z_in(0,dimZ-1,0,0) .lt. z_in(0,dimZ-2,0,0) ) then
> ; We must be interpolating to pressure
> ; This routine needs input field and level in hPa - lets make
> sure of this
> if ( max(z_in) .gt. 2000. ) then
> ; looks like we have Pa as input - make this hPa
> z_in = z_in * 0.01
> end if
> end if
> z = z_in(0,:,:,:)
> else
> if ( z_in(dimZ-1,0,0) .lt. z_in(dimZ-2,0,0) ) then
> ; We must be interpolating to pressure
> ; This routine needs input field and level in hPa - lets make
> sure of this
> if ( z_in(0,0,0) .gt. 2000. ) then
> ; looks like we have Pa as input - make this hPa
> z_in = z_in * 0.01
> end if
> end if
> z = z_in
> end if
>
> ; set vertical cross section
> if (opts) then
> xy = wrf_user_set_xy( z, loc_param(0), loc_param(1), \ ;
> assumes 0-based indexing in v6.5.0
> loc_param(2), loc_param(3), \
> angle, opts )
> else
> xy = wrf_user_set_xy( z, loc_param(0), loc_param(1), \ ;
> 0.0, 0.0, angle, opts )
> end if
> xp = dimsizes(xy)
>
>
> ; first we interp z
> var2dz = wrf_interp_2d_xy( z, xy)
> ; interp to constant z grid
> if(var2dz(0,0) .gt. var2dz(1,0) ) then ; monotonically decreasing
> coordinate
> z_max = floor(max(z)/10)*10 ; bottom value
> z_min = ceil(min(z)/10)*10 ; top value
> dz = 10
> nlevels = tointeger( (z_max-z_min)/dz)
> z_var2d = new( (/nlevels/), typeof(z))
> z_var2d(0) = z_max
> dz = -dz
> else
> z_max = max(z)
> z_min = 0.
> dz = 0.01 * z_max
> nlevels = tointeger( z_max/dz )
> z_var2d = new( (/nlevels/), typeof(z))
> z_var2d(0) = z_min
> end if
>
> do i=1, nlevels-1
> z_var2d(i) = z_var2d(0)+i*dz
> end do
>
> ; interp the variable
> if ( dimsizes(dims) .eq. 4 ) then
> var2d = new( (/dims(0), nlevels, xp(0)/), typeof(var2dz))
> do it = 0,dims(0)-1
> var2dtmp = wrf_interp_2d_xy( var3d(it,:,:,:), xy)
> do i=0,xp(0)-1
> var2d(it,:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i),
> z_var2d)
> end do
> end do
> var2d!0 = var3d!0
> var2d!1 = "Vertical"
> var2d!2 = "Horizontal"
> else
> var2d = new( (/nlevels, xp(0)/), typeof(var2dz))
> var2dtmp = wrf_interp_2d_xy( var3d, xy)
> do i=0,xp(0)-1
> var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d)
> end do
> var2d!0 = "Vertical"
> var2d!1 = "Horizontal"
> end if
>
>
> st_x = tointeger(xy(0,0)) ; + 1 (removed 1-based indexing in 6.5.0)
> st_y = tointeger(xy(0,1)) ; + 1
> ed_x = tointeger(xy(xp(0)-1,0)) ; + 1
> ed_y = tointeger(xy(xp(0)-1,1)) ; + 1
> if (opts) then
> var2d at Orientation = "Cross-Section: (" + \
> st_x + "," + st_y + ") to (" + \
> ed_x + "," + ed_y + ")"
> else
> var2d at Orientation = "Cross-Section: (" + \
> st_x + "," + st_y + ") to (" + \
> ed_x + "," + ed_y + ") ; center=(" + \
> loc_param(0) + "," + loc_param(1) + \
> ") ; angle=" + angle
> end if
> return(var2d)
> end if
>
>
> end
If you want to set your own levels, set z_var2d to the levels you want and
comment out the section of code that computes the 1% increments. You'll
want to rename this function and copy this code in to your script to use it.
Hope this helps,
Bill
On Wed, Sep 19, 2018 at 7:44 AM Trisanu Banik <baniktrisanu at gmail.com>
wrote:
> Hi
>
> I am still struggling with the plot of vertical cross section of
> reflectivity with respect to height and various lat keeping long constant
> and vice versa.
> During the search of cross sectional plots in web, I got few functions
> that might be helpful for generating vertical cross sectional plot with
> respect to height.
>
> But my data has four dimensions time, eta level, lat, long. If I consider
> the values for particular time it will become 3D i.e. level, lat, long. Now
> to convert eta level to pressure I have utilized "q_isobaric =
> int2p_n_Wrap(P_HYD,REFL_10CM,level,linlog,0) " this. Presently I am able to
> plot reflectivity w..r.t pressure for various lat long.
>
> Now the challenge is to convert vertical level to height. Various inbuilt
> functions are there in NCL that convert pressure to height like "wrf_user_vert_interp"
> and gsn_csm_pres_hgt etc. But my data dimension does not satisfy the
> functions requirement. In wrf_out file, the necessary variables related to
> this are, one dimentional vertical eta level (ZNW), Hydrostatic Pressure
> level (P_HYD) and reflectivity (REFL_10CM). No height related 3D parameter
> is available in the wrf_out file. Only terrain height is available which is
> 2D.
>
> For your convenience, I have attached my tentative code. Please guide me
> to plot the vertical cross section of reflectivity with Height in Y axis
> for various lat long.
>
> Thanks & Regards
> Trisanu
>
>
>
> *Trisanu Banik,PhD*
>
> *Research Scientist*
>
> *North Eastern Space Applications Centre (NESAC)*
> *Government of India*
> *Department of Space*
>
> *Umiam-793103, Meghalaya*
> *Mobile-9774837581*
>
>
>
> On Thu, Sep 13, 2018 at 7:58 PM Jim Means <jim at weatherextreme.com> wrote:
>
>> In addition to Soma's suggestion, you might also look at the cross
>> section plot examples here:
>>
>> http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/NCL_examples.htm
>>
>> On 9/13/2018 7:09 AM, Soma Roy wrote:
>>
>> See the different examples at the below link:
>> https://www.ncl.ucar.edu/Applications/height_long.shtml
>>
>> Best,
>> Soma
>>
>> On Thu, Sep 13, 2018, 19:10 Trisanu Banik <baniktrisanu at gmail.com> wrote:
>>
>>> Hi,
>>>
>>> I need help to plot a vertical cross section plot of WRF REFL_10CM
>>> (Reflectivity). The parameter is four dimensional (Time, eta level, lat,
>>> lon). I tried to plot the vertical cross section of reflectivity with
>>> respect to height for a particular time. But the REFL_10CM parameter
>>> contain only eta level. Therefore whenever I generated plot, it always
>>> provide vertical cross section w.r.t eta level only. i need suggestion how
>>> I can convert eta level to pressure level/Height level. So that, I can plot
>>> a vertical cross section of WRF reflectivity for various longitude
>>> corresponding to a single latitude. Here I have added a sample plot. I need
>>> to plot vertical cross section w.r.t height (not w.r.t eta level)
>>>
>>> Please help me to resolve the issue.
>>>
>>>
>>> *Trisanu Banik,PhD*
>>>
>>> *Research Scientist *
>>>
>>> *North Eastern Space Applications Centre (NESAC) *
>>> *Government of India*
>>> *Department of Space*
>>>
>>> *Umiam-793103, Meghalaya *
>>> *Mobile-9774837581*
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>
>>
>> _______________________________________________
>> ncl-talk mailing listncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>> --
>>
>> James D. Means, Ph.D.
>> Senior Atmospheric & Climate Scientist
>> California Office
>> Tele: 619-495-1638 | Fax: 775-636-8430
>> 930 Tahoe Blvd., Suite 802-560
>> Incline Village, Nevada 89451
>> jim at weatherextreme.com | vcard
>> <http://www.weatherextreme.com/vcards/James%20Means.vcf>
>> www.weatherextreme.com
>> <http://www.weatherextreme.com/>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
> _______________________________________________
> 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/20180919/d4b0def2/attachment.html>
More information about the ncl-talk
mailing list