**Archive for the ‘FORTRAN’ Category**

January 16th, 2011 by efloehr

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.