Radio propagation widgets & utilities

Repository:	http://gitorious.org/splat2
Mailing list:	http://groups.google.com/group/qt-ham

Splat! 

The other day I stumbled upon Splat. Nice, and open-source.

Unfortunately, it's rather slow and cannot effectively be parallelized (e.g. with OpenMP) due to static variable usage. It also lacks a GUI: you must start it with command-line parameters, then you get a graphics file as output. This graphics file happens to be unwildly *.ppm format, huge and slow.

Nevertheless, I got hooked and interested. At first I looked into the itm.cpp file that implements the Longley-Rice propagation model, also known as "Irregular Terrain Model" or ITM.

Longley-Rice / ITM 

Then I surfed the web and found some nice documents about ITM.

First I found at http://flattop.its.bldrdoc.gov/itm.html the description and algorithm of ITM.

I used this to heavily document the source code, thus yielding my own version of itm.cpp (this is a different link than the one above!).

Other notes helped me:

• Technical Note TN101 from ITS, dated back to 1967 and written with a type-writer. Therefore no easy read, especially for the formulas, but still it gives good background.
• An introduction from NTIA about Longley-Rice in area prediction.

Some more are listed in my README file.

I still have not comprehended all of itm.cpp, questionable things I marked with "XXX". For example, sometimes the reference documentation use different values than Splat!, e.g. in Algorithm 6.4. I have yet to find out why. One reason might be that the author of Splat! did all calculations in outdated units (miles, yards, ell bows, whatever --- smile ---) and not in SI-units.

Geographical calculations 

Base of every propagation program are geographical (or geodetic) calculations. Some of this I ripped out of Splat!, some I wrote by myself.

As usual I wrote this in Qt4. Therefore it should be dead simple to re-compile this for Windows, MacOSX ... or Linux, which I prefer.

With GeoUtils.cpp / GeoUtils.h you can do:

Great circle distance 

The great circle distance is the shortest distance between two points on a sphere. The difference between the airports of Los Angeles and Nashville is:

const QPointF BNA(-86.67, 36.12);
const QPointF LAX(-118.4, 33.94);

double d = Great_Circle_Distance(BNA, LAX);
qDebug("Great circle distance Nashville -> Los Angeles: %.2f km\n", d/1e3);

And it will calculate 2886.45, the same value as the linked wikipedia article above yields.

You might notice that I store the longitude first, as this corresponds to an X-coordinate. Then comes the latitude, which is the Y-coordinate. I'm used to x/y coordinates and re-use that scheme in all of my code.

Azimuthal angle 

This is the angle that you have to turn from north clockwise to look to the desired point. The angle between the 2m FM repeater DB0FT and my home QTH can be calculated like this:

    const QPointF DB0FT(8.45741, 50.23153);
    const QPointF DH3HS(8.77485, 50.28425);

    double a = AzimuthR(DB0FT, DH3HS);
    qDebug("Azimuth DB0FT - DH3HS: %.2f degrees\n", a / DEG2RAD);

where DEG2RAD is TWO_PI (from math.h> divided by 360. This calculation yields 75.32 degrees.

Make path 

If you have two locations, know the azimuthal angle and how to traverse a great circle, you can make the path between those locations. For example, the path from my home QTH to the relay, in 500 meter steps, can be calculated like this:

    GeoData data;
    MakePath(DH3HS, DB0FT, data, 500);

GeoData will receive the result. It basically consists of bunch QVector member variables, which will get the longitude, latitude and distance from origin for each step.

This will yield 48 points. I'll spare you with a printout of the values :-)

Maidenhead locator 

Some relay lists don't specify correct coordinates, just the coarse Maidenhead approximation. That's good enough for HF propagation models, but not good enough for VHF, UHF, EHF etc propagation models. Nevertheless you'll from time to time be in need to convert those coordinates into each other:

    QPointF lonLat;
    bool ok = fromMaidenhead("JO40JG", lonLat);
    qDebug("ok %d", ok);
    if (ok)
            qDebug("  %f,%f", lonLat.x(), lonLat.y());

    QString s = toMaidenhead(lonLat);
    qDebug("locator: %s\n", qPrintable(s));

This first converts "JO40JG" into 8.791667,50.270833, and than converts this back into a Maidenhead locator.

Map display 

Sooner or later you'll want to show propagation data into some map form. Splat! had only very limited support for this, by reading a data file with coordinates and text, and rendering that into the map.

Slippy map 

So I wrote a simple map widget that can display pre-rendered map data. This are tiles of 256 times 256 dots with map data. Google maps uses them, Openstreetmap uses them, and a of other sources as well. They come usually with zoom-levels between 1 and 18. In this last zoom level you can already see individual houses.

I wrote a simple widget, SlippyMapWidget.h / SlippyMapWidget.cpp, that can efficiently downloads and displays those tiles.

Of course you can interactive zoom, pan, and resize this widget.

Then I wrote SlippyMapParts, that draw stuff over the slippy map. I wrote some cross-hair, zoom-buttons, location display, and great-circle-path display draw parts. Here's an example:

Map sources 

In case you don't like the default OpenStreetMap tiles, I wrote also a map chooser, which can make the slippy map look like this:

It knows about various OpenStreetMap maps like Mapnik, Tiles-at-Home, Cloudmade, ÖPNV, Reit- und Wanderkarte, but it can also display map or satellite imagery from Google. As this is open source, you can add any map source by yourself if you know the URLs for the tiles.

Terrain data 

Without terrain data, you could only calculate free space loss between a transmitter and receiver. Only with terrain data you would know if there is a line-of-sight, if the Fresnel zone is free, if there diffraction or reflection can happen. And a bunch more.

So we need terrain data.

Splat! uses SRTM v3 data. This has the benefit that the data is free for all. However, it comes in an unwieldy format: text files ending with .hgt. If you need to get the elevation of arbitrary points inside such a file, you'd need to open the file, and parse it completely until you find that point. This is what GDAL does, and what makes GDAL access to .hgt files dead slow.

Splat! converts this file into it's own binary file. Basically you get a file with a big array of with 1200 times 1200 short integers. Splat! reads the whole file into memory. There it stays and occupies RAM. It would have been better to use mmap().

I chosed SRTMv4.1 file. This data is free for noncommercial uses, which is true for my usage. It also comes in a nice GeoTIFF format. GeoTIFF can store any rasterized spatial data under the sun, so it's both very versatile ... and complex. GDAL (with the help of LibTIFF can open it out-of-the box. Unfortunately, this combination is still slow when you want to determine the heights of arbitrary points.

However, I detected an interesting property of those SRMVv4.1 GeoTIFF-files: despite the complex GeoTIFF formats, the actual elevation data lies in a consecutive area inside the .TIF file. So I wrote SRTMv4.h / SRTMv4.cpp. It's main job is to parse the TIFF and GeoTIFF headers inside the file, to make sure that it has the right size, correct number of bands, and that the data really is organized contiguously. Once that is determined, accessing the data is very simply:

    // somewhere in the initialization phase:
    map = f.map(dataofs, length * width * sizeof(quint16));

    // and later you can retrieve
    qint16 SRTMv4::get(quint32 x, quint32 y)
    {
            quint32 pos = y*width + x;
            pos *= sizeof(quint16);
            return * ((qint16*) (map+pos));
    }

I omitted all the error checking code in this example, but you get the picture. mmap() uses the file system cache from Linux, so it's usually much faster than reading all the data into memory. First, because it only reads those pages in that are needed. And second, because memory pressure is now reduced.

For you as a user, accessing elevation data is even simpler:

    qDebug("Woellstadt (DH3HS): %d m", srtm.get(DH3HS));
    qDebug("Feldberg (DB0FT):   %d m", srtm.get(DB0FT));

and, provided that you downloaded the SRTM file and put it into the data directory, you'll get back 127 and 866 meters.

Elevation profile 

Now I have together all the basics to make elevation profiles. For example, just 8 km from my QTH is a 2m relay DB0ZAV, which I cannot work. Why is that?

I wrote a simple widget, ProfileWidget.h / ProfileWidget.cpp, which can do the work for me. Here is the code to use it:

GeoData data;
    MakePath(DH3HS, DB0ZAV, data, 100);

    for (int i=0; i<data.posCount(); i++)
            data.e.append(srtm.get( QPointF(data.lon.at(i), data.lat.at(i)) ));

    ProfileWidget *pl = new ProfileWidget(&data);
    pl->show();

And this will yield something like:

On the left is my QTH, on the right is DB0ZAV. And look at that hill. No question that my little CT-270H doesn't come thought.

I must however admit that this picture isn't that informational, but please be patient, it gets better :-)

Terrain classification 

There are several ways to classify terrain. While I still want to implement Longley-Rice in an efficient, I allowed me a sidestep. From some Chinese website you can download almost all ITU recommendations. The P.xxxx recommendations are about radio propagation. An interesting material to read, much nicer than the typewriter-written ITM article I mentioned above.

From those recommendation, ITU-R P.1812 caught my eye. It contains terrain classification material which is from ITU-R P.452. And for this one, there's actually an .XLS file that you can download and which OpenOffice happily digests.

Calculation 

I just took the terrain classification from this example and used it as a basis for writing ITUTerrain.h / ITUTerrain.cpp.

This ITU Algorithm takes the real radius, and multiplies it with a factor k to yield a_e:

    a_e:     9617.8 km    Medium effective earth radius value

Then it adds the antenna heights to the first last elevation data, giving h_ts and h_rs:

    h_ts:      50.0 m     Antenna height AMSL at Tx
    h_rs:     193.0 m     Antenna height AMSL at Rx

AMSL stands for "above mean see level".

Now it calculates some elevation angles theta_d and theta_max. The first is the angle from transmitter to receiver. The second is the maximum angle from the transmitter to every terrain point.

    theta_d: -4.4 mrad
    theta_max: -0.6 mrad

If theta_max is higher than theta_d, then some terrain obstructs the direct view from transmitter to receiver, and we have a trans-horizon condition. theta_max is then also automatically theta_t, the angle from the transmitter to it's horizon:

    Trans horizon? 1
    theta_t:   -0.6 mrad  Horizontal elevation angle (above local horizon) at Tx

If we have a theta_t, we also have a theta_r:

    theta_r:   -1.4 mrad  Horizontal elevation angle (above local horizon) at Rx

And from knowing what obstructs out view, we can calculate d_lt and d_lr, the great-circle distances to the obstructing terrain from Transmitter or Receiver:

    d_lt:      28.0 km    Great circle distance from Tx to it's horizon point
    d_lr:      11.0 km    Great circle distance from Rx to it's horizon point

For some calculations, it's necessary to reduce the real, complex terrain into a straight line approximation. This is called "Smooth earth surface".

We start by calculate h_a, a geometric mean of the real path heights:

    h_a:       38.3 m     Mean of the real path heights

Now we determine the slope m of the line approximation:

    m:          0.6 m/km  slope of least-squares surface relative to sea level

and can the determine the start h_stand end h_sr elevation of this line:

    h_st:       5.4 m     Height of smooth-earth AMSL at Tx
    h_sr:      71.1 m     Height of smooth-earth AMSL at Rx

Now, for every elevation data point, we look if the real elevation is higher than the elevation of this smooth earth approximation. The highest value of this gives us h_m, a terrain roughness value:

    h_m:      119.5 m     Terrain roughness

Visualization 

Here is the terrain from the ITU sample spreadsheet:

You see blue: the two antennas. In black you see the straight line from the transmitter antenna to the receiver. On this line written is the great-circle distance between those locations. You see that near the receiver a hill blocks this line-of-sight.

In brown you see the terrain, and in dark brown the area that is beyond the radio horizon of both the transmitter and sender. That's actually the area between d_lt and d_lr.

In light gray you see the smooth earth approximation line. And to the right you see the terrain roughness value h_m, which equals 120 m.

This is the terrain between my home QTH and the "Wasserkuppe" Relais, DB0WAS.

And here the terrain between the relays on the "Grosser Feldberg", DB0FT, and my home QTH.

And, for reverence, again the more pretty terrain profile between my home QTH and DB0ZAV:

To be continued 

I'm not at the end yet, more is to come. :-)

 

Made on 2010-12-14 with Webber.