Wrong wind direction of vertcross

[zxdawn]

Background

Hi, all,

I'm trying to plot the wind quiver on the crosssection. But, I find the wind direction is wrong. Here're the details.

Read data and vertcross

import numpy as np
from wrf import getvar, interplevel, ALL_TIMES, CoordPair, vertcross, to_np, latlon_coords
from netCDF4 import Dataset
import matplotlib.pyplot as plt

data_dir = '../data/wrfchem/wrfout/20200901/lnox_500/'
file = data_dir+'wrfout_d04_2020-09-01_05:00:00_subset'
ncfile = Dataset(file)

# set cross
start_lon = 118.6
end_lon = 119
start_lat = 32.1
end_lat = 31.7

start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)

# get vars
u = getvar(ncfile, 'ua')
v = getvar(ncfile, 'va')
w = getvar(ncfile, 'wa')
wspd = getvar(ncfile, 'uvmet_wspd_wdir')
dbz = getvar(ncfile, 'dbz')
z = getvar(ncfile, 'z', units='km')

# vertcross and subset to 0-17 km
cross_u = vertcross(u, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True).sel(vertical=slice(0, 17))
cross_v = vertcross(v, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True).sel(vertical=slice(0, 17))
cross_w = vertcross(w, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True).sel(vertical=slice(0, 17))
cross_dbz = vertcross(dbz, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True).sel(vertical=slice(0, 17))
cross_wspd = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True).sel(vertical=slice(0, 17))

# get xlabels
coord_pairs = to_np(cross_dbz.coords["xy_loc"])
xticks = np.arange(coord_pairs.shape[0])
xlabels = [pair.latlon_str() for pair in to_np(coord_pairs)]

fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xticks(xticks[::20])
ax.set_xticklabels(['(' +s.replace(',', '\n') + ')' for s in xlabels[::20]], fontsize=4) # rotation

# plot wind
xs = np.arange(0, cross_w.shape[-1], 1)
ys = to_np(cross_w.coords['vertical'])
step = 2

# plot dbz
m = ax.contourf(xs[::step],
                ys[::step],
                to_np(cross_dbz[::step, ::step]),
                cmap='Reds'
               )

q = ax.quiver(xs[::step],
              ys[::step],
              to_np(cross_u[::step, ::step]+cross_v[::step, ::step]),
              to_np(cross_w[::step, ::step])
             )

plt.colorbar(m, label='dBZ')
plt.savefig('./cross.jpg', dpi=300)

Result:

Wind interpolated to 15 km

# Get the lat/lon coordinates
lats, lons = latlon_coords(u)

zlevel = 15  # km
u_zlevel = interplevel(u, z, zlevel)
v_zlevel = interplevel(v, z, zlevel)

step = 10

fig, ax = plt.subplots(figsize=(10, 6))
q = plt.quiver(to_np(lons[::step,::step]), to_np(lats[::step,::step]),
              to_np(u_zlevel[::step, ::step]), to_np(v_zlevel[::step, ::step])
             )

qk = plt.quiverkey(q, 1, 1.05, 10, r'$10 m/s}$', labelpos='E')
plt.plot([start_lon, end_lon], [start_lat, end_lat], '-r')
plt.title('Wind at 15 km level')
plt.savefig('./zlev.jpg', dpi=300)

Problem

As you can see, the wind at the 15 km level should blow from the start point to the endpoint. But, the 15 km wind in the crosssection is the opposite.

Is there anything wrong with my method? Thanks in advance!

[michaelavs]

Hi @zxdawn, My first thought with this issue is that the two graphs are looking at difference perspectives of the data. The first plot (with dBZ overlay) uses vertcross while the second plot uses interplevel. If my understanding is correct, in the first plot, you will see data that is above the ground (i.e. looking out at the horizon) whereas the second plot you are looking down from above (i.e. looking from the z coordinate down at the xy-axis). So I think the reason for the differences comes from the perspective of the plots and not from an issue with your code.

However, I would definitely be interested to know if you have found a solution to this that isn't the result of what I have outlined above.

[zxdawn]

Hi @michaelavs , thanks for your explanation.

If I understand correctly, the first perspective is looking from the start point to the end point at each z level, right? Then, according to the second plot, the wind component on the cross line at 15 km should blow from the start point to the end point. But, it's opposite in the first plot.

Did I misunderstand anything?

[michaelavs]

Hi @zxdawn, The second plot is looking at the data from a different view point than the first. If you look at your code for the first plot, you use

q = ax.quiver(xs[::step],
              ys[::step],
              to_np(cross_u[::step, ::step]+cross_v[::step, ::step]),
              to_np(cross_w[::step, ::step])
             )

which is adding the u and v components together and plotting them with the w component. I am curious as to why you are adding the u and v components with this diagram. From my understanding, they should be orthogonal to each other and, based on the type of plot used, the v wind vector component would not be viewable from this orientation (this plot is looking in the xz or uw plane). With that in mind, when looking at the next plot, you use

q = plt.quiver(to_np(lons[::step,::step]), to_np(lats[::step,::step]),
              to_np(u_zlevel[::step, ::step]), to_np(v_zlevel[::step, ::step])
             )

which no longer uses the w component of the wind vector field, and further, is being viewed from the w (or z) plane, so the vectors shown will not have the w component effecting their position as we see it did in the first plot.

Now, that is all based on just vector calculus/vector fields. I would like to see if the following changes to your code yield the results you are looking for:

# Plot 1
q = ax.quiver(xs[::step],
              ys[::step],
              to_np(cross_u[::step, ::step]+cross_v[::step, ::step]),
              to_np(cross_w[::step, ::step])
             )

# Plot 2
q = plt.quiver(to_np(lons[::step,::step]), to_np(lats[::step,::step]),
              to_np(u_zlevel[::step, ::step]), to_np(w_zlevel[::step, ::step])
             )

# Or the other version could be

q = ax.quiver(xs[::step],
              ys[::step],
              to_np(cross_u[::step, ::step],
              cross_v[::step, ::step]))
             )

# Plot 2
q = plt.quiver(to_np(lons[::step,::step]), to_np(lats[::step,::step]),
              to_np(u_zlevel[::step, ::step]), to_np(w_zlevel[::step, ::step])
             )

You may have to change your code a little bit to get these two segments to work, but in theory, it should create the matching plots you were looking for. I unfortunately have not been able to test your code fully (due to not having the same wrfout file). But from what I have been able to run so far, and see from above, it looks to be a very solid code so I don't think you have any major issues overall. If this does not help with your issue, I think I may need to see an example of what you are trying to accomplish as an end result, and then I can work backwards from there.

Thank you!

[zxdawn]

Hi @michaelavs ,

Thanks for your suggestions! I suppose you misunderstood the situation.

I want to combine u,v, and w on the crosssection in the first plot. Here's the illustration:

That's why I get the u_cross and v_cross first and then add them together. Then, the u_cross+v_cross should be equal to u_zlevel+v_zlevel. Since the w only affect the angle, the general direction of cross vector should be similar to u_zlevel+w_zlevel.

I checked the u_cross and v_cross in detail:

$ cross_v.sel(vertical=15, method='nearest')

array([-4.5392895, -4.4891686, -4.4133954, -4.3930106, -4.4239078,
       -4.4032164, -4.3887043, -4.4004874, -4.40285  , -4.34502  ,
       -4.271366 , -4.1980615, -4.099139 , -3.989687 , -3.9039516,
       -3.881829 , -3.8909402, -3.9069679, -3.9016984, -3.8796399,
       -3.8316743, -3.7652757, -3.7081583, -3.6756363, -3.5856895,
       -3.4674985, -3.4614208, -3.535829 , -3.5605037, -3.4509377,
       -3.3177335, -3.2450762, -3.1743526, -3.0860896, -2.9951584,
       -2.8961744, -2.8069835, -2.789996 , -2.8981721, -3.0832877,
       -3.0994935, -3.0342689, -3.2060142, -3.497759 , -3.7033718,
       -3.61847  , -3.2701426, -2.781209 , -2.327874 , -2.0972714,
       -2.177136 , -2.4467058, -2.8243783, -3.2627208, -3.5687144,
       -3.6456037, -3.5517495, -3.436864 , -3.3582752], dtype=float32)

$ cross_u.sel(vertical=15, method='nearest')
array([ 1.0334523e+00,  1.0291141e+00,  1.0432535e+00,  1.0845915e+00,
        1.0369910e+00,  9.2963511e-01,  8.2718337e-01,  7.2292691e-01,
        6.1574626e-01,  5.6084716e-01,  6.3430572e-01,  7.6110315e-01,
        8.7388504e-01,  9.0909505e-01,  9.3065763e-01,  1.0016983e+00,
        1.0971006e+00,  1.1544362e+00,  1.1085669e+00,  9.9084026e-01,
        7.7602780e-01,  5.4501384e-01,  3.5534182e-01,  2.2430542e-01,
        1.0720103e-01,  9.0908118e-02,  2.9955885e-01,  5.9632653e-01,
        7.3460537e-01,  6.4204049e-01,  5.4950768e-01,  4.8518804e-01,
        3.6334375e-01,  2.1753691e-01,  8.7367281e-02, -8.0834376e-04,
       -8.3125807e-02, -1.9285528e-01, -2.9180545e-01, -3.2306853e-01,
       -2.7556732e-01, -4.7320965e-01, -8.4909719e-01, -1.1045828e+00,
       -1.1837103e+00, -1.1230041e+00, -1.0139287e+00, -8.8813722e-01,
       -7.6439202e-01, -6.8525493e-01, -6.2830311e-01, -6.6367453e-01,
       -7.2526580e-01, -8.1777006e-01, -9.2741233e-01, -1.0539831e+00,
       -1.1154382e+00, -1.1736290e+00, -1.2590120e+00], dtype=float32)

From the cross-section view, the v_cross should be positive. However, that's negative. If I multiple it by -1, then the 15km wind of the first plot looks fine now: So, is this a bug of vertcross?

In case you want to reproduce the steps, I attach the data by google drive (download link).

Thanks! Xin

[michaelavs]

Hi @zxdawn, Sorry for not being as active on this issue. I recently went over it with another colleague and here's what we determined: Both projections are correct for what data is being read in and how it is being read in. The reason vertcross shows the arrows going a different direction than interplevel is because vertcross is taking the u and v values as just magnitude while interplevel is taking them as the resultant vector of u and v.

The way this was verified is from looking at the print outs of the values of cross_u and cross_v. When you called cross_u + cross_v Python interpreted that as just adding each integer in the array and then plotting the magnitude and, technically, direction of the arrow (based on the sign of the integer). For the very first arrow at the level of interest, we expect to see an arrow that would have a size and direction proportional to (-4.5392895) + (1.0334523) which is going to be about 3 m/s in the negative direction. When you look at vertcross, the arrow is indeed pointing in the negative x direction and with a length proportional to, what we can assume is, about 3 m/s. When looking at interplevel, you now should be looking at this with respect to the arrows being resultant vectors of the u and v components of the wind. So again looking at the first cross_u and cross_v components, we expect to see a resultant vector that is pointing downward with a slight tilt to the right. We expect this because the u component is much smaller in magnitude than the v component, and the v component being negative means the arrow should point towards the negative y-axis of the plane. As we go along the values for u and v, the arrows should then start to rotate in a clockwise direction as the u component becomes more negative along the line of interest. When looking at interplevel, this is exactly what we see.

I hope this cleared up any confusion you may have with this issue!

[zxdawn]

Thanks, @michaelavs .

When you look at vertcross, the arrow is indeed pointing in the negative x-direction and with a length proportional to, what we can assume is, about 3 m/s.

Yes, the direction of arrow is based on the sign of the integer which is negative. But, from the cross line view [from (118.6, 32.1) to (119, 31.7)], it has to be converted into positive to make the correct plot. So, the cross line direction needs to be considered. Then, users have to decide it manually ... Maybe, it's better to deal with it in wrf-python?

[elliotable]

Hi, I'm one year late but I think I understand what's going on here.. You assume that u_cross and v_cross are projections of u and v respectively on the cross section. But I don't think that vertcross is projecting. It's only evaluating the value of u and v by interpolation on the cross section. You need to do the projection of u_cross and v_cross manually if you want the right angle on your first figure. I hope this is clear !

[zxdawn]

@elliotable thanks for your comment. I understand it now ;) So, users should take care when plotting the cross wind. Is it better to mention it in the User Guide?

[zemega]

Hello. What is the solution for this question? Would it be possible to see an example code to plot circulation vector over a cross section?

[erogluorhan]

@elliotable thanks for your comment. I understand it now ;) So, users should take care when plotting the cross wind. Is it better to mention it in the User Guide?

Glad to see that was helpful! Yes, we can consider clearly mentioning this in the docs. Thanks!

[erogluorhan]

Hello. What is the solution for this question? Would it be possible to see an example code to plot circulation vector over a cross section?

Hi, solution is in the comments made since June 4, 2021, i.e. projection of u_cross and v_cross needs to be taken care of manually by the user, not WRF-Python. Hope that helps.

[zemega]

Hello. This is currently my solution to this question. I hope this solution is correct. Basically I make a new function that rotates the wind vectors from east to west, and south to north in the vertcross output into along the vertcross and perpendicular to vertcross. I would like to the function to be checked whether its correct or not. Although personally, it seems to be correct to me in my project analysis.

from math import acos, degrees, pi, cos, sin
import numpy as np
from netCDF4 import Dataset
from wrf import (getvar, to_np, vertcross,CoordPair)
import xarray as xr
from geopy.distance import geodesic

xr.set_options(keep_attrs=True)

def rotate_ua_va_vert_cross(ua_vertcross, va_vertcross):
    '''
    Takes u and v wind component vertcross, rotates them to align with
    transect meteorological direction, from start_point to end_point. 
    '''

    coord_pairs_1 = to_np(ua_vertcross.coords["xy_loc"])
    coord_pairs_2 = to_np(va_vertcross.coords["xy_loc"])
    if (any(coord_pairs_1 != coord_pairs_2)):
        print("u-component and v component does not match")
        return
    coord_pairs = coord_pairs_1
    main_lat = [(x.lat) for x in coord_pairs]
    main_lon = [(x.lon) for x in coord_pairs]
    # Create an emptry transect
    met_dir_transect = []

    point_a = main_lat[0], main_lon[0]
    point_b = main_lat[1], main_lon[1]
    point_c = main_lat[0], main_lon[1]
    A = geodesic(point_a, point_b).km
    B = geodesic(point_b, point_c).km
    C = geodesic(point_c, point_a).km
    degrees_A_C = 90 + degrees(acos((A * A + C * C - B * B)/(2.0 * A * C)))
    met_dir_transect.append(degrees_A_C)

    for point in range(len(main_lat)):
        if point == 0:
            continue
        point_a = main_lat[point-1], main_lon[point-1]
        point_b = main_lat[point], main_lon[point]
        point_c = main_lat[point-1], main_lon[point]
        A = geodesic(point_a, point_b).km
        B = geodesic(point_b, point_c).km
        C = geodesic(point_c, point_a).km
        degrees_A_B = 180 - \
            degrees(acos((A * A + B * B - C * C)/(2.0 * A * B)))
        met_dir_transect.append(degrees_A_B)

    met_dir_transect_2 = np.array(met_dir_transect, ndmin=1)
    # if deg == True:
    a = met_dir_transect_2/180
    # else:
    #    a = met_dir_transect_2/pi
    c = [cos(pi*X) for X in a]
    s = [sin(pi*X) for X in a]
    c_tile = np.tile(c, (len(ua_vertcross.vertical), 1))
    s_tile = np.tile(s, (len(ua_vertcross.vertical), 1))
    # if clockwise == True:
    un = ua_vertcross * c_tile - va_vertcross * s_tile
    vn = ua_vertcross * s_tile + va_vertcross * c_tile
    # else:
    #    un = ua_cross * c_tile + va_cross * s_tile
    #    vn = va_cross * c_tile - ua_cross * s_tile
    return (vn, un)

wrf_file = Dataset(filename)
ua = getvar(wrf_file, 'ua')
va = getvar(wrf_file, 'va')
p = getvar(wrf_file, "pressure")

start_point = CoordPair(lat=30, lon=80)  # Random
end_point = CoordPair(lat=40, lon=90)  # Random

ua_c = vertcross(ua, p, wrfin=wrf_file, start_point=start_point,
                    end_point=end_point, latlon=True, meta=True,
                    autolevels=180)
va_c = vertcross(va, p, wrfin=wrf_file, start_point=start_point,
                    end_point=end_point, latlon=True, meta=True,
                    autolevels=180)

# horizantal wind and perpendicular wind along the transect line.
h_wind_c, p_wind_c = rotate_ua_va_vert_cross(ua_c, va_c)

[ankanGITHUB22]

Dear all, I also want to plot vertical circulation along the cross section for my study domain. I used the function 'rotate_ua_va_vert_cross' provided by @zemega. But, I am getting weird kind of vertical vectors. I tried without using this function too. But, every time I got weird vertical vectors. Seems like the vectors are coming only at particular levels, not in the vertical direction. I don't know why. For your convenience, I am attaching my code and the output image here. '''

import numpy as np import xarray as xr from matplotlib import pyplot import matplotlib.pyplot as plt from matplotlib.cm import get_cmap from matplotlib.colors import from_levels_and_colors from cartopy import crs from cartopy.feature import NaturalEarthFeature, COLORS from netCDF4 import Dataset from math import acos, degrees, pi, cos, sin from geopy.distance import geodesic from wrf import (getvar, to_np, get_cartopy, latlon_coords, vertcross, cartopy_xlim, cartopy_ylim, interpline, CoordPair) path ='/data/wrfout_d02_2018-02-29_00_00_00' wrf_file = Dataset(path) cross_start = CoordPair(lat=16, lon=78) cross_end = CoordPair(lat=28, lon=92) ter = getvar(wrf_file, "ter", timeidx=73) ht = getvar(wrf_file, "z", timeidx=73) w = getvar(wrf_file, "wa", timeidx=73) u = getvar(wrf_file, "ua", timeidx=73) w_cross = vertcross(w, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True) u_cross = vertcross(u, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True) v = getvar(wrf_file2, "va", timeidx=73) v_cross = vertcross(v, ht, wrfin=wrf_file, start_point=cross_start, end_point=cross_end, latlon=True, meta=True) xr.set_options(keep_attrs=True) def rotate_ua_va_vert_cross(ua_vertcross, va_vertcross): ''' Takes u and v wind component vertcross, rotates them to align with transect meteorological direction, from start_point to end_point. '''

coord_pairs_1 = to_np(ua_vertcross.coords["xy_loc"])
coord_pairs_2 = to_np(va_vertcross.coords["xy_loc"])
if (any(coord_pairs_1 != coord_pairs_2)):
    print("u-component and v component does not match")
    return
coord_pairs = coord_pairs_1
main_lat = [(x.lat) for x in coord_pairs]
main_lon = [(x.lon) for x in coord_pairs]
# Create an emptry transect
met_dir_transect = []

point_a = main_lat[0], main_lon[0]
point_b = main_lat[1], main_lon[1]
point_c = main_lat[0], main_lon[1]
A = geodesic(point_a, point_b).km
B = geodesic(point_b, point_c).km
C = geodesic(point_c, point_a).km
degrees_A_C = 90 + degrees(acos((A * A + C * C - B * B)/(2.0 * A * C)))
met_dir_transect.append(degrees_A_C)

for point in range(len(main_lat)):
    if point == 0:
        continue
    point_a = main_lat[point-1], main_lon[point-1]
    point_b = main_lat[point], main_lon[point]
    point_c = main_lat[point-1], main_lon[point]
    A = geodesic(point_a, point_b).km
    B = geodesic(point_b, point_c).km
    C = geodesic(point_c, point_a).km
    degrees_A_B = 180 - \
        degrees(acos((A * A + B * B - C * C)/(2.0 * A * B)))
    met_dir_transect.append(degrees_A_B)

met_dir_transect_2 = np.array(met_dir_transect, ndmin=1)
# if deg == True:
a = met_dir_transect_2/180
# else:
#    a = met_dir_transect_2/pi
c = [cos(pi*X) for X in a]
s = [sin(pi*X) for X in a]
c_tile = np.tile(c, (len(ua_vertcross.vertical), 1))
s_tile = np.tile(s, (len(ua_vertcross.vertical), 1))
# if clockwise == True:
un = ua_vertcross * c_tile - va_vertcross * s_tile
vn = ua_vertcross * s_tile + va_vertcross * c_tile
# else:
#    un = ua_cross * c_tile + va_cross * s_tile
#    vn = va_cross * c_tile - ua_cross * s_tile
return (vn, un)

h_wind_c, p_wind_c = rotate_ua_va_vert_cross(u_cross, v_cross)

To remove the slight gap between the wind contours and terrain due to the

contouring of gridded data, a new vertical grid spacing, and model grid

taggering, fill in the lower grid cells with the first non-missing value

or each column.

Make a copy of the cross section data. Let's use regular numpy arrays for this.

w_cross_filled = np.ma.copy(to_np(w_cross))

For each cross section column, find the first index with non-missing

values and copy these to the missing elements below.

for i in range(w_cross_filled.shape[-1]): column_vals = w_cross_filled[:,i]

Let's find the lowest index that isn't filled. The nonzero function

# finds all unmasked values greater than 0. Since 0 is a valid value
# for w, let's change that threshold to be -200 instead.
first_idx = int(np.transpose((column_vals > -200).nonzero())[0])
w_cross_filled[0:first_idx, i] = w_cross_filled[first_idx, i]

Get the terrain heights along the cross section line

ter_line = interpline(ter, wrfin=wrf_file, start_point=cross_start, end_point=cross_end) lats, lons = latlon_coords(w)

Get the cartopy projection object

cart_proj = get_cartopy(w)

Create the figure

fig = pyplot.figure(figsize=(8,6)) ax_cross = pyplot.axes()

Make the cross section plot for w

w_levels = np.arange(-4E-1, +4E-1, 5E-2) xs = np.arange(0, w_cross.shape[-1], 1) ys = to_np(w_cross.coords["vertical"]) w_contours = ax_cross.contourf(xs, ys, to_np(w_cross_filled), levels=w_levels, cmap='seismic', extend="both")

Add the color bar

cbar = fig.colorbar(w_contours, ax=ax_cross) cbar.ax.tick_params(labelsize=12) cbar.set_label('Vertical Velocity (m.s-1)', rotation=-270, fontsize=12)

Fill in the mountain area

ht_fill = ax_cross.fill_between(xs, 0, to_np(ter_line), facecolor="black")

quiver

ax_cross.quiver(xs[::5], ys[::5], to_np(h_wind_c[::5, ::5]+p_wind_c[::5, ::5]), to_np(w_cross[::5, ::5]*100))

ax_cross.quiver(xs[::5], ys[::5],to_np(u_cross[::5, ::5]+v_cross[::5, ::5]), to_np(w_cross[::5, ::5]))

Set the x-ticks to use latitude and longitude labels

coord_pairs = to_np(u_cross.coords["xy_loc"]) x_ticks = np.arange(coord_pairs.shape[0]) x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]

Set the desired number of x ticks below

num_ticks = 5 thin = int((len(x_ticks) / num_ticks) + .5) ax_cross.set_xticks(x_ticks[::thin]) ax_cross.set_xticklabels(x_labels[::thin], rotation=90, fontsize=8)

Set the x-axis and y-axis labels

ax_cross.set_xlabel("Latitude, Longitude", fontsize=12) ax_cross.set_ylim([0,10000])

pyplot.show()

Did I do any mistake in my code? Can anyone please have a look at these and help me to solve this issue? That will be very helpful for me. I am really looking forward to hearing from you. Thank you for your time and consideration in advance. With regards, Ankan