 | Level: Intermediate Nathan Harrington, Programmer, IBM
23 Sep 2008 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.
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.
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.
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.
Generating KML
You can browse the output image files located in the tiles directory or continue with
Listing 7.
Listing 7. kmlGenerator.pl program
#!/usr/bin/perl -w
# kmlGenerator.pl - create kml for sliced images created by altitudeMapper.pl
use strict;
die "specify starting altitude, layer distance" unless @ARGV==2;
my ( $alt, $distance ) = @ARGV;
my $drawOrder = 1;
print qq{<?xml version="1.0" encoding="UTF-8"?> \n};
print qq{<kml xmlns="http://earth.google.com/kml/2.0"> \n};
print qq{<Folder> \n};
print qq{ <name>kmlGenerator output</name>\n};
while( my $line = <STDIN> )
{
chomp($line);
print " <GroundOverlay>\n";
print " <drawOrder>$drawOrder</drawOrder>\n";
print " <Icon><href>$line</href></Icon>\n";
print " <altitude>$alt</altitude>\n";
print " <altitudeMode>absolute</altitudeMode>\n";
print " <name>$drawOrder</name>\n";
print " <LatLonBox>\n";
print " <north>51.000</north>\n";
print " <south>21.000</south>\n";
print " <east>-67.000</east>\n";
print " <west>-127.000</west>\n";
print " </LatLonBox>\n";
print " </GroundOverlay>\n";
$drawOrder++;
$alt += $distance;
}#while line in
print "</Folder>\n";
print "</kml>\n";
|
The kmlGenerator.pl program wraps the appropriate KML around the sliced images.
Various combinations of starting altitude and distance between layers will produce
different volumetric effects at specific distances and viewing angles. Listing 8
demonstrates the program usage with the parameters used to generate the image shown on
the bottom of Figure 1.
Listing 8. kmlGenerator.pl usage
ls -1 tiles/* | perl kmlGenerator.pl 20000 10000 > kmlOutput.kml
|
Open the kmlOutput.kml file in Google Earth to view your new altitude attribute-enhanced imagery.
Conclusion, further examples
As evidenced, slicing the images into specific components based on imagery data in a
separate image is an effective method for adding a third dimension of data to your
formerly flat visualizations. These color-averaging and layer void-filling techniques
are only one way to achieve similar effects. Consider drawing polygons tracing the
edge of certain layers with heights equal to the cloud tops. Or add transparent
graphics with computed transparencies perpendicular to the Earth plane to create a more
realistic layered effect for all viewing angles.
Download | Description | Name | Size | Download method |
|---|
| Sample code | os-google-earth-altitude-echoTops.zip | 4KB | HTTP |
|---|
Resources Learn
Get products and technologies
- Innovate your next open source development project with IBM trial software, available for download or on DVD.
- Download IBM product evaluation versions, and get your hands on application development tools and middleware products from DB2®, Lotus®, Rational®, Tivoli®, and WebSphere®.
Discuss
About the author | 
| | Nathan Harrington is a programmer at IBM currently working with Linux and resource-locating technologies. |
Rate this page
| |
