Logo

Archive for the ‘Weather’ Category

Wind Observation Calculations in FORTRAN and Python

The National Climatic Data Center (NCDC) is the government organization in charge of preserving and archiving weather observations in the United States. One of the most popular datasets they provide is the Local Climatic Data (LCD) product. It contains daily summaries and raw hourly reports from about 1,600 stations. The daily summary includes high temperature, low temperature, amount of precipitation, atmospheric pressure, and wind summaries. An example can be found here.

The daily summary has several wind observations. There are maximum five-second and two-minute speeds, and the average wind direction during those gusts. And then there are three aggregate observations: “Resultant wind speed”, “Resultant wind direction”, and “Average wind speed”. What do these observations mean, and what is the difference between “resultant” and “average”?

The average wind speed is a simple average of wind velocity, regardless of direction, over the day. If the wind blew from the north at 10 knots for half the day, and from the south at 20 knots for the other half, the average wind speed that day would be 15 knots.

Resultant wind speed and resultant wind direction are the vector average (the resultant vector) wind speed and direction. Imagine the wind observation as a vector starting from the origin (the observing station) pointing in the direction the wind is originating from and with a length equal to the wind speed. If you take all those vectors and add them up, the resulting direction of the vector would be the resultant wind direction, and the magnitude of the resulting vector divided by the number of vectors is the resultant wind speed.

Adding all of the hourly wind vectors and averaging them results in the values expressed in the daily resultant wind speed and resultant wind direction observations. That is how the NCDC calculates those values, and they do it in FORTRAN. I am going to walk through the specific code that generates those daily averages from hourly data. I’ll also show the same code in Python. Finally, using the complex math library in the Python standard library, I’ll show how the vector average can be calculated with far less code. The full FORTRAN example can be downloaded here, while the Python source can be found here.

The following FORTRAN code was supplied by the NCDC, and is the actual code they use to compute daily resultant wind direction and speed. The Python code was created by me to exactly match the FORTRAN code, to compare and contrast the differences between FORTRAN and Python.

1
2
3
4
5
6
7
8
9
10
11
12
13
      SUBROUTINE R_WIND (N_RWD, RRDD, RWND, RWDIR, RWSPD)
 
      IMPLICIT NONE
 
      REAL*8 RRDD(31), RWND(31), R_DIR, R_SPD, TOTRX, TOTRY
 
      INTEGER X, WI1, WI2, N_RWD
 
      CHARACTER*2 RWDIR
      CHARACTER*5 RWSPD, WC1
 
      TOTRX = 0.0
      TOTRY = 0.0

The subroutine R_WIND takes three parameters, and returns two. N_RWD is the number of hourly wind observations. RRDD is the array of wind directions, while RWND is the array of wind speeds. The wind direction angles passed in are degrees plus 180 converted into radians.

The subroutine returns two results, RWDIR, which is the resultant wind direction, and RWSPD, which is the resultant wind speed. These values are returned as character strings.

The Python equivalent is much shorter in code, because in Python you don’t need to explicitly declare variable types. You do have to import the math library, since the sine and cosine functions that will be used are not implicitly available, as with FORTRAN. There are also no explicitly declared output variables. Because of this dynamic nature, documentation is even more important. The text between triple-quotes immediately following the function definition is called a “docstring.” A docstring is used in Python to document the code, and is used by documentation generators like Sphinx to build web-based documentation for your code.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import math
 
def r_wind(n_rwd, rrdd, rwnd):
    """
    Takes an array of houry wind directions and speeds and returns the
    resultant average wind speed and direction, or None for speed and direction
    if there are fewer than 16 hourly observations.
 
    n_rwd -- number of hourly wind observations
    rrdd  -- array of wind directions plus 180 degrees, in radians
    rwnd  -- array of wind speeds
 
    Returns a tuple of resultant wind direction and wind speed.  Both are
    returned as strings.  Wind direction is in tens of degrees.
    """"
    rwdir = rwspd = None
    totrx = 0.0
    totry = 0.0

Once the variables have been declared, the subroutine then checks to see if there are at least 16 hourly observations. If there are not at least 16 hourly observations out of 24 (or two-thirds), it won’t try to calculate a daily resultant wind vector. Then, it loops through the direction and speed arrays, converting into Cartesian coordinates. Converting polar (angle and magnitude, in this case direction and speed) coordinates into Cartesian coordinates to sum them up, then average, is the normal way to add vectors.

15
16
17
18
19
20
21
      IF (N_RWD .GT. 15) THEN
         X = 1
         DO WHILE (X .LE. N_RWD)
            TOTRX = (RWND(X)*(DSIN(RRDD(X)))) + TOTRX  ! TOT X & Y COOR.
            TOTRY = (RWND(X)*(DCOS(RRDD(X)))) + TOTRY
            X = X + 1
         END DO

There are a couple of unique things happening in this code. First, adding 180 degrees to the wind directions prior to passing in the array. Second, using sine for the X coordinate conversion and cosine for the Y coordinate conversion is the opposite of the general way to convert from polar to Cartesian coordinates (see here). The reason this is done, even though it isn’t required in this subroutine, is to create a Cartesian plane that matches a map. When wind is measured at 0 degrees, it means it is coming out of the north. The visualized vector would be an arrow pointing down (the negative y-axis in a standard Cartesian plane). Additionally, the angle would be measured clockwise (a 90 degree wind would be coming from the east, and thus a vector pointing left along the negative x-axis). Using the standard Cartesian coordinate conversion, the zero angle vector would point right along the positive x-axis, with angles measured counter-clockwise.

20
21
22
23
24
25
    if n_rwd > 15:
        x = 0
        while x < n_rwd:
            totrx = (rwnd[x]*(math.sin(rrdd[x]))) + totrx # tot x & y coor.
            totry = (rwnd[x]*(math.cos(rrdd[x]))) + totry
            x = x + 1

The Python version is quite similar. The biggest difference is that Python arrays are zero indexed, while FORTRAN has arrays indexed starting at one. Another difference is that Python arrays use brackets to indicate the index, with parentheses reserved for function calls. Also, no “end while” statement is required, because in Python, indentation has syntactical meaning. Specficially, anything indented past the while statement are considered part of the while block. Once code returns to the same indentation as the while, it is not longer part of the while block. Finally, other differences are that in Python, parentheses are not required around conditional expressions, and comments in Python are indicated with a hash character rather than a bang.

The next section of the FORTRAN code finds the angle of the resultant x,y coordinates and averages the resultant length (which is speed).

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
         R_DIR = DATAN(TOTRX/TOTRY) ! AVG R_WIND DIR
         R_DIR = R_DIR / (3.14156/180)   ! convert radians back to degrees
 
         IF ((TOTRX*TOTRY) .LT. 0) THEN
            IF (TOTRX .LT. 0) THEN
               R_DIR = 180 + R_DIR
            ELSE
               R_DIR = 360 + R_DIR
            END IF
         ELSE
            IF (TOTRX .GT. 0) R_DIR = 180 + R_DIR
         END IF
 
         R_SPD = (((TOTRX**2) + (TOTRY**2)) / N_RWD**2)**.5
 
         IF (R_DIR .EQ. 0 .AND. R_SPD .NE. 0) R_DIR = 360
         IF (R_SPD .EQ. 0) R_DIR = 0

Lines 23 and 24 take the arctangent of the x,y coordinate to compute the angle to convert back into polar coordinates. Line 27-35 convert the angle to the original solution space, based on the quadrant the resultant x,y coordinate resides in. Line 37 is simply the distance between the origin and the averaged x,y coordinate, which is the vector length of the averaged resultant vector, which is resultant wind speed. This is calculated by realizing that x and y are just legs of a right triangle, with hypotenuse being the distance between the origin and the point. Line 39 and 40 are just clean up. The subroutine will return an angle of 360 for a north wind, and zero only if there is resultant wind. Line 40 is required because the DATAN function will return NaN (Not A Number) if x and y both sum to zero, which is the only time the speed would be zero.

The Python version is quite similar. The only major difference is that there needs to be a check that totry isn’t zero, as that would result in a ZeroDivisionError Python exception. When a division by zero is encountered in Python, it will interrupt code flow and throw an exception. FORTRAN on the other hand returns the discrete value NaN which is then passed to the DATAN function, which knows how to handle NaN. There actually is a NaN in Python, but Python chooses to throw an exception rather than return NaN on a division by zero.

35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
        if totry == 0.0:
            r_dir = 0
        else:
            r_dir = math.atan(totrx/totry) # avg r_wind dir
        r_dir = r_dir / (3.14156/180) # convert radians back to degrees
 
        if (totrx*totry) < 0:
            if totrx < 0:
                r_dir = 180 + r_dir
            else:
                r_dir = 360 + r_dir
        else:
            if totrx > 0:
                r_dir = 180 + r_dir
 
        r_spd = ((totrx**2 + totry**2) / n_rwd**2)**.5
 
        if r_dir == 0 and r_spd != 0:
            r_dir = 360
 
        if r_spd == 0:
            r_dir = 0

The final section of the subroutine just rounds and formats. It rounds the degrees to tens of degrees, and returns in tens of degrees. Thus, 260 degrees would be returned as 26. Speed is returned as tenths of a knot. The values are returned as strings, as they will be ending up in a CSV file, not used in further processing.

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
         WRITE (WC1, FMT='(F4.0)') R_DIR
 
         READ (WC1, 10) WI1, WI2
 10      FORMAT (I2, I1)
 
         IF (WI2 .GT. 4) THEN
            WI1 = WI1 + 1
         END IF
 
         WRITE (RWDIR, FMT='(I2)') WI1
         WRITE (RWSPD, FMT='(F5.1)') R_SPD
      END IF
 
      RETURN
      END

The Python version is noticeably different as formatting and substring extraction is quite different between FORTRAN and Python. Python uses C-style format strings. The FORTRAN code is also using CHARACTER arrays rather than STRING arrays (perhaps when this was written the FORTRAN dialect used did not support STRING), which doesn’t support substring notation. Python strings are character arrays and provide all the same functionality as other arrays, including sub-array indexing, which is used here. Finally, in Python don’t declare the return variables in the function declaration, you return them explicitly, and anonymously, which gives the programmer a lot of power, but also adds to the responsibility of the programmer to document what is returned in docstrings.

58
59
60
61
62
63
64
65
66
67
68
        wc1 = '%#4.0f' % r_dir
        wi1 = int(wc1[0:2])
        wi2 = int(wc1[2:3])
 
        if wi2 > 4:
            wi1 = wi1 + 1
 
        rwdir = '%2i' % wi1
        rwspd = '%#5.1f' % r_spd
 
    return rwdir, rwspd

That’s it. Of course, if we were to rewrite this function today in Python, we wouldn’t need to do all of this calculation in our code, but we can instead rely on higher level functions in the standard library to help us out. This is probably also true of FORTRAN today, but I’ve not actively written FORTRAN in many years, so don’t know.

Python has a complex math library, called cmath that we can import in order to give us the functionality we need. We can convert polar coordinates into a Cartesian complex number, add up all the complex numbers and average them, then convert the complex number back to polar coordinates. This will accomplish the task of the original function, which is to return the resultant wind vector angle and speed.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def r_wind_new(n_rwd, rrdd, rwnd):
    wind_vector_sum = None
 
    if n_rwd > 15:
        for i in range(0, n_rwd):
            wind_polar = cmath.rect(rwnd[i], rrdd[i] - math.pi)
            if wind_vector_sum is None:
                wind_vector_sum = wind_polar
            else:
                wind_vector_sum += wind_polar
 
        r, phi = cmath.polar(wind_vector_sum / n_rwd)
 
        rwdir = int(round(int(round(math.degrees(phi) % 360)) / 10.0))
        rwspd = int(round(r*10))/10.0
    else:
        return None, None
 
    return '%2i' % rwdir, '%#5.1f' % rwspd

This new Python function, r_wind_new is a drop-in replacement for the original r_wind function. It uses cmath.rect() to convert a vector angle and length to a Cartesian complex number. These complex numbers are added together and then averaged. Then, cmath.polar() is used to convert the complex number into an angle and length, which is the resultant wind vector.