[Developer]

Enhance image overlay maps in Google Earth with altitude attributes


Use Perl and the Imager module to enhance mapping applications by extracting and applying height information based on color to reveal a third dimension of data, showing more information in the same space.

Recent applications have greatly increased the ease of development and ubiquity of 2-D maps. Tools like Microsoft® Live Search Maps and Google Maps offer a wealth of tools for enhancing these single-plane maps, but often ignore altitude as the third dimension of information. This article presents tools and code to allow you to extract height information based on pixel colors, and apply that height information across the mapping context. The end result is a third dimension of data, showing more information in the same space and opening up new methods of visualization for your map users.
This article demonstrates using the Imager Perl module to read and process images to create per-pixel level altitude slices of source images. Combine these image slices with Keyhole Markup Language (KML), and Google Earth will render an approximate volumetric visualization of the new data set.
Hardware and software requirements
You'll need a modern machine to run Google Earth, including 3-D acceleration for usable results. Processing large images with Imager can be time-consuming, so a fast processor and lots of RAM is useful.
In addition to Google Earth, you'll need the Imager Perl module (see Resources), as well as Perl itself. Note that all of the code presented here is cross-platform and should run on any platform Google Earth and Perl will run on.


Back to top

General approach used in this article
To demonstrate the concepts described in the introduction, NOAA weather radar data for the contiguous United States is used. The associated color gradients are determined for each altitude, height data for each pixel is extracted, then the secondary image is sliced to produce the computed altitude images. Figure 1 shows an example of what this can look like.

Figure 1. Example 2-D and 3-D precipitation visuals


In the example shown above, the top image shows the standard 2-D echo returns for water in the atmosphere. The bottom image shows the same echo returns with their altitude mapped according to a second image that displays the cloud-top heights. As viewed from miles away and near ground level, the Google Earth image overlay rendering provides a much broader range of information to the viewer.
Although focusing on radar return data for "Echo Tops" and "Merged Reflectivity Composite" maps, these techniques are applicable to any set of planar images.
Extracting altitude information based on color
NOAA provides outstanding resources for acquiring its imagery products, but not for the data itself. See Resources for links to the KML and imagery files. In this article, the Echo Tops image will provide the height data, but will not produce any visualized pixels. All of the viewable imagery will come from the "Merged Reflectivity Composite" map. Figure 2 shows an example of cloud-top visualization and the associated height data from the key.

Figure 2. Example Echo Top image and key with colors


The amorphous blob in the bottom left is pixel data indicating the cloud-top heights from a selected portion of the full image. Height information in thousands of feet is presented in the top portion of the image, along with magnified sections of the color bar. Note how multiple colors make up the graduated color bar and that their index values vary in a relatively inconsistent manner. This inconsistency and lack of an easily processed table to translate color to height is a relatively common issue that must be dealt with when extracting data from images. This is most common where the original data used in generating the image is unavailable.
Creating an average color value for each "color area" is performed by first extracting each color area as shown in the bottom right of the image. These color chunks are extracted from the Echo Tops key image and saved as a series of files in the color.tiles/ directory. If your selected data does not have a single color for each level of altitude, you'll need to mimic the process. Make sure you save the files with filenames that will be ordered correctly by the ls -1 color.tiles/ command — for example: c00.png, c01.png, etc.


Back to top

Altitude mapping with Perl and Imager
With the color tiles created, the first step in the program is to read these color tiles and develop an average color value for each tile. Create a file called altitudeMapper.pl and insert the contents of Listing 1.

Listing 1. altitudeMapper.pl header, color tiles read
#!/usr/bin/perl -w
# altitudeMapper.pl - create slices of images based on computed height from
# color data in separate image
use strict;
use Imager;

die "specify altitude image, slicee image " unless ( @ARGV == 2 );
my( $fNameAlt, $fNameSli ) = @ARGV;

my %cls = (); # color averages
my %hgt = (); # height data for each pixel

######## compute the average color for each tile ########
for my $fName ( `ls -1 color.tiles/` )
{
chomp($fName);
my $tileImage = Imager->new;
$tileImage->read(file=>"color.tiles/$fName") or die "cannot open $fName: ";

my( $rAvg, $gAvg, $bAvg ) = 0;

for( my $colN = 0; $colN <= ($tileImage->getwidth -1); $colN++ )
{
for( my $rowN = 0; $rowN <= ($tileImage->getheight -1); $rowN++ )
{
my $col2 = $tileImage->getpixel( x=>$colN, y=>$rowN );
my ( $r,$g,$b ) = $col2->rgba();

# count all of the colors for averaging later
$rAvg += $r; $gAvg += $g; $bAvg += $b;

}#row
}#column

$cls{ $fName }{ r } = $rAvg / ($tileImage->getwidth * $tileImage->getheight);
$cls{ $fName }{ g } = $gAvg / ($tileImage->getwidth * $tileImage->getheight);
$cls{ $fName }{ b } = $bAvg / ($tileImage->getwidth * $tileImage->getheight);

}#for each color.tiles

After variable declaration and usage checks, each color-area file in the color.tiles directory is read. Using the Imager getpixel subroutine, the red, green, and blue color components for each pixel are recorded to be averaged after all the data points are collected. This step is critical to finding a closest match later on in the program.
The next section of code will read the Echo Tops image and record the appropriate setting of altitude based on the color of each pixel in the image.

Listing 2. altitudeMapper.pl read image-height data
########### Read height data ##############
my $altImage = Imager->new;
$altImage->read( file=> $fNameAlt ) or die "can't open $fNameAlt: ";
my $maxCol = $altImage->getwidth -1; # getwidth is a very slow read
my $maxRow = $altImage->getheight -1; # getheight is a very slow read

for( my $colN = 0; $colN <= $maxCol; $colN++ )
{
for( my $rowN = 0; $rowN <= $maxRow; $rowN++ )
{
my $currColor = $altImage->getpixel( x=>$colN, y=>$rowN );
my ($r,$g,$b ) = $currColor->rgba();

# skip if white or black pixel
next unless ( "$r,$g,$b" ne "0,0,0" && "$r,$g,$b" ne "255,255,255" );

my $totalDev = 10000;
my $closest = "none";

for my $key( keys %cls)
{
my $devR = abs( $cls{$key}{r} - $r );
my $devG = abs( $cls{$key}{g} - $g );
my $devB = abs( $cls{$key}{b} - $b );

# skip if deviation is greater than currently selected
next unless ( ($devR + $devG + $devB) < $totalDev );

$totalDev = $devR + $devG + $devB;
$closest = $key;

}#for key comparison

$hgt{ $closest }{ $colN }{ $rowN } = $currColor;

}#row

print "$colN\n" if( $colN % 100 == 0 ); # progress indicator

}#column

Images used in this example are 6,000 x 3,000 pixels, so reading from the getwidth and getheight subroutines is done before the for loops for a significant speed increase. Each pixel is then read, and its red, green, and blue deviations are measured from each of the average colors recorded in the previous step. Once the closest match has been found, that pixels' altitude setting is stored in the %hgt hash variable. Now the appropriate altitude information has been recorded, allowing the creation of image slices from a secondary image as shown below.

Listing 3. altitudeMapper.pl image slicer based on height
############ Apply height data to image ###########
my %layer = (); # pixel x,y data for that layer
my %avgLev = (); # average colors for that layer
my $lev = keys %cls; # number of image components to process
my $compImage = Imager->new;
$compImage->read( file=>$fNameSli ) or die "Cannot load $fNameSli: ";
$maxCol = $compImage->getwidth -1; # getwidth is a very slow read
$maxRow = $compImage->getheight -1; # getheight is a very slow read

# process the layers from top to bottom
for my $key( reverse sort keys %cls )
{
my $imgOut = Imager->new;
$imgOut->read( file=>"images/blank6000x3000.png") or die "blank image borked";

for( my $colN = 0; $colN <= $maxCol; $colN++ )
{
for( my $rowN = 0; $rowN <= $maxRow; $rowN++ )
{
# if layer data exists, and not on the top level
if( exists($layer{ $colN }{ $rowN }) && ($lev < (keys %cls)) )
{
# the color average from the previous layer
my $newColor = Imager::Color->new( $avgLev{ $lev+1 }{ r },
$avgLev{ $lev+1 }{ g },
$avgLev{ $lev+1 }{ b } );
$imgOut->setpixel( x=>$colN, y=>$rowN, color=>$newColor );

}#set to previous layer if exists

After declaring variables and opening the image to be "sliced," each layer in the %hgt hash is processed. The imgOut variable holds a transparent 6,000 x 3,000 canvas suitable for pixel placement and will be written to disk at the end of each pass. If layer data exists for the above layer (and there is a layer above), this average color data for the previous layer is written to the output image. This ensures that there are no "holes" in the image, where a lower layer would not contain the appropriate imagery from higher layers. As you can see in Figure 1, the bottom layers are a consistent color, even though higher layers contain different shades of yellow and green. Filling in these holes is one of the requirements of translating planar maps into visualizations more suitable for a 3-D environment.
Listing 4 overwrites these previous layer pixels if data for that layer exists in the $fNameSli (Merged Reflectivity Composite) image.

Listing 4. altitudeMapper.pl copy new slice pixels
# copy pixels from the original image if height data exists for this layer
if( exists( $hgt{ $key }{ $colN }{ $rowN } ) )
{
my $currColor = $compImage->getpixel( x=>$colN, y=>$rowN );
my($r, $g, $b) = $currColor->rgba();

$imgOut->setpixel( x=>$colN, y=>$rowN, color=>$currColor );
$layer{ $colN }{ $rowN } = $currColor;

$avgLev{ $lev }{ r } += $r;
$avgLev{ $lev }{ g } += $g;
$avgLev{ $lev }{ b } += $b;
$avgLev{ $lev }{ count } ++;

}# if height data exists for that pixel

}#row
}#column

In addition to copying the pixel from the $fNameSli image, the color data for that pixel is recorded for averaging. This ensures a uniform color background for each slice where it is filling in the holes from higher layers. Listing 5 writes out the final image and averages the color values for each layer.

Listing 5. altitudeMapper.pl write image out, compute averages for layer
print "write tiles/$key\n";
$imgOut->write( file=>"tiles/$key", type=>"png") or die "can't write $key";

# compute color averages for this level
if( exists( $avgLev{ $lev }{ count } ) )
{
$avgLev{ $lev }{ r } = $avgLev{ $lev }{ r } / $avgLev{ $lev }{ count };
$avgLev{ $lev }{ g } = $avgLev{ $lev }{ g } / $avgLev{ $lev }{ count };
$avgLev{ $lev }{ b } = $avgLev{ $lev }{ b } / $avgLev{ $lev }{ count };

}#if layer exists

$lev--;

}#for each key

Note that the average color of the $fNameSli (Merged Reflectivity Composite) pixels is recorded in the %avgLev hash for each image coordinate. Certain pixel colors may be present in the color-area tiles that are not present in the $fNameSli image, and this averaging helps smooth out the look of the image where the colors do not match precisely.


Back to top

altitudeMapper.pl usage
With the altitudeMapper.pl code in place, you need to create a directory for the sliced and output images, as well as create source transparent image. The first three lines in Listing 6 show the required commands.

Listing 6. Usage commands
mkdir images
mkdir tiles
convert -size 6000x3000 xc:none images/blank6000x3000.png
# or create with The Gimp

perl altitudeMapper.pl \
images/EchoTop_18_20080725-061737.png \
images/MergedReflectivityComposite_20080725-061750.png

The last three lines in Listing 6 show the usage of the altitudeMapper.pl program itself. Note that the \ characters are to indicate a line break for formatting purposes. The first image is the cloud-tops data, and the second image is the source to split into tiles.