How to calculate wind direction from U and V wind components in R

Question

I have U and V wind component data and I would like to calculate wind direction from these values in R.

I would like to end up with wind direction data on a scale of 0-360 degrees, with 0° or 360° indicating a wind blowing to the north, 90° indicating a wind blowing to the east, 180° indicating a wind blowing to the south and 270° indicating a wind blowing to the west.

Below is some example data:

wind <- read.table(header = TRUE, text = "
u_ms     v_ms     u_rad        v_rad
-3.711   -1.471   -0.064769155 -0.025673788
-2.2417  -1.6118  -0.039125038 -0.028131211
-1.8188  -1.6613  -0.031744042 -0.028995149
-1.6164  -1.7037  -0.028211496 -0.029735168
-1.3941  -1.7388  -0.02433163  -0.030347779
-1.0682  -1.8748  -0.018643603 -0.032721426
-0.57611 -1.8359  -0.010055014 -0.032042493
-1.5698  -1.6766  -0.027398173 -0.029262184
-1.4976  -1.6994  -0.026138045 -0.029660119
-1.3537  -1.7505  -0.023626517 -0.030551982
-1.0901  -1.4947  -0.01902583  -0.026087431
-0.60403 -0.96283 -0.01054231  -0.01680455
-0.70812 -1.1194  -0.012359023 -0.019537212
-0.49045 -0.6849  -0.008559966 -0.011953758
-0.39849 -0.7847  -0.006954961 -0.013695596
0.17875  -0.80349 0.003119775  -0.014023543
0.48356  -0.19352 0.008439712  -0.00337756
1.5082   -0.97815 0.02632305   -0.017071935
1.4219   -1.0835  0.024816831  -0.018910639
2.5881   -0.81666 0.045170857  -0.014253403
")

I have used the following code to try and obtain wind direction (column td), but I am not convinced that the returned angles are those that I want (i.e. 0°/360° indicating a wind blowing to the north, 90° indicating a wind blowing to the east etc…).

u = wind$u_rad # u component in radians
v = wind$v_rad # v component in radians

d = (180/pi)*(atan2(u,v))
td = as.matrix(d + 180)
df = cbind(wind, d, td)

> df
       u_ms     v_ms        u_rad       v_rad         d        td
1  -3.71100 -1.47100 -0.064769155 -0.02567379 -111.6228  68.37716
2  -2.24170 -1.61180 -0.039125038 -0.02813121 -125.7164  54.28357
3  -1.81880 -1.66130 -0.031744042 -0.02899515 -132.4087  47.59129
4  -1.61640 -1.70370 -0.028211496 -0.02973517 -136.5062  43.49379
5  -1.39410 -1.73880 -0.024331630 -0.03034778 -141.2788  38.72124
6  -1.06820 -1.87480 -0.018643603 -0.03272143 -150.3269  29.67308
7  -0.57611 -1.83590 -0.010055014 -0.03204249 -162.5780  17.42199
8  -1.56980 -1.67660 -0.027398173 -0.02926218 -136.8842  43.11576
9  -1.49760 -1.69940 -0.026138045 -0.02966012 -138.6118  41.38819
10 -1.35370 -1.75050 -0.023626517 -0.03055198 -142.2844  37.71557
11 -1.09010 -1.49470 -0.019025830 -0.02608743 -143.8963  36.10365
12 -0.60403 -0.96283 -0.010542310 -0.01680455 -147.8980  32.10204
13 -0.70812 -1.11940 -0.012359023 -0.01953721 -147.6830  32.31699
14 -0.49045 -0.68490 -0.008559966 -0.01195376 -144.3939  35.60607
15 -0.39849 -0.78470 -0.006954961 -0.01369560 -153.0774  26.92258
16  0.17875 -0.80349  0.003119775 -0.01402354  167.4578 347.45783
17  0.48356 -0.19352  0.008439712 -0.00337756  111.8112 291.81121
18  1.50820 -0.97815  0.026323050 -0.01707193  122.9656 302.96561
19  1.42190 -1.08350  0.024816831 -0.01891064  127.3077 307.30771
20  2.58810 -0.81666  0.045170857 -0.01425340  107.5128 287.51279

I would appreciate any advice on whether my method is correct, and if not how I could correctly obtain the desired wind direction values. While Calculating wind direction from U and V components of the wind using lapply or ifelse was helpful, the code did work with my data, and I am sure that there is an easier away to obtain wind direction. Many thanks!

Insurgent

Answer 1

There are three problems with this:

1. You cannot convert m/s to radians. In order to input wind components into atan2, you must normalize them, but you don't do this by multiplying m/s by pi/180 (which you did to get u_rad and v_rad). You should make a column of absolute windspeed (sqrt(u_ms^2 + v_ms^2)) and take atan2(u_ms/wind_abs, v_ms/wind_abs). (also note that atan2 takes y component first - make sure that's what you want)
2. atan2 will give you an answer in the unit circle coordinates, which increase counterclockwise and have a zero on the x-axis. You want an answer in cardinal coordinates which increase clockwise and have a zero on the y-axis. To convert unit circle to cardinal coordinates, you must subtract the unit circle angle from 90.
3. You must know whether the wind info refers to the direction the wind is coming from (standard for cardinal coordinates) or the direction the wind is blowing (standard for trig/vector operations)

If you are given u_ms = = -3.711 and v_ms = -1.471 (on the unit circle it is blowing down and slightly to the left, so it is coming from the northeast), then:

wind_abs = sqrt(u_ms^2 + v_ms^2)
wind_dir_trig_to = atan2(u_ms/wind_abs, v_ms/wind_abs) 
wind_dir_trig_to_degrees = wind_dir_trig_to * 180/pi ## -111.6 degrees

Then you must convert this wind vector to the meteorological convention of the direction the wind is coming from:

wind_dir_trig_from_degrees = wind_dir_trig_to_degrees + 180 ## 68.38 degrees

Then you must convert that angle from "trig" coordinates to cardinal coordinates:

wind_dir_cardinal = 90 - wind_dir_trig_from_degrees
[1] 21.62284 #From the northeast.

Iranian

@ Senor O: Thank you for your explanation. It is very helpful. My question is, why do you need the final step? I have discovered that using (270-(atan2(v_ms, u_ms)*(180/pi)))%%360 gives the same result as the first two steps above (i.e. 68.38). However, does the 270 not already convert from the unit circle into meteorological wind direction? (see eol.ucar.edu/content/wind-direction-quick-reference). I am also looking for the direction that the wind is blowing towards (td) not from (fd). This can presumably be obtaining using `ifelse(fd>180, fd-180, (fd-180)+360)? Thanks. – Insurgent

Doing 270 - atan... is sufficient - I just broke it into two steps to explain what's going on (add 180 to change from to 'to', and subtract angle from 90 to change trig to cardinal). If you are looking for the direction the wind is blowing to, then skip my step where you add 180, Or replace 270-atan2... with 90-atan2... – Robles

Point 1. is certainly nonsense; atan2(x,y)=atan2(ax,ay) for any positive a by definition (numerical errors may apply for huge values, but this is rather unplausible in this use case). Also u and v may not necessarily be on EW NS directions; for instance many numerical models will output them on their internal grid (and they cannot be easily converted without knowing how precisely this grid is defined). – Butterscotch

In your example -u is larger than -v, so your wind direction should be between 45 and 90 degrees (if -u becomes very large compared to -v we should have a 90 degrees wind direction). I think in your example 68.38 is the right answer: I think the last step "90 - wind_dir_trig_from_degrees" should not happen – Anthe

Answer 2

While the accepted answer has the right idea it has a flaw. As mentioned in the comments one does not need to normalize the u and v component in order to use atan2 on them. The flaw comes when u == v == 0 and wind_abs becomes 0. In C# the two divisions will return infinity (in compliance to IEEE 754) and atan2 will return NaN. When not normalizing the components atan2(0,0) happily returns 0. So not only is normalizing not necessary, it also introduces an error.

Please also be aware that the most common function signature of atan2 is atan2(y, x) -- Microsoft Excel being an exception.

Encage

Answer 3

In Python:

Dir=np.mod(180+np.rad2deg(np.arctan2(U, V)),360)

This leads to a Southerly wind (180 degrees) for a [0,1] vector, a Northerly wind (0 degrees) for a [0,-1] vector, a South Westerly wind (225 degrees) for a [1,1] vector:

U
Out[86]: 
array([[ 0.  ],
       [ 0.  ],
       [ 1.  ],
       [-3.47]])

V
Out[87]: 
array([[ 1.  ],
       [-1.  ],
       [ 1.  ],
       [-1.47]])

np.mod(180+np.rad2deg(np.arctan2(U, V)),360)
Out[89]: 
array([[180.        ],
       [  0.        ],
       [225.        ],
       [ 67.04097233]])

Anthe