[Contents]
Copyright © 2011 jsd

Getting Started with Geographic Information Systems

1 Introduction

There are a number of open-source Geospatial Information System packages available. In recent years they have become powerful enough and reliable enough to do some interesting things.

One such package is GRASS, which stands for "Geographic Resources Analysis Support System”. Another is QGIS, which stands for Quantum GIS. Yet another is GDAL, which stands for Geospatial Data Abstraction Library. Here is a brief summary of the good news and bad news about these packages.

++ At long last, GRASS has become reliable enough and efficient enough to be useful. You can create some useful and beautiful maps with it.
++ GRASS is rather cleanly divided into a “front end” that provides a graphical user interface (GUI), and a “back end” that computes the images. If you don’t want to use the front end, the back end can be controlled using a command-line interface (CLI).
–– Some features you might like GRASS to provide have never been implemented, and are not likely to be implemented anytime soon.
–– Some GRASS features that have allegedly been implemented are quite buggy. Some features don’t do anything at all. Some features run wildly amok. See section ‍7. The front end is conspicuously buggier than the back end.
–– Although some GRASS features are quite efficient, others are still suprisingly inefficient. For example, panning or scrolling the image even 1mm causes GRASS to recompute and repaint the entire image, which can take several seconds. This very inferior in comparison to the state of the art, for instance computer games and flight simulators which can smoothly pan and scroll through very complex scenes.
++ QGIS is basically just a frontend with no backend of its own. QGIS is much more reliable than the GRASS frontend. QGIS graphics updates are orders of magnitude faster. QGIS can reproject vector data on-the-fly.
++ QGIS knows how to talk to the GRASS backend. So the recommended practice in most cases is to install GRASS, install QGIS with the QGIS-GRASS plugin, and then use QGIS (ignoring the GRASS frontend).
++ QGIS and GRASS are licensed under the GPL. They are available for almost any platform (Linux, Mac, Windows).
–– The GRASS documentation is what I call wizard-friendly. Alas, it is quite unfriendly to non-experts.

In particular, if you happen to be looking over the shoulder of an expert, you can see what commands he is using, and then look into the reference manual to find out, roughly, what those commands are supposed to do. The manual is a mapping from commands to ideas.

However, in the real world, most people cannot get started by looking over the shoulder of an expert. They start with an idea, and then they want to find out which commands, if any, allow them to carry out the idea. This requires searching through the entire manual, which is very laborious and inefficient. It is something of a chicken-and-egg problem: if you knew which command to use you could look it up, but since you don’t, you can’t.

What we need is a good tutorial. I benefitted from the data conversion and import notes at reference ‍1. There is also the famous newsletter article reference ‍2. Some documentation on the Lambert Conic Conformal Projection can be found in reference ‍3. There is something of a tutorial at reference ‍4 but it didn’t do me much good.

++ Important: You can find out hints about what commands are relevant by using the gis.m GUI and then reading the "GIS.m" output window to see what it thinks it is doing. Some corrections may be needed, but at least you have a hint as to where to start. This partially solves the chicken-and-egg problem.
++ The GDAL package is mostly just a backend. That is, it contains programs (and library functions) to do GIS computations, but doesn’t do much to display the results. There are some application programs (such as gdalwarp) that can be controlled via a command-line interface.

*   Contents

2 Some Examples

Figure ‍1 shows a shaded relief map. Esthetically speaking, I find this rather attractive. Practically speaking, it is good for portraying slope, but not particularly good at portraying height. That is, it shows you what’s steep and what’s not, but it doesn’t directly show you what’s high and what’s not.

shaded
Figure ‍1: Shaded Relief Map

Figure ‍2 shows a color-coded elevation map. This is a relatively good way of portraying height.

colored
Figure ‍2: Color-Coded Elevation Map

We can combine these two ideas as shown in figure ‍3.

draped
Figure ‍3: Shaded and Color-Coded Elevation Map

If we are particularly interested in certain elevations, we can add contour lines, as shown in figure ‍4. The 3000-foot, 6000-foot, and 9000-foot levels are displayed.

draped-con
Figure ‍4: Shaded and Color-Coded Elevation, with Contour Lines

We can indicate the names and locations of some of the mountain peaks, as shown in figure ‍5. In this diagram, the shaded relief is turned off, so we just see the color coded elevations, with no direct indication of slope. Turning off the shaded relief makes the labels easier to read.

colored-con-peaks
Figure ‍5: Color-Coded Elevation, with Contour Lines and Peak Labels

We can plot the identifiers and locations of the local airports, along with the direct route of flight from KAVQ (Marana) to E60 (Eloy) as shown in figure ‍6. Also indicated by the light magenta shading is the route corridor, i.e. everything within 4 nautical miles of the route of flight. Also, the terrain is not plotted if it is below 3000 feet.

Important point: It is not necessary to show everything on a single map. One of the nice things about a GIS system is that you can conveniently turn off one layer and turn on another. This keeps the map from becoming too cluttered.

This stands in contrast to paper maps, where people often try to maximize the amount of information while minimizing the number of maps that they need to carry ... resulting in maps that are hard to interpret because they are too cluttered.

colored-apt-direct
Figure ‍6: Color-Coded, with Airports and Direct Route

In case you were wondering, the point of the exercise is that we can plan a route that goes between the peaks rather than over them, as shown in figure ‍7. Also in this map, shown in blue, are the named navigation fixes known to the FAA ... plus a fix that I defined myself, namely the “pass” fix, which I use to specify the route of flight.

colored-fix-pass
Figure ‍7: Color-Coded, with Fixes and Route Through Pass

3 Fundamental Concepts

We begin by reviewing some basic concepts, and then gradually develop some more advanced concepts.

Note the following contrast:

When dealing with vector data, you don’t need a full-blown SRS. In particular you don’t need to specify any definite cartographic projection. If you know that a certain point is at such-and-such latitude, longitude, and elevation (relative to some specified datum) that’s all you need to know. That point is where it is.   When dealing with raster data, the data is essentially a flat, rectangular array of pixels. In order to know what those pixels mean, you need a way of projecting them onto the surface of the globe. Therefore, to interpret raster data you do need a full-blown SRS, including a projection. That is, you need to project the pixels onto the surface of the earth.

Meanwhile, there is another very different, very important contrast to be made: A GIS package has both a back-end (for storing data and doing computations) and a front-end (for making maps and displaying them to humans).

The projection used by the front-end
does not need to be the same as
the projection (if any) used by the back-end.
 ‍ ‍ ‍ ‍ ‍

Various Universal Transverse Mercator (UTM) projections are widely used for representing raster data in memory and in files. This is perfectly reasonable, except at high latitudes. There are standards for this, which allow for compatibility when storing and communicating data.   There is no excuse for using any Mercator-like projection for making actual maps that will be seen by human eyes (except maybe for small areas near the equator). Just because your GIS system makes this easy to do doesn’t mean it is an acceptable practice.

Let’s be clear: For storage and computation involving raster data, UTM is often OK.   It is also perfectly OK to do storage and computation using non-Mercator projections. For example, very nice Sectional Aeronautical Charts are available on line, as raster images using Lambert Conformal Conic projections. See reference ‍5.

For storage and computation involving vector data, you don’t need any projection at all.   For making maps that will be seen by humans, Mercator projections are almost never OK. You need to find a way to reproject the data into some better-behaved projection.

You don’t want to run around reprojecting raster data back and forth unnecessarily, because it is computationally expensive and smudges the data a little bit.

4 Reprojection

GRASS has commands to reproject raster data. That means converting data from one projection to another. An example is converting from a UTM projection to a Lambert Conformal Conic projection. You can also reproject data using the gdalwarp program provided by the GDAL package.

Using these commands is a royal pain. However, taking shortcuts is not recommended, because the tricks that are easy aren’t accurate, and the tricks that are accurate aren’t particularly easier than using the GRASS reprojection commands.

One drawback of gdalwarp is that the rasters it produces measure positions in meters relative to some weird datum, rather than using degrees of latitude and longitude. This is inconvenient when trying to query the map. Also: I don’t see why this is necessary. It should be possible to find your way around a Lambert projection using lat/lon.

This means that if you are using any reasonable projection, i.e. not Mercator, and you want to read in a kml file or anything else that uses lat/lon, GRASS insists that you switch to a UTM "location", read in the file, switch back to your original "location", and reproject the data.

The script krunch does this. It’s a mess, but since the script takes care of it, I can live with it.

QGIS knows how to reproject vector data on-the-fly, so you don’t need to mess with GRASS for this. (As discussed in section ‍3, vector data doesn’t really have an associated projection anyway, so technically QGIS is just projecting (not reprojecting) the data.) To make this work, you need to define a custom Coordinate Reference System for your project (menu -> Edit -> Custom CRS...). Give it a name and a definition, then push the save button (second-to-last button on the right).

Here is how to define the projection used for the examples in section ‍2 :

+proj=lcc +lat_1=31 +lat_2=33 +lat_0=0 +lon_0=-111
 ‍ ‍ ‍ ‍ ‍ +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

Then go to (menu -> File -> Project Properties –> Coordinate Reference System). Check the box that says ‘Enable on-the-fly CRS transformations.’ Then go to the list-box and select the desired transformation.

There are also some options under (menu -> Edit -> Options) concerning prompting for a CRS when needed, when reading raster data. Some raster data (e.g. GTiff format) is internally tagged with metadata describing the projection it is using, but plenty of other data is not.

Also: You don’t need to fuss with reprojections when using google-earth as your frontend. It knows how to read vector data from .kml files and project it according to whatever projection it is using at the moment. The free versions of google-earth won’t read in raster data at all, so you don’t need to worry about reprojecting that.

5 Other GIS Systems

Some things that are hard to do in GRASS are easier to do in google-earth. For example, google-earth does a lot of nice things with 3D views. The data in GRASS is of course 3D, but the viewer views it from directly above. I don’t know of any plans to support viewing at an angle.

Conversely, some things that are hard to do in google-earth are relatively easy in GRASS. Example: computing contour lines.

There are commercial versions of google-earth (starting at several hundred dollars per license) that are more capable than the freeware version. For example, you could compute a contour line in GRASS, export it, and then import it into the commercial version of google-earth.

ArcGIS is a long-established widely-used highly-capable commercial GIS package.

Mapnik is a free toolkit for developing mapping applications. See reference ‍6.

QGIS is free and open-source. It has a good front-end for displaying things. It has a plugin for communicating with GRASS, so you can use the GRASS back-end for computing things.

6 How To Create Some Maps

6.1 Preliminary Setup

The GRASS commands to create the maps shown in section ‍2 are so complicated that I wrote a perl script to do the work. It is shown in krunch.

  1. Before you can even get to the grass-shell the first time, you need to create a working directory ("grassdata") because grass will ask you for it.
  2. GRASS uses the word "location" with a peculiar idiomatic meaning. What GRASS calls a "location" is more properly called a subdirectory or mapset or some such. In addition to containing zero or more maps, a GRASS "location" is associated with a specific Spatial Reference System (SRS) aka Coordinate Reference System (CRS). The concept of SRS is discussed in section ‍3. Do not confuse a GRASS "location" with a GRASS "region" (which is something else entirely). Also not confuse a GRASS "location" with an ordinary location i.e. an ordinary position in space.

    When starting GRASS, before you can do anything else, GRASS requires you to choose some GRASS "location". If there are no pre-existing locations, you are required to create one, which implies choosing a SRS (aka CRS) for the new "location". (This requirement is partially but not entirely logical, for reasons discussed in section ‍3.) You are forced to choose something, but the choice is not super-important, because you can change it later, using the "mapset ..." command.

    Suggestion: If this is a first-time startup in an empty directory, push the EPSG button and use code 4326. That’s the EPSG identifier for UTM (Universal Transverse Mercator). That CRS is used by Shuttle Radar Terrain Mapping (SRTM) data files and lots of other things. This is a reasonable choice for data storage, communication, and computation ... but not for map-making, as discussed in section ‍3.

  3. GRASS has a notion of "region". This is basically a clipping region. You need to set the region correctly, or all sorts of things will fail, and it won’t be obvious what the problem is. Some (but not all) commands act only on data within the region. Hint: to find out about the region: g.region -p

    This is one of the reasons why I knew I had to write a script. The script systematically calculates the region and tells GRASS about it. Keep in mind that the GRASS “region” is a very different concept from the GRASS “location”.

  4. Once GRASS is well and truly started, there is GUI and also a shell command prompt, in the window that you used to start GRASS. It has pushed a subshell. The script krunch is normally run in this GRASS shell. (It can be run elsewhere if you set the PATH and LD_LIBRARY_PATH properly.)

    There are a number of commands that will fail in the GRASS shell, including most of the commands that start with “d.” such as d.font and d.rast. However, see also next item.

  5. The GUI can provide a window called “Output - GIS.m” although this does not initially appear. You can cause it to appear by pushing certain buttons on the GUI. Perhaps the easiest way is to go to the main “GIS Manager” GUI window and push Raster –> Map Calculator –> Run.

    This “Output - GIS.m” window is not just for output; it has a little input window near the bottom where you can type commands, rather like in the aformentioned GRASS shell, except that the display is defined in this environment, so some of the “d.” commands will work, such as d.info. Most of them still don’t work, as discussed in section ‍7, but at least some of them do.

    For example, d.mon -L in the GRASS shell window reports that all possible monitors are “not running”. However, the same command in the “Output GIS.m” window reports them all as running.

  6. You need some data.

    In the examples in this document, I started with elevation data from SRTM i.e. the Shuttle Radar Topography Mission, reference ‍7. The script knows how to download this data in some cases.

    I also used data for airports, navaids, and navigation fixes from the X-plane databases, reference ‍8. You need to download this by hand, if you don’t already have it lying around.

6.2 GRASS Commands Used by the Script

Here is a concise listing of all the commands used by the krunch script to generate the map layers used in the examples shown in section ‍1.

Script started on Thu 22 Sep 2011 08:13:58 PM MST
+ ./krunch -f
: +++; touch N32W111.force
: +++; mosaic='mosaic' tiles='N32W111 N32W112 N31W111 N31W112' make -f tiles.make
wget -c http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/N32W111.hgt.zip
--2011-09-22 20:13:58--  http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/N32W111.hgt.zip
Resolving dds.cr.usgs.gov... 152.61.128.95
Connecting to dds.cr.usgs.gov|152.61.128.95|:80... connected.
HTTP request sent, awaiting response... 200 OK

    The file is already fully retrieved; nothing to do.

rm -f N32W111.force
: +++; r.in.gdal input=mosaic.tif output=mosaic  --q
: +++; g.region rast=mosaic
: +++; g.proj --q -c georef=mosaic.tif
: +++; r.shaded.relief --q map=mosaic zmult=10
: +++; r.contour --q input=mosaic output=mosaic_con3 levels=914.4
: +++; r.contour --q input=mosaic output=mosaic_con6 levels=1828.8
: +++; r.contour --q input=mosaic output=mosaic_con9 levels=2743.2
: +++; v.out.kml --q -l mosaic_con3 size=2 color=255:255:0
/usr/bin/ogr2ogr --formats
WARNING: Default driver / database set to:
         driver: dbf
         database: $GISDBASE/$LOCATION_NAME/$MAPSET/dbf/
WARNING: Vector map <mosaic_con3> is 3D. Use format specific layer creation
         options (parameter 'lco') to export in 3D rather than 2D (default)
 -f KML mosaic_con3.kml.raw mosaic_con3.sh
: +++; g.mapset mapset=PERMANENT location=UTM_4326_SRTM
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r
/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history
: +++; v.in.ogr --q -z dsn=some-apt.kml out=kapts
WARNING: Width for column Name set to 255 (was not specified by OGR), some
         strings may be truncated!
WARNING: Width for column Description set to 255 (was not specified by
         OGR), some strings may be truncated!
WARNING: Width for column Name set to 255 (was not specified by OGR), some
         strings may be truncated!
WARNING: Width for column Description set to 255 (was not specified by
         OGR), some strings may be truncated!
WARNING: Width for column Name set to 255 (was not specified by OGR), some
         strings may be truncated!
WARNING: Width for column Description set to 255 (was not specified by
         OGR), some strings may be truncated!
WARNING: Cleaning polygons, result is not guaranteed!
: +++; g.mapset mapset=PERMANENT location=z111
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r /home/jsd/grassdata/z111/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/z111/PERMANENT/.bash_history
: +++; v.proj --q location=UTM_4326_SRTM input=kapts
: +++; g.mapset mapset=PERMANENT location=UTM_4326_SRTM
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r
/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history
: +++; v.in.ogr --q -z dsn=some-track.kml out=ktrack
WARNING: Width for column Name set to 255 (was not specified by OGR), some
         strings may be truncated!
WARNING: Width for column Description set to 255 (was not specified by
         OGR), some strings may be truncated!
: +++; g.mapset mapset=PERMANENT location=z111
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r /home/jsd/grassdata/z111/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/z111/PERMANENT/.bash_history
: +++; v.proj --q location=UTM_4326_SRTM input=ktrack
: +++; g.mapset mapset=PERMANENT location=UTM_4326_SRTM
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r
/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history
: +++; v.in.ogr --q -z dsn=some-corridor.kml out=corridor
WARNING: Width for column Name set to 255 (was not specified by OGR), some
         strings may be truncated!
WARNING: Width for column Description set to 255 (was not specified by
         OGR), some strings may be truncated!
WARNING: Cleaning polygons, result is not guaranteed!
: +++; g.mapset mapset=PERMANENT location=z111
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r /home/jsd/grassdata/z111/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/z111/PERMANENT/.bash_history
: +++; v.proj --q location=UTM_4326_SRTM input=corridor
: +++; g.mapset mapset=PERMANENT location=UTM_4326_SRTM
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r
/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/UTM_4326_SRTM/PERMANENT/.bash_history
: +++; v.in.ogr --q -z dsn=other-corridor.kml out=othercorridor
WARNING: Width for column Name set to 255 (was not specified by OGR), some
         strings may be truncated!
WARNING: Width for column Description set to 255 (was not specified by
         OGR), some strings may be truncated!
WARNING: Cleaning polygons, result is not guaranteed!
WARNING: 119 areas represent more (overlapping) features, because polygons
         overlap in input layer(s). Such areas are linked to more than 1
         row in attribute table. The number of features for those areas is
         stored as category in layer 2
: +++; g.mapset mapset=PERMANENT location=z111
Erasing monitors...
Cleaning up temporary files...
WARNING: Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r /home/jsd/grassdata/z111/PERMANENT/.bash_history;
HISTFILE=/home/jsd/grassdata/z111/PERMANENT/.bash_history
: +++; v.proj --q location=UTM_4326_SRTM input=othercorridor
Reprojecting primitives: 
Script done on Thu 22 Sep 2011 08:14:15 PM MST

krunch.log  ‍  ‍ (DL)

6.3 A Script to Produce Maps Using GRASS

#! /usr/bin/perl -w

use strict;
use Symbol;
use List::Util qw[min max];
use POSIX;      ## for floor and ceil
use File::Basename;

# To do:  FIXME:
# r.fillnulls
# Note: color=rainbow is (for now) favored.
# May need to remove gdalwarp output file;  see manpage.

########################################
# some constants:
my $nm = 1853.248;          ## nautical mile (not nanometer)
my $pi = 4*atan(1);
my $deg = $pi/180.;
my $inch = 0.0254;          ## in meters
my $foot = 12*$inch;
my $Q = '"';                ## double quote

##########
# some global settings and mode flags:

my $goodloc = 'z111';

my $recompute = 0;
my $verbosity = 1;
my $verboflag = '--q';
my $vbf = '-q';           # shorter version of verboflag

# current bounding box, in degrees:
  my $left = 9e99;
  my $right = -9e99;
  my $top = -9e99;
  my $bot = 9e99;
  my $lrmid = 0;

# bounding box, in meters:
##??? [not used]  my ($north, $south, $east, $west);

######################################################################
# functions

sub xystem {
  my ($cmd, $vlevel) = @_;
  if (! defined $vlevel) {
    $vlevel = 1;
  }
  if ($vlevel <= $verbosity) {
    print ": +++; $cmd\n";
  }
  system $cmd;
}

sub g_gisenv_get {
  my ($var) = @_;
  my $cmd = "g.gisenv get=$var";
  my $pipe = Symbol::gensym;
  open($pipe, "-|", $cmd)
    || die "Could not pipe from '$cmd'\n";
  my $rslt = <$pipe>;
  chomp $rslt;
  close $pipe;
  return $rslt;
}

sub map_exists {
  my ($mapname) = @_;
##  my $cmd = "g.mlist -m type=all";
  my $cmd = "g.list -m type=all";
  my $pipe = Symbol::gensym;
  open($pipe, "-|", $cmd)
    || die "Could not pipe from '$cmd'\n";
  liner: while (my $line = <$pipe>){
    chomp $line;
    $line =~ s'@.*'';
    if ($line eq $mapname) {
      return 1;
    }
  }
  return 0;
}

sub read_kml {
  my ($fname, $table) = @_;
  my $utm_location = "UTM_4326_SRTM";
  my $loc = g_gisenv_get('LOCATION_NAME');
  if ($loc ne $utm_location) {
    if (! -d $utm_location){
      print STDERR "Location does not exist: $utm_location\n";
      return;
    }
    xystem "g.mapset mapset=PERMANENT location=$utm_location";
  }
  if (map_exists($table)) {
    xystem "g.remove -f $verboflag type=vector name=$table";
  }
#?? # set the region before doing the import:
#??   xystem "g.region w=$left s=$bot e=$right n=$top";
  xystem "v.in.ogr $verboflag input=$fname output=$table";
  if ($verbosity >= 2) {
    xystem "g.proj -p";              # info about current projection
    xystem "g.region -p";            # basic info about region
                                     # including projection
##    xystem "g.mlist -t -m type=all"; # list available map layers
    xystem "g.list -t -m type=all"; # list available map layers
                # in the format "type/layer@MAPSET"
    xystem "v.info map=$table";      # info about particular map
  }
  print "..............\n";

  if (! -d $goodloc){
    print STDERR "Location does not exist: $goodloc\n";
    return;
  }
  xystem "g.mapset mapset=PERMANENT location=$goodloc";

  if (map_exists($table)) {
    xystem "g.remove -f $verboflag type=vector name=$table";
  }
  xystem "v.proj $verboflag location=UTM_4326_SRTM input=$table"
}

sub mk_contour {
  my ($mosaic, $suffix, $levels) = @_;
  my $table = "${mosaic}_$suffix";
  if (map_exists($table)) {
    xystem "g.remove -f $verboflag type=vector name=$table";
  }
  xystem "r.contour $verboflag input=$mosaic output=$table levels=$levels";
}

# Compute an overall bounding box to cover all tiles:
sub get_box_alltiles {
  my @tiles = @_;

  for my $tile (@tiles) {
    my $cmd = "gdalinfo -nogcp -nomd -noct $tile.hgt";
    my $pipe = Symbol::gensym;
    open($pipe, "-|", $cmd)
      || die "Could not pipe from '$cmd'\n";
    liner: while (my $line = <$pipe>){
      chomp $line;
      $line =~ s'[(),]' 'g;
      my @stuff = split (' ', $line);
      if (! defined $stuff[1]) {
        next liner;
      }
      if ($stuff[0] eq 'Upper') {
        $top = max($top, $stuff[3]);
      }
      if ($stuff[0] eq 'Lower') {
        $bot = min($bot, $stuff[3]);
      }

      if ($stuff[1] eq 'Right') {
        $right = max($right, $stuff[2]);
      }
      if ($stuff[1] eq 'Left') {
        $left = min($left, $stuff[2]);
      }
    }
  }
# round to nearest whole degree:
  $top   = floor($top + .5);
  $bot   = floor($bot + .5);
  $left  = floor($left + .5);
  $right = floor($right + .5);
  $lrmid = ($left + $right) / 2;
  print "alltiles: $top > $bot ... $left < $right : $lrmid\n";
}

main: {
  my @tiles = ();

  for my $arg (@ARGV) {
    if ($arg eq '-v') {
      $verbosity++;
    } elsif ($arg eq '-f') {
      $recompute++;
    } elsif ($arg =~ '^-') {
      die "Unrecognized option '$arg'\n";
    } else {
      push @tiles, $arg;
    }
  }

  if (!@tiles) {
    @tiles = ("N32W111", "N32W112", "N31W111", "N31W112");
  }
  my $alltiles = join(' ', @tiles);
  if ($verbosity > 1) {
    print "Tiles are: $alltiles\n";
  }

  if ($verbosity > 2) {
    $verboflag = '--v';
    $vbf = '';
  }

  if (basename(getcwd) ne 'grassdata') {
    print STDERR "Please cd to 'grassdata'\n";
    exit(1);
  }

  my $loc = g_gisenv_get('LOCATION_NAME');
  if ($loc ne $goodloc) {
    print STDERR "Warning: Location ($loc) should be $goodloc\n";
  }

## Data sources:
## http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/N32W111.hgt.zip
## Bigger, not necessarily better:
## http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N32W112.hgt.zip

my $america = "http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America";
my $eurasia = "http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia";

  my $mosaic = "mosaic";

  if ($recompute) {
    xystem "touch $tiles[0].force";
  }
  xystem "mosaic='$mosaic' tiles='$alltiles' make -f tiles.make";

  {
# read it into grass
    my $newloc = '';
    if (-d "$goodloc/PERMANENT") {
      # location exists; switch to it
      my $loc = g_gisenv_get('LOCATION_NAME');
      if ($loc ne 'z111') {
        xystem "g.mapset mapset=PERMANENT location=$goodloc";
      }
    } else {
      # location doesn't exist; create it from file
      $newloc = "location=$goodloc";
    }
    if ($recompute && map_exists($mosaic)){
      xystem "g.remove -f $verboflag type=raster name=$mosaic";
    }
    my $slurp = "r.in.gdal input=$mosaic.tif output=$mosaic";
    $slurp .= " $newloc $verboflag";
    if ($recompute || $newloc ne '' || !map_exists($mosaic)) {
      xystem "$slurp";
    }

# and set the region accordingly:
# this is important:
    xystem "g.region rast=$mosaic";
# Also set the projection accordingly.
# This happens automatically when the region gets created,
# but not if we read a new file into an old region, so
# we really need to do this explicitly.
# Also note that the -c is important and non-obvious:
    xystem "g.proj $verboflag -c georef=$mosaic.tif";

# Create shaded-relief representation.
# Note that the contour routine reports that
# the initial color table is set to "grey"

    my $sname = "$mosaic.shade";
    my $shadex = map_exists($sname);
    if ($recompute && $shadex) {
      xystem "g.remove -f $verboflag type=raster name=$sname";
      $shadex = 0;
    }
    if (!$shadex) {
      xystem "r.shaded.relief $verboflag map=$mosaic zmult=10";
    }

    # make kml file containing contour lines
    if (0) {
      mk_contour($mosaic, "con3", 3000*$foot);
      mk_contour($mosaic, "con6", 6000*$foot);
      mk_contour($mosaic, "con9", 9000*$foot);
      xystem ("v.out.kml $verboflag -l ${mosaic}_con3 size=2 color=255:255:0");
    }
  }

  {     # create shaded maps
    my $cut = 3000*$foot;
    xystem "r.mapcalc --overwrite $Q${mosaic}_3k "
      . " = if($mosaic >= $cut, $mosaic, null())$Q";
    xystem "r.mapcalc --overwrite $Q${mosaic}_3k.shade "
      . " = if($mosaic >= $cut, ${mosaic}.shade, null())$Q";
  }

  read_kml("some_fix.kml", "some_fix");
  read_kml("some_apt.kml", "some_apts");
  read_kml("some_track.kml", "some_track");
  read_kml("some_corridor.kml", "some_corridor");
  read_kml("direct_corridor.kml", "direct_corridor");
}

krunch  ‍  ‍ (DL)

And here is a makefile that fetches data, builds a mosaic, removes nulls, and reprojects from UTM to Lambert.

# Key variables, defined in environment:
#   $tiles (input), $mosaic (principal output)

# Other useful variables
#   $projection, $verbosity, $CRS (i.e. GRASS "location")

ifndef CRS
  CRS = UTM_4326_SRTM
endif

ifndef tiles
  tiles = N31W111 N31W112 N32W111 N32W112
endif

ifndef projection
  projection = $(shell ./gis-bbox $(tiles))
endif

# for reference:
lcc_projection = +proj=lcc +lat_1=31 +lat_2=33 +lat_0=0 +lon_0=-111   \
     +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

ifndef mosaic
  mosaic = mosaic
endif
m = $(mosaic)

ifndef verbosity
  verbosity = -q
endif

america = http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America
eurasia = http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia

# see if location already exists;
# if not, construct phrase for creating it:
setloc = $(shell if ! test -d $(CRS)/PERMANENT ; then \
         echo location=$(CRS) ; fi)

hgt_files = $(tiles:%=%.hgt)

zip_files = $(hgt_files:%=%.zip)

tif_files = $(tiles:%=%.tif)

force_files = $(tiles:%=%.force)

%.tif : %.hgt
	gdalwarp $(verbosity) -of GTiff -s_srs EPSG:4326 \
           -t_srs EPSG:4326 $< $@
	#? rm $<        # .hgt file no longer needed

%.hgt : %.hgt.zip
	unzip -o $<
	touch -r $@ $@.creation_time
	touch $@        # time of unzipping, not time of creation
	#? rm $<        # .zip file no longer needed

# FIXME: should check whether america or eurasia:
%.hgt.zip : %.force
	wget -c $(america)/$@
	rm -f $<

.SECONDARY : $(hgt_files) $(zip_files) $(force_files)

.PHONY : all

all : $(mosaic).tif

$(mosaic)_utm.tif : $(tif_files)
	gdalwarp $(verbosity) -of GTiff -s_srs EPSG:4326 \
 -t_srs EPSG:4326 $^ $@

$(mosaic)_filled.tif : $(mosaic)_utm.tif
	if test -z "$(newloca)" \
  -a "_$$(g.gisenv get=LOCATION_NAME)" != "_$(CRS)" ; then \
  g.mapset mapset=PERMANENT location=$(CRS) ; fi
	r.in.gdal $(verboflag) input=$< output=$(mosaic)_utm \
  $(setloc) --overwrite
	g.region rast=$(mosaic)_utm
	r.mapcalc "$m_utm = if($m_utm == 0, null(), $m_utm)"
	r.fillnulls input=$(mosaic)_utm output=$(mosaic)_filled --overwrite
	r.out.gdal format=GTiff input=$(mosaic)_filled output=$@

$(mosaic).tif : $(mosaic)_filled.tif
	gdalwarp $(verbosity) -of GTiff -s_srs EPSG:4326 \
  -t_srs '$(projection)' \
  $^ $@

tiles.make  ‍  ‍ (DL)

6.4 A Program to Parse apt.dat and fix.dat Files

// Parse fix.dat and apt.dat files
// select the fixes of interest
// calculate track (leg course and distance ... and corridor)
// and format the results.

// To do:
// *) http://code.google.com/apis/kml/documentation/kml_tut.html
// *) use distance to TRACK rather than distance to target points
// *) clean up box selector
//      -- parse 12e 34w
//      -- swap depending on which is left, which is right
// *) seen_it would be better as a list of Features (not strings).
// *) use actual commas (and quotes) in .csv output
// *) better color map?  custom rules file?
// *) Get rid of restrictors;
//   it's easier and much quicker to write all file-types every time.
// *) 3D corridor, at altitude.

using namespace std;
#include <iostream>
#include <fstream>
#include <GeographicLib/Geodesic.hpp>
#include <cmath>
#include <stdio.h>      /* for popen() */
#include <sstream>
#include <algorithm>    /* for sort */
#include <vector>
#include <set>
#include <map>
#include <iomanip>      /* for setw */
#include <sys/types.h>  /* for stat */
#include <sys/stat.h>
#include <unistd.h>

using namespace GeographicLib;
double a = Constants::WGS84_a<double>();
double f = Constants::WGS84_f<double>();

const Geodesic geod(a, f);
#include "GeoPoint.h"

static const double nm(1853.248);          // nautical mile (not nanometer)
static const double pi(M_PI);
static const double deg(pi/180.);
static const double inch(0.0254);          // measured in SI meters
static const double foot(12.*inch);
static const double bogus(-9e99);

#define X(aaa) aaa,
enum catter {
#include "catter.h"
};
#undef X

map<catter,string> catter_name;

#define X(aaa) catter_name[aaa] = #aaa;
void init_catter_names() {
#include "catter.h"
}
#undef X

class Runway{
public:
  GeoPoint end[2];
  double width;         // in meters
  int sfc;
};

class Feature {
public:
  string ident, fullname;
  GeoPoint loc;
  double dist;
  double fwd;           // in radians
  int has_tower;
  vector<Runway> rwy;
  catter cat;
// Note that category code "feature" indicates
// a generic feature with no definite category.

// constructors:
  Feature(const string _ident, const double _lat, const double _lon)
  : ident(_ident), fullname(""), loc(_lat, _lon, -9e99),
    has_tower(0), cat(feature)
  {}
  Feature()
  : ident("x-x"), fullname(""), loc(0, 0, -9399),
    has_tower(0), cat(phantom)
  {}
};


typedef vector<Feature> VP;

// A Bucket is list of Features ... all with the same ident!
// This is necessary, because fix names are not unique.
// Here are the multiplicity numbers for some examples:
//      9 DELTA
//      8 TANGO
//      8 BRAVO
//      7 TOMAS
//      7 FRANK
//      7 CORAL
//      7 ALPHA
//      6 STONE
//      6 ROCKY
//      6 ROBIN
//      6 GOLFO
//      6 CRANE
//      6 CANDY
//      6 APPLE
//
// VOR names are similarly non-unique.
//      5 SIE
//      4 TRN
//      4 TAN
//      4 SAV
//      4 SAM
//      4 MGA
//      4 CMA
//      4 BNA
//      4 BAL
//
// Whenever a non-unique name is mentioned, we go to a lot
// of trouble to find the one that is nearest the previous fix.
class Bucket : public VP {
public:
  Feature first_for_debugging_only() const {
    Feature rslt;
    if (size() == 1) {
      return front();
    }
    if (!size()) {
      cerr << "? How did we get a bucket of zero size in the dbase?" << endl;
      rslt.cat = phantom;
      return rslt;
    }
///    if (size() > 1) {...}
    { // must be ambiguous
      if (0) cerr << "? Identifier '" << front().ident << "' is "
           << size() << "-fold ambiguous."
           << endl;
      rslt.cat = ambiguous;
      rslt.has_tower = size();
      return rslt;
    }
  }

  Feature nearest_to(const Feature& vicinity) const {
    Feature rslt;
    if (size() == 1) {
      return front();
    }
    if (size() == 0) {
      cerr << "? How did we get a bucket of zero size in the dbase?" << endl;
      rslt.cat = phantom;
      return rslt;
    }
    if (vicinity.cat == phantom){   // ambiguous and not fixable
// don't print error message; let caller do it:
      if (0) cerr << "? Identifier '" << front().ident << "' is "
           << size() << "-fold ambiguous."
           << endl;
      rslt.cat = ambiguous;
      rslt.has_tower = size();
      return rslt;
    }
    double dist_min(9e99), fwd, rev;
    for (VP::const_iterator ptr = begin(); ptr != end(); ptr++) {
      double dist = vicinity.loc.geoDistTo(ptr->loc, fwd, rev);
      if (dist < dist_min) {
        dist_min = dist;
        rslt = *ptr;
        if (0) cerr << "Hit at " << dist / nm
           << "  vlat : " << vicinity.loc.latRad/deg
           << "  vlon : " << vicinity.loc.lonRad/deg
           << endl;
      }
    }
    return rslt;
  }
};

//////////////////////////////////////////////////////////////////////
// various useful functions

// remove \n or \r\n from the end of a string:
string chomp(const string sss){
  string rslt(sss);
  int where;
  where = rslt.length();
  if (where > 0) {
    where--;
    if (rslt[where] == '\n') rslt = rslt.substr(0, where);
  }
  where = rslt.length();
  if (where > 0) {
    where--;
    if (rslt[where] == '\r') rslt = rslt.substr(0, where);
  }
  return rslt;
}

string ltrim(const string sss){
  string rslt(sss);
  size_t where = rslt.find_first_not_of(" \t");
  if (where == string::npos) return rslt;
  return rslt.substr(where);
}

// Compare two featuress according to their "dist"
// field, which presumably indicates how close they
// are to some reference point.
class distcomp{
  int origin;
public:
  bool operator()(const Feature& a, const Feature& b) const{
    return a.dist < b.dist;
  }
  distcomp(int _origin)
  : origin(_origin)
  {
    // planning for the future
  }
};

string toupper(const string sss){
  string rslt(sss);
  for (string::iterator ptr=rslt.begin(); ptr != rslt.end(); ptr++){
    *ptr = toupper(*ptr);
  }
  return rslt;
}

class CaseFold{
public:
  int operator() (const string& a, const string& b){
    return toupper(a) < toupper(b);
  }
};

FILE* open_fix_file(const string _fn, const int oknone=0){
  string fn(_fn);
  string tested("");
  struct stat statbuf;

// This stat business is important, because it is hard
// to catch errors in the zcat command in the pipe.
// Also this makes it easy to add an implicit ".gz" if required.
  int rslt = stat(fn.c_str(), &statbuf);
  tested = fn;
  if (rslt < 0) {
    if (fn.find(".gz")+3 == fn.length()) {
      // already an explicit ".gz"
    } else {
      fn += ".gz";
      rslt = stat(fn.c_str(), &statbuf);
      tested += " " + fn;
    }
  }

  if (rslt < 0) {
    if (oknone) return 0;
    cerr << "Could not stat file(s) '" << fn << "'\n"
        << endl;
    exit(1);
  }

  FILE* pipe(0);
  if (fn.find(".gz")+3 == fn.length()) {
    pipe = popen(("zcat " + fn).c_str(), "r");
  } else {
    pipe = fopen((fn).c_str(), "r");
  }
  if (!pipe) {
    cerr << "Could not open input file '" << fn << "'\n";
    exit(1);
  }
  return pipe;
}

// Fix file version 601 is the same as version 600
// with the _elevation_ inserted as the third field
//
void scan_fixes(const string fn, const int oknone,
        vector<Feature>& foo, const catter setcat = fix){

  FILE* pipe = open_fix_file(fn, oknone);
  if (!pipe) return;

  size_t size(100);
  char* buf = (char*)malloc(size);
  if (getline(&buf, &size, pipe) < 0) {
    cerr << "Oops.\n";
    exit(1);            // can happen on file-not-found
  }
  string sbuff(chomp(buf));

  if (sbuff != "I") {
    cerr << "Unexpected format in fix file '" << fn << "'" << endl;
    exit(1);
  }

// read second line:
  if (getline(&buf, &size, pipe) < 0) {
    cerr << "Can't read second line.\n";
    exit(1);
  }

  int fileversion(0);
  sscanf(buf, "%d", &fileversion);
  if (fileversion != 600 && fileversion != 601) {
    cerr << "Unexpected version: " << fileversion
         << " in fix file '" << fn << "'"
         << endl;
    cerr << "Line was: " << chomp(buf) << endl;
    exit(1);
  }

// skip rest of header,
// i.e. scan for empty line at end of header:
  for (;;) {
    if (getline(&buf, &size, pipe) < 0) {
      break;
    }
    if (chomp(buf).length() == 0) break;
  }

  Feature ppp;
  for (;;){
    if (getline(&buf, &size, pipe) < 0) {
      break;
    }
    //fprintf (stdout, "%s...\n", buf);
    stringstream parse(buf);
    string tmp;

    parse >> tmp;
    tmp = ltrim(tmp);
// Can happen at end of file:
    if (tmp.length() == 0) continue;
    if (tmp[0] == '#') continue;
    ppp.loc.latRad = atof(tmp.c_str()) * deg;
    parse >> tmp;
    ppp.loc.lonRad = atof(tmp.c_str()) * deg;
    if (fileversion == 601) {
      parse >> ppp.loc.elev;
      ppp.loc.elev *= foot;     // convert from feet to meters
    }
    parse >> ppp.ident;
    ppp.cat = setcat;
    foo.push_back(ppp);
  }
  fclose(pipe);
}

// Declare this locally, because
// it adds things in Mercator space,
// which is not safe unless things are very nearby.
// It's OK for present purposes, though.
GeoPoint& operator+=(GeoPoint& a, const GeoPoint& b){
  a.latRad += b.latRad;
  a.lonRad += b.lonRad;
  return a;
}

GeoPoint& operator/=(GeoPoint& a, const double denom){
  a.latRad /= denom;
  a.lonRad /= denom;
  return a;
}

void scan_apts(const string fn, const int oknone,
        vector<Feature>& foo){

  FILE* pipe = open_fix_file(fn, oknone);
  if (!pipe) return;

  size_t size(100);
  char* buf = (char*)malloc(size);
  getline(&buf, &size, pipe);
  string sbuff(chomp(buf));
  if (sbuff != "I") {
    cerr << "Unexpected format in airport file '" << fn << "'" << endl;
    exit(1);
  }

// read second line:
  if (getline(&buf, &size, pipe) < 0) {
    cerr << "Can't read second line.\n";
    exit(1);
  }

  int fileversion(0);
  sscanf(buf, "%d", &fileversion);
  if (fileversion != 810 && fileversion != 850) {
    cerr << "Unexpected version: " << fileversion
         << " in apt file '" << fn << "'"
         << endl;
    cerr << "Line was: " << chomp(buf) << endl;
    exit(1);
  }

// skip rest of header,
// i.e. scan for empty line at end of header:
  for (;;) {
    if (getline(&buf, &size, pipe) < 0) {
      break;
    }
    if (chomp(buf).length() == 0) break;
  }

// for 810
  double lat_deg, lon_deg, elev_foot;
  string fullname, rwylabel;
  double heading, rwylength;
  double displaced, stopway;
  int bldgs;
  int lighting_code;

  int shoulderType;
  double smoothness;
  int centerLights, edgeLights, distanceSigns;
// 850, each end of the runway:
  double wid_ft, overrun;
  int markingType, approachLights, tdzLights, REIL;

  int phase(0);         // 0 means pre-record, looking for record
                        // 1 means we have started but not finished
                        //  processing a record.
  GeoPoint rwyAccum;
  int nRwyEnds;
  Feature ppp;
  for (;;){             // loop over all rows in file
    if (getline(&buf, &size, pipe) < 0) {
      break;
    }
    stringstream parse(chomp(buf));
    int rowtype;
    parse >> rowtype;
    if (rowtype == 1) {
      if (phase > 0) {
        foo.push_back(ppp);
      }
      ppp = Feature();          // start a new Feature
      parse >> elev_foot;
      parse >> ppp.has_tower;
      parse >> bldgs;
      parse >> ppp.ident;

      getline(parse, fullname);
// trim off leading spaces:
      while (fullname.length() && fullname[0] == ' ')
        fullname = fullname.substr(1);
      ppp.fullname = fullname;
      ppp.cat = apt;
      rwyAccum = GeoPoint(0,0);
      nRwyEnds = 0;
      phase = 1;
    } else if (rowtype == 10) { // occurs in fileversion 810
                        // inherited by 850 but used only for taxiways
      if (phase != 1) { // should never happen
        cerr << "Phase error; bad format in file '" << fn << "'" << endl;
        exit(1);
      }

      Runway rwy;
      parse >> lat_deg;
      parse >> lon_deg;
      parse >> rwylabel;
      if (rwylabel == "xxx") continue;   // indicates a taxiway
      parse >> heading;
      parse >> rwylength;
      parse >> displaced;
      parse >> stopway;
      parse >> wid_ft;          // width in feet
      parse >> lighting_code;
      parse >> rwy.sfc;
      GeoPoint center(lat_deg*deg, lon_deg*deg);
      GeoPoint firstend = center.displace(heading*deg, rwylength*foot/2.);
      rwyAccum += firstend;
      nRwyEnds++;
      GeoPoint otherend = center.displace(pi + heading*deg, rwylength*foot/2.);
      rwyAccum += otherend;
      nRwyEnds++;

      ppp.loc = rwyAccum;
      ppp.loc /= nRwyEnds;
      ppp.loc.elev = elev_foot*foot;
      rwy.end[0] = firstend;
      rwy.end[1] = otherend;
      rwy.width = wid_ft * foot;
      ppp.rwy.push_back(rwy);
    } else if (rowtype == 100) { // occurs in fileversion 850
      if (phase != 1) { // should never happen
        cerr << "Phase error; bad format in file '" << fn << "'" << endl;
        exit(1);
      }

      Runway rwy;
      parse >> rwy.width;           // in meters
      parse >> rwy.sfc;
      parse >> shoulderType;
      parse >> smoothness;
      parse >> centerLights;
      parse >> edgeLights;
      parse >> distanceSigns;

      for (int rend = 0; rend < 2; rend++){
        parse >> rwylabel;
        parse >> lat_deg;
        parse >> lon_deg;
        //...
        parse >> displaced;
        parse >> overrun;
        parse >> markingType;
        parse >> approachLights;
        parse >> tdzLights;
        parse >> REIL;
        rwyAccum += rwy.end[rend] = GeoPoint(lat_deg*deg, lon_deg*deg);
        nRwyEnds++;
        ppp.loc = rwyAccum;
        ppp.loc /= nRwyEnds;
        ppp.loc.elev = elev_foot*foot;
      }
      ppp.rwy.push_back(rwy);
      // that's all for this type-100 runway-record
    } else {
      // ignore other record types
    }
  } // end of file
  if (phase > 0) foo.push_back(ppp);
  fclose(pipe);
}

int prefix(const string shorter, const string longer){
  return shorter == longer.substr(0, shorter.length());
}

// returns the difference in the range -180 to +180.
double diffangle(const double a, const double b){
  double rslt(a-b);
  rslt = fmod(rslt, 360.);
  if (rslt > 180.) rslt -= 360.;
  return rslt;
}

void usage() {
  cerr <<
"Usage:  getafix [options] TARGET\n"
"Options include:\n"
"  Selectors:\n"
"       -range R        Show all fixes within a radius R.\n"
"       -count N        Show the nearest N fixes.\n"
"       -box x1 y1 x2 y2  Show all fixes in box with corners (x1,y1) (x2,y2)\n"
"                         [all specified in degrees]\n"
"  Restrictors:\n"
"       -azimuth A W    Limit radial search to direction A +- W [in degrees]\n"
"       -apt            Show only records of type 'apt'\n"
"       -fix            Show only records of type 'fix'\n"
"       -obstruction    Show only records of type 'obst'\n"
"       -duplicates     Include duplicate records in output\n"
"  Special features:\n"
"       -track          Ignore selectors;\n"
"                       just spit out the lat/lon of the targets\n"
"                       along with the true course and distance\n"
"                       (and reverse course) to (from) the next target.\n"
"       -kml            ????\n"
"       -paved          ????\n"
"  Settings and modes:\n"
"       -corridor hwid  In track mode: halfwidth of corridor,\n"
"                       i.e. margin on each side of centerline.\n"
"                       Default=4 is relevant to FAR 91.177.\n"
"       -meters         Output all distances in meters\n"
"                       (as opposed to some in nm).\n"
"       -near ident     Locate first target in this vicinity if there is ambiguity.\n"
"\n"
"Notes: *) On the cmd line, all distances are in nautical miles.\n"
"       *) On the cmd line, bearings are in degrees, CW from true north.\n"
"       *) Selectors are cumulative; e.g. specifying a range and a\n"
"         box gives you everything within that range plus everything\n"
"         within that box.\n"
"\n"
"On the command line, each TARGET or ident must already be in the database;\n"
"there is not yet any provision for entering lat/lon on the fly.\n"
"Your best bet is to edit ./fix2.dat and define suitable fixes.\n"
"\n"
;
}

// some globals ... mostly mode settings
  double az_center;
  double az_width(360.);
  double track_elevation(0);
  int elevation_relative(1);
  int track_mode(0);
  int kml_mode(0);
  int need_paved(0);
  double corr_hwid(4. * nm);
  int include_dupes(0);
  double max_range(0);
  unsigned int max_count(0);
  catter need_cat(feature);
  double output_unit(nm);
  string vicinity_name("");
  GeoPoint box1(9e99, 9e99), box2(9e99, 9e99);
  string outfile_prefix("some");

string protect(const string xxx){
  string rslt;
  for (string::const_iterator ptr = xxx.begin();
        ptr != xxx.end(); ptr++){
    char ch = *ptr;
    if (ch == '&') rslt += "&amp;";
    else rslt += ch;
  }
  return rslt;
}

string fixup(const string xxx){
  string rslt(xxx);
  for (string::iterator ptr = rslt.begin();
        ptr != rslt.end(); ptr++){
    if (*ptr == '$') *ptr = '"';
  }
  return rslt;
}

string kml_header(fixup(
"<?xml version=$1.0$ encoding=$UTF-8$?>\n"
"<kml xmlns=$http://www.opengis.net/kml/2.2$>\n"
"<Document>\n"
));

string kml_trailer(fixup(
"</Document>\n"
"</kml>\n"
));

string kml_styleblock(fixup(
"    <Style id=$fat-blue$>\n"
"        <LineStyle>\n"
"          <color>ffFF0000</color>\n"
"          <width>4</width>\n"
"        </LineStyle>\n"
"    </Style>\n"
"\n"

// track centerline
"    <Style id=$track$>\n"
"      <LineStyle>\n"
"        <color>ff0000FF</color>\n"
"        <width>4</width>\n"
"      </LineStyle>\n"
"    </Style>\n"
"\n"

// route corridors:
"    <Style id=$corridor$>\n"
"      <LineStyle>\n"
"        <color>000000FF</color>\n"
"        <width>0</width>\n"
"      </LineStyle>\n"
"      <PolyStyle>\n"
"        <color>800000FF</color>\n"
"      </PolyStyle>\n"
"    </Style>\n"
"\n"

// airport style includes runways and the
// circular point-marks to which the names attach:
"    <Style id=$poly$>\n"
"      <IconStyle>\n"
"         <color>ffFF4040</color>\n"
"         <Icon>\n"
"   <href>http://maps.google.com/mapfiles/kml/shapes/donut.png</href>\n"
"         </Icon>\n"
"         <scale>0.5</scale>\n"
"      </IconStyle>\n"
"      <LineStyle>\n"
"        <color>ffFF0000</color>\n"
"        <width>4</width>\n"
"      </LineStyle>\n"
"      <PolyStyle>\n"
"        <color>40FF0000</color>\n"
"      </PolyStyle>\n"
"    </Style>\n"
"\n"
"    <Style id=$triangle$>\n"
"      <IconStyle>\n"
"         <color>ff000000</color>\n"
"         <Icon>\n"
"   <href>http://maps.google.com/mapfiles/kml/shapes/triangle.png</href>\n"
"         </Icon>\n"
"         <scale>0.5</scale>\n"
"      </IconStyle>\n"
"      <LabelStyle>\n"
"        <color>ff000000</color>\n"
"        <colorMode>normal</colorMode>\n"
"        <scale>0.6</scale>\n"
"      </LabelStyle>\n"
"    </Style>\n"
"\n"
"    <Style id=$green-poly$>\n"
"      <LineStyle>\n"
"        <color>ff00FF00</color>\n"
"        <width>1</width>\n"
"      </LineStyle>\n"
"      <PolyStyle>\n"
"        <color>8000FF00</color>\n"
"      </PolyStyle>\n"
"    </Style>\n"
"\n"
));

class bad_thing: public std::exception{
  const char* msg;
  virtual const char* what() const throw() {
    return msg;
  }
public:
  bad_thing(const char* _msg)
  : msg(_msg) {}
};

// If indent is less than zero, don't do anything.
class Bracket{
public:
  ostream& xout;
  vector<string> trailer;

// Constructor

  Bracket(ostream& _xout) : xout(_xout), trailer(0) {}

  int operator() (const int indent,
        const string header, const string new_trailer)
  {
    trailer.push_back(string(indent, ' ') + new_trailer);
    xout << string(indent, ' ') << header;
    return indent;
  }

// nnl is the number of newlines to add to the trailer
  int operator() (const int indent, const string header, int nnl = 0)
  {
    string new_trailer = header;
    size_t where = new_trailer.find('<');
    if (where == string::npos) throw bad_thing("Bracket: missing '<'");
    new_trailer.insert(1+where, "/");
    new_trailer += string(nnl, '\n');
    trailer.push_back(string(indent, ' ') + new_trailer);
    xout << string(indent, ' ') << header;
    return indent;
  }


// Destructor:
  ~Bracket(){
    for (vector<string>::const_reverse_iterator sss = trailer.rbegin();
                sss != trailer.rend(); sss++)
      xout << *sss;
  }

};

typedef map<string,Bucket,CaseFold> MSB;

void secchi(ostream& xout, const GeoPoint center, const double phase){
// radius of secchi disk, in meters:
    const double secchi_rad(10.);

    int indent(2);
    Bracket gorp(xout);
    indent = 2 + gorp(indent, "<Polygon>\n");
    xout << string(indent, ' ')
       << "<altitudeMode>absolute</altitudeMode>" << endl;
    indent = 2 + gorp(indent, "<outerBoundaryIs>\n");
    indent = 2 + gorp(indent, "<LinearRing>\n");
    indent = 2 + gorp(indent, "<coordinates>\n");

    center.dump(xout);
    GeoPoint offcenter_EW(center), offcenter_NS(center);
    offcenter_EW.lonRad += 1.;        // add one radian east
    offcenter_NS.latRad += 1.;        // add one radian north
// one radian of longitude, measured in SI meters:
// roughly equal to the effective radius of curvature
// of the (non-great) circle of constant latitude
    double radcurve_EW = center.geoDistTo(offcenter_EW);
    double radcurve_NS = center.geoDistTo(offcenter_NS);
    for (int ii = 0; ii <= 90; ii += 10) {
      double theta = - ii*deg + pi/4. + phase;
      GeoPoint circum(center);
      circum.lonRad += secchi_rad*sin(theta) / radcurve_EW;
      circum.latRad += secchi_rad*cos(theta) / radcurve_NS;
      circum.elev   += secchi_rad;
      circum.dump(xout);
    }
    center.dump(xout);
}

class plain_handler_base {
public:
  int placemark_isopen;
  string placemark_trailer;
  string mystyle;

  plain_handler_base()
  : placemark_isopen(0), placemark_trailer(""), mystyle("????")
  {}

// the do_folder routine is not even virtual:
  void do_folder(const vector<Feature> incount, ostream& xout,
                const string foldername = "");
  virtual void setup(){};
  virtual void within_placemark(ostream& xout, const Feature&){};
  virtual void do_feature(ostream& xout, const Feature&){};
  virtual void end_placemark(ostream& xout){};
  virtual void do_runway(ostream& xout, const Feature&, const Runway&){};
};

class kml_handler_base : public plain_handler_base{
public:
  virtual void within_placemark(ostream& xout, const Feature&);
  virtual void do_feature(ostream& xout, const Feature&){};
  virtual void end_placemark(ostream& xout);
  virtual void do_runway(ostream& xout, const Feature&, const Runway&){};
};

class kml_handler_ascii : public plain_handler_base{
public:
  virtual void do_feature(ostream& xout, const Feature&);
};

class kml_handler_rectangle : public kml_handler_base{
public:
  virtual void do_runway(ostream& xout, const Feature&, const Runway&);
  virtual void setup();
};

class kml_handler_secchi : public kml_handler_base{
public:
  virtual void do_runway(ostream& xout, const Feature&, const Runway&);
  virtual void setup();
};

class kml_handler_centerline : public kml_handler_base{
public:
  virtual void do_runway(ostream& xout, const Feature&, const Runway&);
  virtual void setup();
};

class kml_handler_points : public kml_handler_base{
public:
  virtual void do_feature(ostream& xout, const Feature&);
  virtual void setup();
};

void kml_handler_ascii::do_feature(ostream& xout,
        const Feature& feat){
      int xwid(0);
      xout << fixed << setprecision(3)
          << right << setw(xwid) << feat.fwd/deg
          << "|" << setw(xwid) << feat.dist / output_unit
          << setprecision(9)
          << "|" << setw(xwid) << feat.loc.lonRad/deg
          << "|" << setw(xwid) << feat.loc.latRad/deg
          << "|" << setw(xwid) << left << feat.ident
          << setprecision(1);
          if (feat.loc.elev <= -9e99)
            xout << "|" << scientific << feat.loc.elev;
// output elevation in meters;
// may or may not be the smart thing to do:
          else xout << "|" << fixed << feat.loc.elev;
          xout << fixed
          << "|" << feat.has_tower
          << "|" << feat.fullname
          << "|"
          << endl;
}

void kml_handler_secchi::do_runway(ostream& xout,
        const Feature& feat, const Runway& rrr){
  within_placemark(xout, feat);
  for (int rend = 0; rend < 2; rend++){
    GeoPoint elevpoint;
    elevpoint = rrr.end[rend];
    elevpoint.elev = feat.loc.elev;
    secchi(xout, elevpoint, 0.);
    secchi(xout, elevpoint, pi);
  }
}

void kml_handler_rectangle::do_runway(ostream& xout,
        const Feature& feat, const Runway& rrr){
  within_placemark(xout, feat);
  double fwd, rev;
  double len_not_used = rrr.end[0].geoDistTo(rrr.end[1], fwd, rev);
  if (len_not_used) {}

  xout <<
"        <Polygon>\n"
"            <altitudeMode>clampToGround</altitudeMode>\n"
"            <outerBoundaryIs>\n"
"                <LinearRing>\n"
"                    <coordinates>\n";

  GeoPoint corner = rrr.end[0].displace(fwd-pi/2, rrr.width/2.);
  corner.dump(xout);
  rrr.end[1].displace(rev+pi/2, rrr.width/2.).dump(xout);
  rrr.end[1].displace(rev-pi/2, rrr.width/2.).dump(xout);
  rrr.end[0].displace(fwd+pi/2, rrr.width/2.).dump(xout);
  corner.dump(xout);

xout <<
"                    </coordinates>\n"
"                </LinearRing>\n"
"            </outerBoundaryIs>\n"
"        </Polygon>\n";

}

void kml_handler_centerline::do_runway(ostream& xout,
        const Feature& feat, const Runway& rrr){
  within_placemark(xout, feat);
  int indent = 6;
  Bracket gorp(xout);
  indent = 2 + gorp(indent, "<LineString>\n");
  xout << string(indent, ' ')
        << "<altitudeMode>clampToGround</altitudeMode>\n";
  indent = 2 + gorp(indent, "<coordinates>\n");
  rrr.end[0].dump(xout);
  rrr.end[1].dump(xout);
}

void kml_handler_points::do_feature(ostream& xout, const Feature& feat){
  within_placemark(xout, feat);
  xout << "  <Point>" << endl;
  xout << "  <altitudeMode>clampToGround</altitudeMode>" << endl;
  xout << "    <coordinates>" << endl;
  GeoPoint tmp(feat.loc);
  tmp.elev = 0;
  tmp.dump(xout);
  xout << "    </coordinates>" << endl;
  xout << "  </Point>" << endl;
}

string std_placemark_trailer("  </MultiGeometry>\n"
                             "  </Placemark>\n");

void kml_handler_base::within_placemark(ostream& xout, const Feature& feat){
    if (placemark_isopen) return;
    xout << "<Placemark>" << endl;
    xout << "  <name>" << feat.ident << "</name>" << endl;
    xout << "  <description>" << protect(feat.fullname)
            << "</description>" << endl;
#if 0
    if (need_cat == fix)
      xout << "  <styleUrl>#triangle</styleUrl>\n";
    else
      xout << "  <styleUrl>#poly</styleUrl>\n";
#endif
    xout << "  <styleUrl>#" << mystyle << "</styleUrl>\n";
    xout << "  <MultiGeometry>" << endl;
    placemark_isopen = 1;
    placemark_trailer = std_placemark_trailer;
}

// using green coloration for the rectangles
void kml_handler_secchi::setup() {
  mystyle = "green-poly";
}

// using triangles for the airports
void kml_handler_points::setup() {
  mystyle = "triangle";
}

void kml_handler_rectangle::setup() {
  mystyle = "poly";
}

void kml_handler_centerline::setup() {
  mystyle = "poly";             // may not be correct
}

void kml_handler_base::end_placemark(ostream& xout){
    if (!placemark_isopen) return;
    xout << placemark_trailer;
    placemark_isopen = 0;
}

void plain_handler_base::do_folder(const vector <Feature> grand,
        ostream& xout, const string foldername) {
  setup();
  placemark_isopen = 0;
  if (foldername.length()) {
     xout << "<Folder><name>" << foldername << "</name>" << endl;
     xout << kml_styleblock;
  }
  for (vector<Feature>::const_iterator feat = grand.begin();
      feat != grand.end(); feat++) {

    if (need_paved) {
      int is_paved(0);
      for (vector<Runway>::const_iterator rptr = feat->rwy.begin();
          rptr != feat->rwy.end(); rptr++) {
        int code = rptr->sfc;
        if (code == 1 || code == 2) {
          is_paved = 1;
          break;
        }
      }
      if (!is_paved) continue;
    }

    do_feature(xout, *feat);

// loop over all runways belonging to this Feature:
    for (vector<Runway>::const_iterator rptr = feat->rwy.begin();
        rptr != feat->rwy.end(); rptr++) {
      do_runway(xout, *feat, *rptr);
    }

    end_placemark(xout);

  } // all features in featureset
  if (foldername.length()) {
    xout << "</Folder>" << endl;
  }
}

///////////////////////////////
// the big output routine
//
// called by collect_and_show_features
//
void output_featureset(const vector<Feature>& grand, ostream& xout) {
  if (kml_mode){
    if (1) {
// Areas have to be in layer #1, for GRASS
// (as discussed in grass-intro.htm).
      kml_handler_rectangle().do_folder(grand, xout, "runways");
// Labels don't have to be in layer #1:
      kml_handler_points().do_folder(grand, xout,
                        catter_name[need_cat] + "-names");
      kml_handler_secchi().do_folder(grand, xout, "elevations");
    } else {
      kml_handler_centerline().do_folder(grand, xout);
    }
  } else {
    kml_handler_ascii().do_folder(grand, xout);
  }
}

void collect_and_show_features(vector<Feature>& foo,
        vector<Feature> target_features) {
  ofstream xout;
  string fname = outfile_prefix;
  if (fname.length()) fname += "_";
  fname += catter_name[need_cat];
  if (kml_mode) fname += ".kml";
  else fname += ".csv";
  cerr << "Using fname " << fname << endl;
  xout.open(fname.c_str());
  if (!xout.good()) {
    cerr << "Could not open output file '" << fname << "'" << endl;
    exit(1);
  }

  Bracket gorp(xout);
  if (kml_mode) gorp(0, kml_header, kml_trailer);
  Feature tarfeat;
  set<string> seen_it;
  vector<Feature> grand_featureset;
  int first_target(1);

// outer loop over all targets:
  for (vector<Feature>::const_iterator tp = target_features.begin();
        tp != target_features.end(); tp++) {
    tarfeat = *tp;

// Preliminary inner loop: calculate distance and azimuth for each point,
// and cache the results.
// This is an optimization, based on the idea that calculating
// the distance is expensive, and we are going to need the distance
// more than once.
// Indeed std::set is going to use the distance again and again
// when we're not looking.
    for (vector<Feature>::iterator ptr = foo.begin();
          ptr != foo.end(); ptr++) {
      double rev;
      ptr->dist = tarfeat.loc.geoDistTo(ptr->loc, ptr->fwd, rev);
    }

    vector<Feature> inrange;
// Using std::set accomplishes two things:
// It keeps the points sorted, and eliminates duplicates.
    set<Feature,distcomp> incount(7);
// inner loop over all fixes in the world:
    for (vector<Feature>::const_iterator ptr = foo.begin();
          ptr != foo.end(); ptr++) {
      string ident(ptr->ident);
      double dist(ptr->dist);
      if (need_cat != feature && ptr->cat != need_cat) continue;
// range-based selector:
      if (dist <= max_range){
        if ( (az_width >= 180.)  // don't even need to compute azimuth
            || (dist < 1.)   // azimuth is ill-determined when range is small.
            || fabs(diffangle(ptr->fwd, az_center)) <= az_width) {
          inrange.push_back(*ptr);
        }
      }
// count-based selector:
      if (max_count > 0) {
        if (incount.size() < max_count
             || dist < incount.rbegin()->dist) {
          incount.insert(*ptr);
        }
        if (incount.size() > max_count) {
          incount.erase((++incount.rbegin()).base());
        }
      }
// box-based selector
// FIXME:  should be independent of target
// ... but if we don't have a target, what would
// we use for the distance and azimuth fields
// that appear in the output?
// For now, latch onto first target:
      if (first_target) {
        if (     ptr->loc.lonRad >= box1.lonRad
              && ptr->loc.latRad >= box1.latRad
              && ptr->loc.lonRad <= box2.lonRad
              && ptr->loc.latRad <= box2.latRad ) {
          inrange.push_back(*ptr);
        }
      }
    } // inner loop over all fixes
    first_target = 0;

#if 1
    cerr << "For target " << tarfeat.ident
         << " found " << inrange.size()
         << " plus " << incount.size() << endl;
#endif

// combine results from various selectors,
// in such a way that they become sorted by distance:
    for (vector<Feature>::const_iterator xxx = inrange.begin();
        xxx != inrange.end(); xxx++){
      incount.insert(*xxx);
    }
// and accumulate the results:
    for (set<Feature>::const_iterator xxx = incount.begin();
        xxx != incount.end(); xxx++){
      if ((!include_dupes)
        && seen_it.find(xxx->ident) != seen_it.end()) continue;
      seen_it.insert(xxx->ident);
      grand_featureset.push_back(*xxx);
    }
  } // target points
  output_featureset(grand_featureset, xout);
}

// Plot the track centerline.
// (Note this does not include the corridor;
// use show_corr for that.)
void show_track_grass_fmt(vector<Feature> target_features) {

  cout << "L"
        << " " << fixed << setprecision(0) << target_features.size()
        << " " << 0     // number of categories
        << endl;
  Feature here, there;
  for (vector<Feature>::const_iterator tp = target_features.begin();
        tp != target_features.end(); )
  {
    here = *tp;
    tp++;
    if (tp == target_features.end()) there = here;
    else there = *tp;

    double fwd, rev, dist;
    dist = here.loc.geoDistTo(there.loc, fwd, rev);
    cout << fixed << setprecision(9)
               << here.loc.lonRad/deg
        << " " << here.loc.latRad/deg
        << " " << here.ident
        << setprecision(3)
        << " " << fwd/deg               /* not used by grass */
        << " " << dist / output_unit
        << " " << rev/deg
        << endl;
  }
}


// draw a half circle, using a total of Nseg segments, i.e.
// i.e. Nseg+1 points
//
// Special case:  if Nseg==0, just output a single point,
// namely the first point that would have belonged to the
// half circle.  This is useful for closing contours.
void endcap(ostream& xout, const Feature& here,
        const double facing, const double flip, const int Nseg) {
  double denom(Nseg);
  if (Nseg == 0) denom = 1;         // special case
  for (int ii = 0; ii <= Nseg; ii++) {
    double theta = facing - pi/2. - flip * pi * ii / denom;

// A point walking around the circle:
    GeoPoint walker = here.loc.displace(theta, corr_hwid);
    walker.dump(xout);
  }
}

class Segment_base{
public:
  string myname;
  int polygonal;
  void all_segments(vector<Feature> target_features);
  virtual void show_segment_coords(ostream& xout,
        const Feature here, const Feature there, int routine){};
  virtual void setup(){}
  virtual void show_segment_kml(
        const int indent, ostream& xout,
        const Feature here, const Feature there, int routine);
};

class Tracker : public Segment_base{
  virtual void show_segment_coords(ostream& xout,
                const Feature here, const Feature there, int routine){

    if (!routine) return;
    here.loc.dump(xout);
    there.loc.dump(xout);
  }
  virtual void setup() {
    myname = "track";
    polygonal = 0;
  }
};

class Corridor : public Segment_base{
  virtual void show_segment_coords(ostream& xout,
                const Feature here, const Feature there, int routine);
  virtual void setup(){
    myname = "corridor";
    polygonal = 1;
  }
};

void Corridor :: show_segment_coords(ostream& xout,
                const Feature here, const Feature there,
                const int routine){
    int Nseg(18);
    if (!routine){
      endcap(xout, here, 0, 2., 2*Nseg);
    } else {
      double dist, fwd, rev;
      dist = here.loc.geoDistTo(there.loc, fwd, rev);
      if (dist) {}      // suppress compiler warning`
      endcap(xout, here,  fwd, -1., Nseg);
      endcap(xout, there, rev, -1., Nseg);
      endcap(xout, here,  fwd, -1., 0);
    }
}

void Segment_base:: show_segment_kml(
        const int _indent, ostream& xout,
        const Feature here, const Feature there, int routine){

  int indent(_indent);
  Bracket gorp(xout);
  indent = 2 + gorp(indent, "<Placemark>", 1);
// We are better off without  a description, so we don't get balloons:
//xxx  xout << "<description>corridor</description>" << endl;
// but we ought to provide some sort of name:
  xout << "<name>" << here.ident;
  if (routine) xout << "_" << there.ident;
  xout << "</name>" << endl;
//  Changing displayMode doesn't fix the clickable problem:
//xxx  xout << "<displayMode>hide</displayMode>" << endl;
// xout << "<clickable>0</clickable>" << endl;

  xout << string(indent, ' ')
        << "<styleUrl>#" << myname << "</styleUrl>\n";

// not needed at the moment:
//xxxx  indent = 2 + gorp(indent, "<MultiGeometry>\n");

  double fwd, rev, dist;
  dist = here.loc.geoDistTo(there.loc, fwd, rev);
  if (dist) {}          // suppress compiler warning

  Bracket bork(xout);
  if (polygonal) indent = 2 + bork(indent, "<Polygon>\n");
  else           indent = 2 + bork(indent, "<LineString>\n");
  if (elevation_relative) xout << string(indent, ' ')
        << "<altitudeMode>clampToGround</altitudeMode>\n";
  else xout << string(indent, ' ')
        << "<altitudeMode>absolute</altitudeMode>\n";

  if (polygonal) {
    indent = 2 + bork(indent, "<outerBoundaryIs>\n");
    indent = 2 + bork(indent, "<LinearRing>\n");
  }
  indent = 2 + bork(indent, "<coordinates>\n");
  show_segment_coords(xout, here, there, routine);
}


// Plot the track centerline or corridor
//
// The first part of this code runs parallel to collect_and_show_features.
void Segment_base::all_segments(vector<Feature> target_features) {

  setup();
  ofstream xout;
  string fname = outfile_prefix;
  if (fname.length()) fname += "_";
  fname += myname;
  fname += ".kml";
  cerr << "Using fname " << fname << endl;
  xout.open(fname.c_str());
  if (!xout.good()) {
    cerr << "Could not open output file '" << fname << "'" << endl;
    exit(1);
  }

  int indent = 0;
  Bracket gorp(xout);
  indent = 2 + gorp(indent, kml_header, kml_trailer);
  indent = 2 + gorp(indent, "<Folder>", 1);
  xout << "<name>" << myname << "</name>" << endl;
  xout << kml_styleblock;

// special case: first target:
  Feature here, there;
  vector<Feature>::const_iterator tp = target_features.begin();
  there = *tp++;         // initial point
  there.loc.elev = track_elevation;
  here = there;
  show_segment_kml(indent, xout, here, there, 0);

// loop over all routine segments ...
// where a segment lies between two targets.
  for (; tp != target_features.end(); tp++) {
    here = there;       // new start-point = old end-point
    there = *tp;
    there.loc.elev = track_elevation;
// the segment:
    show_segment_kml(indent, xout, here, there, 1);
// the disk between segments:
    show_segment_kml(indent, xout, there, there, 0);
  }
}

void suffix(const string arg, string& head, string& sfx){
  int ii = arg.length();
  for (; ii > 0;){
    ii--;
    char ch(arg[ii]);
    if (ch == ' ') continue;
    if (isalpha(ch)) continue;
    break;
  }
  head = arg.substr(0, ii+1);
  sfx = arg.substr(ii+1);
}

int main(int _argc, const char* _argv[]) {
  init_catter_names();
  int argc(_argc);
  const char** argv(_argv);
  string progname(*argv++); argc--;
  vector<string> target_names(0);

  for (;argc > 0;) {
    string arg(*argv++); argc--;
    if (arg.substr(0,2) == "--") arg = arg.substr(1);
    if (prefix(arg, "-help")) {
      usage();
      exit(0);
    } else if (prefix(arg, "-azimuth")) {
      if (!argc) {
        cerr << "Option -azimuth requires two arguments" << endl;
        return 1;
      }
      az_center = atof(*argv++); argc--;

      if (!argc) {
        cerr << "Option -azimuth requires two arguments." << endl;
        return 1;
      }
      az_width = atof(*argv++); argc--;
    } else if (prefix(arg, "-range")) {
      if (!argc) {
        cerr << "Option -range requires an argument" << endl;
        return 1;
      }
      max_range = atof(*argv++) * nm; argc--;
    } else if (prefix(arg, "-count")) {
      if (!argc) {
        cerr << "Option -count requires an argument" << endl;
        return 1;
      }
      max_count = atoi(*argv++); argc--;
    } else if (prefix(arg, "-box")) {
      if (!argc) {
        cerr << "Option -box requires four arguments" << endl;
        return 1;
      }
      box1.lonRad = atof(*argv++)*deg; argc--;

      if (!argc) {
        cerr << "Option -box requires four arguments" << endl;
        return 1;
      }
      box1.latRad = atof(*argv++)*deg; argc--;

      if (!argc) {
        cerr << "Option -box requires four arguments" << endl;
        return 1;
      }
      box2.lonRad = atof(*argv++)*deg; argc--;

      if (!argc) {
        cerr << "Option -box requires four arguments" << endl;
        return 1;
      }
      box2.latRad = atof(*argv++)*deg; argc--;
#if 0
      dumpPoint(&cerr, box1);
      dumpPoint(&cerr, box2);
      double fwd, rev;
      double dist = box1.geoDistTo(box2, fwd, rev);
      cerr << "Diameter of box: " << dist/nm << " nm"
          << "  on heading " << fwd/deg
          << endl;
#endif
    } else if (prefix(arg, "-corridor")) {
      if (!argc) {
        cerr << "Option -corridor requires an argument" << endl;
        return 1;
      }
      corr_hwid = atof(*argv++) * nm; argc--;
    } else if (prefix(arg, "-elevation")) {
      if (!argc) {
        cerr << "Option -elevation requires an argument" << endl;
        return 1;
      }
      string head, sfx;
      suffix(*argv++, head, sfx); argc--;
      track_elevation = atof(head.c_str());
      if (sfx == "ft") track_elevation *= foot;
      else if (sfx == "m") {}
      else {
        cerr << "-elevation requires units: ft or m" << endl;
        exit(1);
      }
      elevation_relative = 0;           // elevation is now absolute
    } else if (prefix(arg, "-near")) {
      if (!argc) {
        cerr << "Option -near requires an argument" << endl;
        return 1;
      }
      vicinity_name = *argv++; argc--;
    } else if (prefix(arg, "-prefix")) {
      if (!argc) {
        cerr << "Option -prefix requires an argument" << endl;
        return 1;
      }
      outfile_prefix = *argv++; argc--;
    } else if (prefix(arg, "-meters")) {
      output_unit = 1.0;
    } else if (prefix(arg, "-apt") || prefix(arg, "-airport")) {
      need_cat = apt;
    } else if (prefix(arg, "-fix")) {
      need_cat = fix;
    } else if (prefix(arg, "-obstruction")) {
      need_cat = obst;
    } else if (prefix(arg, "-include")) {
      include_dupes++;
    } else if (prefix(arg, "-track")) {
      track_mode++;
    } else if (prefix(arg, "-kml")) {
      kml_mode++;
    } else if (prefix(arg, "-paved")) {
      need_paved++;
    } else if (arg[0] == '-') {
      cerr << "Unrecognized option: '" << arg << "'" << endl;
      return 1;
    } else {
      target_names.push_back(arg);
    }
  }

// FIXME : in box mode, only reason to require at least one target
// is not at all a good reason.
// FWIW, in non-box mode, no targets would result in no output.
  if (target_names.size() == 0) {
    cerr << "Need at least one target;  try\n";
    cerr << "   getafix -help\n";
    exit(1);
  }

// there were 140,000 known fixes last time I looked
  vector<Feature> foo;
  foo.reserve(250000);
  scan_fixes("fix.dat", 0, foo);
  scan_fixes("fix2.dat",   1, foo);
  scan_fixes("obst.dat",   1, foo, obst);
  scan_apts("apt.dat", 0, foo);
  scan_apts("apt2.dat", 1, foo);

//////////////////////////////////////////////////////////////////////
// all the data has been read in.

// Build the database, to facilitate looking things up by name:
  MSB dbase;
  for (vector<Feature>::const_iterator ptr = foo.begin();
      ptr != foo.end(); ptr++) {
    dbase[ptr->ident].push_back(*ptr);
  }

  Feature vicinity;
  if (vicinity_name.length()) {
    Bucket bbb(dbase[vicinity_name]);
    if (bbb.size() == 0) {
      cerr << "Vicinity fix '" << vicinity_name << "' not found." << endl;
      // continue, with bated breath
    }
    else if (bbb.size() > 1) {
      cerr << "Vicinity fix '" << vicinity_name << "' is "
        << bbb.size() << "-fold ambiguous"
        << endl;
      // also continue, with bated breath
    } else {
      vicinity = bbb[0];
    }
  }

  vector<Feature> target_features;
  for (vector<string>::const_iterator nnn = target_names.begin();
        nnn != target_names.end(); nnn++){
    string name = *nnn;
    Feature tarfeat = dbase[name].nearest_to(vicinity);
    if (tarfeat.cat == phantom) {
      cerr << "Fix '" << name << "' not found." << endl;
      continue;
    }
    if (tarfeat.cat == ambiguous) {
      cerr << "Fix '" << name << "' is "
        << tarfeat.has_tower << "-fold ambiguous"
        << endl;
      continue;
    }
    vicinity = tarfeat;
    target_features.push_back(tarfeat);
  }

  if (track_mode){
    Tracker().all_segments(target_features);
    Corridor().all_segments(target_features);
  } else {
// non-track mode ... find stuff in box and/or near targets
// and then send it to the output files:
    collect_and_show_features(foo, target_features);
  }
}  // end of main

getafix.c  ‍  ‍ (DL)

6.5 Hand Work

After the script has created all the map layers, we still need to tell the system which ones we want to display. I haven’t figured out how to automate this yet.

Go to the“GIS Manager” window and click the “Add raster layer button” (for raster data, such as was input via r.in.gdal) or the “Add vector layer” button (for vector data, such as was input via v.in.ascii).

Then tell it what map you want to use. The various options are not very self-explanatory but eventually you can figure them out. Beware that some of the more advanced options don’t work.

Whenever you change an option, if you want to see the effect you need to go to the “Map Display” window and push the redraw button, the second button from the left.

Important: If you have an opaque layer, such as any sort of elevation map or shaded relief map, make sure any vectors you draw are layered on top of the opaque map; otherwise they won’t be visible and you won’t know what happened to them. Specifically, in the “GIS Manager” window, the order in which the maps are listed is important. The ones higher in the list get drawn on top of the ones lower in the list.

Choosing the font to be used for labels on the graph apparently needs to be done by hand, every time you start the GUI. Evidently it is not part of the configuration stored in the .grc file.

7 Bugs in the GRASS

Here are only a small fraction of the bugs I’ve noticed:

8 Miscellaneous Notes

You can restart the GUI:
It restarts quickly, so this is a good way to work around GUI bugs:

: +++;  gis.m dmrc=LCC.grc

 1a) d.where seems to be grossly broken.  Runs amok.
 1b) Calling v.what or d.what.vect from the "run"
      window generates no events.
 2) v.what seems to sorta work *if* called by the button
   on the gui, but then only applies to one map at a time.

I wanted to be able to query the map, to identify features et
 cetera.  I've found the following dirty workaround: On the ``Map
 Display'' window, click the ``NVIZ flythrough path'' button, and then
 push the ``mouse'' button, then click on the map.  Go back to the
 ``NVIZ flythrough path'' popup and copy the numbers from the window
 next to the ``mouse'' button.

 This is useful because I can't get any of the mouse-related ``d.''
 modules to work.

List available database files (aka map layers and mapsets) in the
format "type/layer@MAPSET"
: +++; g.mlist -t -m type=all

less convenient format:
: ???;  g.list  type=rast | cat

find name of current location
: +++;  g.gisenv get=LOCATION_NAME

dump *all* gis env variables:
: +++;  g.gisenv --q

Find out about current "region" i.e. clipping region.
If this is not right lots of things will fail insidiously!
: +++;  g.region -p

Find out about current projection:
: +++;  g.proj -p

Find out about projection used by file:
: +++;  g.proj -p -d georef=N32W111.tif

Read in a file ... and create location accordingly ...
assuming you have put bounding-box and projection info into the file:
: +++; r.in.gdal input=N32W111.tif output=m111 --overwrite location=z111

Read it in ... without messing with the location:
: +++; r.in.gdal input=N32W111.tif output=m111 --overwrite

Note that if location already exists, you can switch
to it using g.mapset.

To find out what columns a map layer has
: +++;  v.info -c map=$foo layer=$bar

There does not appear to be a corresponding "r.info -c"
for raster maps.

Dump all the "attributes" of a vector map layer, i.e.  all the rows
and columns of data including e.g. the names of the data points (but
excluding the latitude and longitude, which are somehow not considered
"attributes")
: +++;  db.select table=$foo

Again, there does not appear to be any corresponding way
of looking at raster files.

If the map "$foo" has layers, then the table names
will have suffixes:
: +++; db.select table=${foo}_1
: +++; db.select table=${foo}_2              # et cetera

When reading .kml files into grass, kml "folders" produce grass
"layers".

Grass apparently won't plot an area unless the boundary is /closed/.
So at the end, repeat the first point.  In contrast, almost every
other graphics program in the world will plot the interior of a
rectangle given three of the four sides.

Grass apparently won't plot the labels of areas, just points.  So if
you want to label areas, you need to provide points (as well as
areas) in your data file.

This is tricky, because the gis.m GUI won't let you control the fill
color of areas for any layer other than layer #1.  So put the areas
in layer #1, and put the labels in layer #2.  (The GUI gives you a
choice as to which layer the labels come from, but no corresponding
choice for the areas.)

This gets even trickier if you have multiple areas that you might want
to fill (runways, route corridor, et cetera).  For google-earth it is
convenient to put them in multiple folders in the same document, but
for grass AFAICT you are stuck putting them in separate files.

You can display an image on the screen without fussing with the gis.m
GUI:

: +++; d.mon stop=x0      # kill old x0, if any, so it forgets the region
: +++; g.region rast=$tile -p      # set region BEFORE starting x0
: +++; d.mon start=x0
: +++; d.rast $tile

Similarly, you can create an image on disk without fussing with the
gis.m GUI:

: +++; d.mon start=PNG
: +++; d.vect map=$mapname@PERMANENT color=0:0:0 lcolor=255:0:0 \
        fcolor=255:0:0 display=shape,attr type=point,line,boundary,area \
        icon=basic/circle size=30 layer=2 attrcol=Descriptio \
        lsize=15 xref=left yref=center llayer=2
: +++; d.mon stop=PNG
: +++; display map.png

Grass stores "locations" in subdirectories in the grassdata directory
... and indeed it thinks every subdirectory is a location, which gets
ugly in some cases, e.g. if you have a set of ESRI data that takes the
form of a directory full of files.

To get a list of /possible/ names of locations that
have already been defined, try (in the grassdata directory)
: +++;  ls -d */PERMANENT
Switch to a different location:
: +++; g.mapset mapset=PERMANENT location=z111
Switch back (if ever needed):
: +++; g.mapset mapset=PERMANENT location=UTM_4326_SRTM

find out basic stuff about some map layer:
: +++; r.info m111    # (for raster map)
: +++; v.info apts    # (for vector map)

r.in.gdal allows a user to create a (binary) GRASS raster map layer,
or imagery group, from any GDAL supported raster file format, with an
optional title. The imported file may also be optionally used to
create a new location.

manpage: file:///usr/lib/grass64/docs/html/r.in.gdal.html
formats: http://www.gdal.org/formats_list.html

#Spatial Reference Set
#for local Lambert Conformal Conic projection
#covering the whole bounding box (i.e. all the tiles)

: +++; /usr/local/gdal/bin/gdalwarp -of VRT -s_srs EPSG:4326 -t_srs \
 '+proj=lcc +lat_1=33n +lat_2=45n +lon_0=97w' mcd12_sep2109.asc \
 lambert_mcd12_sep2109.vrt

#Reference:
#http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg06041.html

9 More Notes

These are largely notes made before writing the script. Many of the ideas noted here were incorporated into the script. Some of them have since been changed. In particular, a bevy of small but annoying changes changes were necessary in order to convert from Mercator to Lambert projection.

## Some useful grass commands:

 TILE=N37E141  # mark Fukushima #1 and #2
 echo -e '141.015|37.233|"F1"\n141.025556|37.316389|"F2"' | \
 v.in.ascii output=marks --overwrite

# Some typical data fetch commands:
 TILE=N32W112
 wget -c http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/$TILE.hgt.zip
 unzip $TILE.hgt.zip

 TILE=N37E141
 wget -c http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/$TILE.hgt.zip
 unzip $TILE.hgt.zip

# read in the data directly, in UTM projection,
# without bothering to reproject it using gdal:
 r.in.srtm $TILE.hgt --overwrite

r.info map=$TILE

# One way to set the region;
# not so good if you have multiple tiles:
  g.region rast=$TILE -p
#
## works: g.region rast=$TILE -p
## fails: d.rast map=$TILE catlist=915-9999
## manpage:  file:///usr/lib/grass64/docs/html/d.mon.html

## Kludge: if the map window is cockeyed,
## set the region the way you want it,
## then delete the map window;
##   when it gets recreated it will be focused on the right region

# create shaded relief map:
r.shaded.relief map=$TILE  --overwrite
r.colors map=$TILE.shade color=gray  # probably the default

# Note: color=rainbow is (for now) favored for the elevation data
#  ... as opposed to the shaded relief.
#  haxby and haxby -e are both OK
#               don't use red or magenta
#  elevation is OK, kinda orange, but so are the FAA sectionals
#               doesn't use red or magenta
#  terrain is OK, kinda green, somewhat subdued impression
#               doesn't use red or magenta
#  etopo2 == same as terrain

#  byg = unusual, could be useful, takes some getting used to
#               doesn't use red or magenta

#  corine = boring, mostly white
#  srtm = boring, brown
#  ramp = too stark
#  bcyr = OK
#  rainbow -e and bcyr -e have much more yellow and red


# Get data from the map, assuming it uses lat/lon
# i.e. not Lambert:
 r.what input=$TILE east_north=-110.791,32.44 # Mount Lemmon
## result:
##  -110.791|32.44||2778

## read a KML file:
   g.mapset mapset=PERMANENT location=UTM_4326_SRTM
   v.in.ogr -z --overwrite  dsn=foo.kml output=demo
   g.mapset mapset=PERMANENT location=z111
# get a list of maps that can be pulled over:
   v.proj location=UTM_4326_SRTM -l
   v.proj location=UTM_4326_SRTM input=demo --overwrite

## Note that r.what can read from standard input
## and claims that the standard field-separator is '|',
## but in fact it insists on fs=space, and it cannot be changed.

## cat apt.csv | v.in.ascii x=3 y=4 fs=, output=apts --overwrite
## v.info map=apts

## g.region rast=$TILE1,$TILE2 -p
## cat track.asc | v.in.ascii -n out=track format=standard --overwrite

## cat fix.csv | v.in.ascii fs=, output=fixes x=3 y=4 \
## columns="range double precision, az double precision, x double precision, y double precision, ident varchar(20)"
## v.info -c map=fixes

## g.remove rast=some_old_thing

## TILE=N32W111
## g.region rast=$TILE -p
## r.contour --overwrite input=$TILE output=${TILE}_con3 levels=914,2743
## r.contour --overwrite input=$TILE output=${TILE}_con6 levels=1829

##** You can waste huuuuge amounts of time if the region isn't set right:
##** g.region rast=$TILE -p
##** g.region rast=$TILE1,$TILE2 -p

## Other options, not necessarily recommended:
## wget ??? http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/$TILE.hgt.zip

## r.colors $TILE.shade rast=shaded112
## r.shaded.relief map=$TILE shadedmap=$TILE.shade --overwrite

## d.everything fails:  allegedly no monitor:
##//  d.rast map=$TILE catlist=915-9999
##//  d.font font="DejaVuSans-Bold"

## Works from the type-in window of the "Output - GIS.m" popup.
##  d.font path=/usr/share/fonts/truetype/ttf-dejavu/DejaVuSans-Bold.ttf

##  Documentation
##    file:///usr/lib/grass64/docs/html/d.mon.html

## converting from SRTM (which is UTM 4326) to Lambert:
## http://www.osgeo.org/pipermail/gdal-dev/2010-March/023793.html

10 Yet More Notes

// First MUST set up datum and ellipsoid; see
// http://grass.osgeo.org/newsletter/GRASSNews_vol3.pdf
//
// If this is a first-time startup in an empty directory,
// push the EPSG button with code 4326

// http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/N32W111.hgt.zip
// Bigger, not necessarily better:
// http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N32W112.hgt.zip


   TILE=N37E141  # Fukushima #1
   echo -e '141.015|37.233|"F1"\n141.025556|37.316389|"F2"' | v.in.ascii output=marks --overwrite

   TILE=N32W112
   wget -c http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/$TILE.hgt.zip
   wget -c http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/$TILE.hgt.zip
   unzip $TILE.hgt.zip
   r.in.srtm $TILE.hgt --overwrite
## Color table for raster map <N32W111> set to 'srtm'
   r.info map=$TILE
   g.region rast=$TILE -p
// --> delete the map window if it is cockeyed;
//   when it gets recreated it will be focused on the right region
   r.shaded.relief map=$TILE  --overwrite
   r.colors map=$TILE.shade color=gray

   g.region rast=$TILE -p
// fails: d.rast map=$TILE catlist=915-9999
// file:///usr/lib/grass64/docs/html/d.mon.html

// ...
   r.what input=$TILE east_north=-110.791,32.44 # Mount Lemmon
## -110.791|32.44||2778

  cat apt.csv | v.in.ascii x=3 y=4 fs=, output=apts --overwrite
  v.info map=apts

 g.region rast=$TILE1,$TILE2 -p
 cat track.asc | v.in.ascii -n out=track format=standard --overwrite

 cat fix.csv | v.in.ascii fs=, output=fixes x=3 y=4 columns="range double precision, az double precision, x double precision, y double precision, ident varchar(20)" --overwrite
 v.info -c map=fixes

 g.remove rast=some_old_thing

 TILE=N32W111
 g.region rast=$TILE -p
 r.contour --overwrite input=$TILE output=${TILE}_con3 levels=914,2743
 r.contour --overwrite input=$TILE output=${TILE}_con6 levels=1829

//** You can waste huuuuge amounts of time if the region isn't set right:
   g.region rast=$TILE -p
   g.region rast=$TILE1,$TILE2 -p

// Other options, not necessarily recommended:
???   wget http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/$TILE.hgt.zip

// r.colors $TILE.shade rast=shaded112
// r.shaded.relief map=$TILE shadedmap=$TILE.shade --overwrite

// d.everything fails from GRASS shell window:  allegedly no monitor:
???       d.rast map=$TILE catlist=915-9999
???       d.font font="DejaVuSans-Bold"

// Works from the type-in window of the "Output - GIS.m" popup.
    d.font path=/usr/share/fonts/truetype/ttf-dejavu/DejaVuSans-Bold.ttf

//  Documentation
//    file:///usr/lib/grass64/docs/html/d.mon.html

// converting from SRTM (which is UTM 4326) to Lambert:
// http://www.osgeo.org/pipermail/gdal-dev/2010-March/023793.html

11 References

1.
“Import of multiple SRTM tiles in one step”
http://grass.osgeo.org/wiki/HOWTO_import_SRTM_elevation_data#Import_of_multiple_SRTM_tiles_in_one_step
2.
Markus Neteler, “SRTM and VMAP0 data in OGR and GRASS SRTM data acquisition”
(GRASS-News newsletter, Volume 3, June 2005)
http://grass.osgeo.org/newsletter/GRASSNews_vol3.pdf
3.
Lambert Conic Conformal Projection
http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html
4.
“GRASS Tutorial”
http://www.ing.unitn.it/~grass/docs/tutorial_62_en/index.html
5.
Sectional Aeronautical Charts – Online
www.av8n.com/physics/sectional.htm
6.
Mapnik – free toolkit for developing mapping applications.
http://mapnik.org/faq/
7.
“Shuttle Radar Topography Mission’ (SRTM)
http://www2.jpl.nasa.gov/srtm/
8.
“X-Plane Airport & Navigation Data” http://data.x-plane.com/get_data.html
[Contents]
Copyright © 2011 jsd