Channels ▼
RSS

Database

Rotating a Weather Map

Source Code Accompanies This Article. Download It Now.


Jun99: Rotating a Weather Map

Robert is a Staff Member of the Air Traffic Surveillance group at MIT's Lincoln Laboratory. He can be contacted at grappel@ ll.mit.edu.


I was recently part of a project developing a system for aircraft pilots to access the national ground weather-radar database while in flight. This weather-radar graphical database is generated from the outputs of the FAA and National Weather Service network of radars covering the continental United States and is updated every five minutes. Each pixel in the database covers a square measuring two kilometers (about one nautical mile) on a side. The content of each data pixel is a measure of the radar reflectivity measured at that location -- radar reflectivity is proportional to the water content in the atmosphere (the precipitation rate).

This graphical database is available through several commercial vendors -- it's what you see displayed on The Weather Channel or during typical TV weather reports. Our system, on the other hand, provides a low-speed digital datalink connection from an FAA ground computer to an avionics computer/display located in the aircraft cockpit. The pilot can request the uplink of a portion of the weather database centered on a specified location (the aircraft's current position, a particular airport, and so on) and with a range of up to 200 nautical miles from the centerpoint. (The pilot can also request images from past databases to observe storm motion at a particular location.) The actual graphical data uplinked to the aircraft for a given map image consists of an array of 256×256 two-bit pixels, compressed to about 3500 total bits using a proprietary, lossy technique developed for the FAA. The aircraft computer/display avionics (effectively a 25-MHz, 486-based embedded PC running DOS) decompresses the uplinked image and displays it with the weather intensities color-coded to parallel an airborne weather radar (light precipitation in green, medium precipitation in yellow, and heavy precipitation in red). Figure 1 shows the instrument panel of the test aircraft (Cessna 172 Skyhawk). The datalink display/keyboard is the ARNAV MFD (multifunction display) 5100 CRT in the radio stack (just right of center). The ARNAV MFD 5010 avionics computer is located at the far right, just behind the right control wheel. (In a normal installation, the MFD 5010 would be installed out of sight.) The datalink "modem" function is performed by the Bendix/King KT70 transponder just below the MFD 5010. The weather display shown here is an actual North-up weather image center on the Atlanta, Georgia airport. The weather image is showing a 50 nautical mile radius, with some severe storm cells in the area.

Figure 2, on the other hand, is a close- up of the datalink demonstration unit which groups the three datalink components together. The weather display is the same as in Figure 1. The soft-key labels on the right side of the display indicate the Traffic Information Service (TIS) and Weather Request (WXREQ) pilot inputs. The yellow TIS ALERT indication tells the pilot that there is another aircraft nearby that may be on a conflicting course. (The pilot presses the TIS button to get a display showing where other aircraft are in relation to the pilot's aircraft.) Note the PCMCIA slot on the ARNAV MFD 5010.

The weather map graphical image to be uplinked in response to a pilot's request is simply windowed from the selected national database. The image is oriented North-up, as is conventional for maps. This map-like display is desirable when pilots specify an airport or other landmark as the display centerpoint, but it is undesirable when pilots want the map centered on the aircraft's current position. Pilots would like the map display to match what is visible out of the aircraft windshield -- rotated so that the aircraft's current heading points to the top of the display. However, the ground weather computer doesn't know the heading of every aircraft that might make a request, and aircraft can and do maneuver after the request has been made and processed. What we needed was a way for the avionics computer/display to rotate the uplinked North-up weather image to correspond with the aircraft's current heading.

Clearly, performing such a rotation can take a lot of processing. There are 65,536 pixels to be recomputed for each weather map image. We decided that the avionics display would need to be refreshed about once per second -- and the relatively slow airborne computer is likely to be busy doing other tasks (graphical weather display is only one of its functions). Performing a map image rotation would have to take only a fraction of a second to be practical. All the software had to be written in standard C and be as processor/operating system independent as possible. In this article, I'll describe the algorithm we developed to efficiently perform this rotation of graphical weather maps. I'll also suggest some techniques and approaches that you could use to optimize other time-limited computer applications.

Warning...Trigonometry Ahead!

The first step in developing the map rotation algorithm is to convert each map pixel's Cartesian row index (y-coordinate) and column index (x-coordinate) into polar coordinates. In polar coordinates, the location of each pixel is defined by its distance (R) from the center point of the map and an angle (A) around the map center point. Exactly how the angle A is to be measured is determined by convention. Math books define the polar angle as measured counterclockwise from the positive x-axis (due East). Map makers, however, measure the angle clockwise from due North (the positive y-axis). Since we're doing map rotation here, I'll use the map-maker's convention for the polar angle A. Hence, the Cartesian to polar conversion equations are:

X=Rsin(A)

Y=Rcos(A)

Rotating the weather map around its center point by the angle B leaves R unchanged, while the angle A changes to A+B. The equations for the rotated pixel row index (Yrot) and the rotated pixel column index (Xrot) are simply:

Xrot=Rsin(A+B)

Yrot=Rcos(A+B)

Applying the standard formulas for the sine and cosine function of the sum of two angles and a bit of algebra yields:

Xrot={Rsin(A)cos(B)}+{Rcos(A)sin(B)}

Yrot={Rcos(A)cos(B)}-{Rsin(A)sin(B)}

At first, these rotation equations don't look too promising. Performing all that math for each pixel in the rotated map will take a lot of processing time. We'll have to simplify this before having an efficient implementation.

Simplify...Simplify...Simplify!

The first simplification we can perform on the rotation equations is to recognize that the rotation angle B is a constant for each pixel in the map and needs to be calculated only once for the entire map. Let S=sin(B) and C=cos(B). The second simplification is to notice that the rotation equations are actually operating on the initial, unrotated, Cartesian X and Y coordinates. Hence, we can rewrite the rotation equations for each pixel as:

Xrot=(XC)+(YS)

Yrot=(YC)-(XS)

At this point, we're down to four multiplications, one addition, and one subtraction per map pixel. However, we can do better. Let's attack the multiplications, since there are more of these than any other operation (and multiplication tends to be a more time-consuming operation than addition or subtraction).

The values of X and Y that we will see are just the indices of the map pixels. Since the weather map is square (symmetric in X and Y), we can precompute the results of the multiplications for one edge of the map and store them in two tables. If you assume that the map has N pixels on a side, the precomputation code looks like Figure 3(a) and the rotation equations are reduced to Figure 3(b).

That's looking much better. We've replaced all the multiplications with table lookups. Note, however, that each of the rotation equations contains one table lookup based on row index (Y) and a second table lookup based on column index (X). Since the overall rotation algorithm is going to iterate over all map rows and all map columns, we can simplify the table lookups even further. We don't need to redo the row lookups for each pixel on a given row. Figure 4 presents pseudocode for the basic loop structure of the improved rotation algorithm.

It's Never Quite that Simple!

Unfortunately, things aren't quite as easy as that. First of all, consider the square input map rotated by some angle B that is not a multiple of 90 degrees. The corners of the input map stick out over the edges. There are some input map pixels that will not get mapped to a rotated output pixel. Also, the corners of the rotated map stick out beyond the edges of the input map. You can deal with the first situation by simply initializing the output map ahead of time to a fixed null value. (C's memset function performs this whole-array initialization very efficiently.) In order to deal with the second problem, you need to insert some tests in the basic rotation algorithm to check that the rotated coordinates (Xrot, Yrot) do, in fact, denote valid map pixels. We'll also need a couple of if statements inserted before the last line of the pseudocode:

IF (Xrot < 0) OR (Xrot >= N) CONTINUE IF (Yrot < 0) OR (Xrot >= N) CONTINUE

There's one more problem lurking in the basic pseudocode. The calculations for the rotated pixel coordinates Xrot and Yrot are done in floating-point math, yet the map pixel coordinates must be integers. We need to be careful about rounding these floating-point-to-integer conversions. Normally, C would truncate the floating-point value -- we appear to need rounding code inserted into the rotation algorithm. Adding this code into the algorithm's inner loop will be costly in terms of execution time, but it appears that this can't be helped. Each rotated map pixel (a little square) might cover parts of multiple input pixels -- we need to be careful how we determine which output map pixel we choose for each input map pixel. We could introduce distortions into the rotated map -- or we could generate holes in the rotated map.

Turn the Problem on Its Head

Sometimes, the best way to attack an optimization problem is to look at the algorithm from another direction. When I had just about conceded that rounding code was a necessary evil, a friend suggested that I look at the rotation problem differently. Instead of rotating the input map into the output map through a rotation angle B, try rotating each pixel of the output map back to a pixel in the input map by the angle -B. This has no effect on the derivation of the rotation equations mentioned earlier -- we don't care what the value of B is. Now, the aforementioned last line in the Figure 4 pseudocode reads

output_map[I][J] = input_ map[Xrot][Yrot].

At first, this doesn't appear to be an improvement over the original form. We still have to convert the rotated pixel coordinates from floating point to integers. The point to notice is that now every output map pixel gets set from some input map pixel. We can skip the rounding and every output pixel will still be set to some input pixel value -- no holes. In fact, some input pixels might be used more than once, depending on the rotation angle. We can let C do its normal truncation and the result will be free of holes.

The Final Rotation

Listing One is the C code for the weather map rotation function. As you can see, it consists of little more than the pseudocode I've already discussed. The actual weather maps in our system do not use the 0th row or column, and an additional "border" row and column is added to each map. This results in a border of zero-value pixels around the map, and also impacts the pixel indexing.

Also, the use of C pointer "magic" to perform a minor additional optimization reduces the doubly indexed store into the output map to a singly indexed store.

Conclusion

This algorithm implementation executes a weather map rotation in about one third of a second for a 256×256-pixel weather map. The requirements of my application were met. Also, I got a chance to apply several types of optimization strategies that may help you in your computer applications. Opinions, interpretations, conclusions, and recommendations are those of the author and are not necessarily endorsed by the FAA.

DDJ

Listing One

/* Function to rotate a weather map through an arbitrary angle.
   Arguments: pointer to square input map.
              pointer to square rotated output map.
              size of maps in pixels-per-side.
              rotation angle (in degrees clockwise from North).
   Note: both maps do not use the zero'th row or column.
         both maps sized a maximum of MSIZE * MSIZE */

#include <stdio.h>
#include <math.h>
#include <string.h>

#define MSIZE 258
#define PI 3.141592653

void rotate_map(inmap,outmap,nij,angle) int nij; double angle;
     char inmap[MSIZE][MSIZE],outmap[MSIZE][MSIZE];
{
   char * optr;
   int i,j,nijh,x,y;
   double n0,x0,s,c,cc0,ss0,s0[MSIZE],c0[MSIZE];

   memset(outmap,0,MSIZE*MSIZE);    /* clear output map */ 
   nijh = (nij / 2) + 1;            /* midpoint of the map */
   angle *= (-PI / 180.0);          /* convert neg. angle to radians */
   s = sin(angle);  
   c = cos(angle);
   n0 = 0.5 - nijh;
   for (i=1; i <= nij; i++)
   { /* precompute sine and cosine tables */
       x0 = (double)i + n0;
       c0[i] = x0 * c;  
       s0[i] = x0 * s;
   }
   for (i=1; i <= nij; i++)
   {
       cc0 = nijh + c0[i];  
       ss0 = nijh - s0[i];
       optr = &outmap[i][0];
       for (j=1; j <= nij; j++)
       {
           x = cc0 + s0[j];
           if ((x < 0) || (x > nij)) continue;
           y = ss0 + c0[j];
           if ((y < 0) || (y > nij)) continue;
           optr[j] = inmap[x][y];
       }
   }
}
        




Back to Article


Copyright © 1999, Dr. Dobb's Journal

Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.
 

Video