Shaded relief in AS3

Shaded relief map of Molokai generated in AS3

Did you know that basic relief shading is fairly simple to accomplish? I didn’t until a few days ago when I was flipping through the Slocum et al. cartography textbook, Thematic Cartography and Geographic Visualization, wherein a four-step process is presented. It sounded easy, so it was time for further adventures in raster-based cartography in Flash! I swore off code projects outside my day job because they tend to be more headache- and hair-loss-inducing than enjoyable, but an easy one that produces pretty pictures is worth an exception.

I previously tried raster map projection in AS3, and the same caveats apply. Flash shouldn’t be your first choice for hill-shading, but it’s good to demonstrate that it can be done on the fly in AS3, which could occasionally come in handy. Here follows a tutorial-like explanation of how simple shaded relief can be achieved in AS3. Code is only shown in bits and pieces, so to see it all in one place (or if you want to skip the details), you can download the source AS3 file and the source image file.


1) THE DATA SOURCE

Any shaded relief map is born of elevation data, typically from a digital elevation model (DEM), a.k.a. digital terrain model. Ideally the relief map would be generated from the actual elevation values in the DEM, but a DEM is often visually displayed as a grayscale image, which is what my example uses as a source. (See the bottom of this post for a bit more info on DEM data.) Darker means lower elevation; lighter means higher. The price of using the image as a starting point is precision, as there are only 256 different gray values (elevations) in the image.

My source, an elevation map of the Hawaiian island of Moloka’i, is below.

Molokai DEM


2) LOADING AND PREPARING DATA

To generate the shaded relief it’s necessary to go through the DEM image pixel-by-pixel, calculate a result pixel for the relief, and draw that pixel to an an output image. So the setup for this whole thing is to load the image using a Loader, draw it to a BitmapData, and create another BitmapData for the relief map.

public class ShadedRelief extends Sprite
{
	private var loader : Loader = new Loader();
	private var sourceBitmap : BitmapData;
	private var reliefBitmap : BitmapData;
	public function ShadedRelief()
	{
		loader.load( new URLRequest("molokai.png") );
		loader.contentLoaderInfo.addEventListener(Event.COMPLETE,onLoaded);
	}
	private function onLoaded( event : Event ) : void
	{
		sourceBitmap = new BitmapData( loader.width, loader.height );
		sourceBitmap.draw( loader );
		reliefBitmap = new BitmapData( loader.width, loader.height,true, 0 );
	}
}


3) GET YOUR MATH ON

Now I start following the mathematical steps laid out in the Slocum et al. text, which are three. For each pixel in the grid:

  1. Calculate the slope of the land.
  2. Calculate the aspect of the land (the direction it faces).
  3. Calculate a reflectance value.

These calculations are made by observing the data in a 3×3 pixel window, at the center of which is the pixel of interest. So to start, loop through the rows and columns of the source image. There are, I think, more efficient ways to do pixel-level manipulation than what I spell out below, but I’m keeping it simple here.

private function drawMap() : void
{
	for ( var i : int = 0; i < sourceBitmap.width; i ++ ) {
		for ( var j : int = 0; j < sourceBitmap.height; j ++ ) {
			// magic will occur here
		}	
	}
}

For reference, here’s a diagram of the 3×3 window for any given pixel (i,j):

3x3 window of elevation values

The Z values represent elevation values for each pixel. To find the slope of the cell, do a simple rise over run division in both the x and y directions.

x slope equation
y slope equation

From the x and y slopes follows an overall slope:

overall slope equation

The grayscale image no longer has any real elevation values, so the 0-255 gray value will stand in. That can be obtained with the BitmapData.getPixel method and from that grabbing just one of the RGB channels, e.g. source.getPixel(i+1,j) & 0xFF. The D value (the distance between pixels) isn’t known in meaningful units anymore either, so D could just be 1. However, for some reason I found that using a slightly larger number produced a nicer result. My slope code looks like this:

var topValue : Number = sourceBitmap.getPixel( i, Math.max(j-1,0) ) & 0xff;
var leftValue : Number = sourceBitmap.getPixel( Math.max(i-1,0), j ) & 0xff;
var rightValue : Number = sourceBitmap.getPixel( Math.min(i+1,sourceBitmap.width-1), j ) & 0xff;
var bottomValue : Number = sourceBitmap.getPixel( i, Math.min(j+1,sourceBitmap.height-1) ) & 0xff;
 
var slx : Number = (rightValue - leftValue)/3;
var sly : Number = ( bottomValue - topValue )/3;
var sl0 : Number = Math.sqrt( slx*slx + sly*sly );

Next, the aspect. First, get a local angle between the x slope and the overall slope.

facet angle

To get an azimuth (0º to 360º, where north is 0º, east 90º, etc.), this table is provided:

Converting angle to azimuth

Code, then. Remember that AS3 works in radians, not degrees.

var phi : Number = Math.acos( slx/sl0 );
if ( sl0 == 0 ) { // account for division by zero trouble
	phi = 0;
}
var azimuth : Number;
if ( slx > 0 ) {
	if ( sly > 0 ) azimuth = phi + 1.5*Math.PI;
	else if ( sly < 0 ) azimuth = 1.5*Math.PI - phi;
	else phi = 1.5*Math.PI;
} else if ( slx < 0 ){
	if ( sly < 0 ) azimuth = phi + .5*Math.PI;
	else if ( sly > 0 ) azimuth = .5*Math.PI - phi;
	else azimuth = .5*Math.PI;
} else {
	if ( sly < 0 ) azimuth = Math.PI;
	else if ( sly > 0 ) azimuth = 0;
}

And now, a reflectance value, which will be the brightness of the final output pixel. A Lambertian reflectance is suggested, and the formula provided ends up looking like this in code:

var sunElev : Number = Math.PI*.25;
var sunAzimuth : Number = 1.75*Math.PI;
 
var L : Number = Math.cos( azimuth - sunAzimuth )*Math.cos( Math.PI*.5 - Math.atan(sl0) )*Math.cos( sunElev ) + Math.sin( Math.PI*.5 - Math.atan(sl0) )*Math.sin( sunElev );

The sunAzimuth and sunElevation values can be freely chosen and will affect the final appearance of shadows. 315º (northwest) and 45º respectively are pretty typical values. For some reason that I haven’t figured out, I ended up with some negative reflectance values, which I guess indicate black holes on the surface of the earth. I don’t understand the math going on here, but setting those to zero seemed to make the result look okay. (If anyone can explain or help me here, it’d be greatly appreciated!)


4) DRAW THE RELIEF

The reflectance value just needs to be translated into a shade of gray. The reflectance is basically a proportion of white, so to get RGB values it needs to be multiplied by 255 and assigned to each of the three channels. Finally, the pixel can be drawn on the output image.

var grayValue : Number = int(255 * L);
 
if ( grayValue < 0 ) {
	grayValue = 0;
}
// doing ARGB instead of RGB because transparency comes in handy elsewhere
reliefBitmap.setPixel32(i,j,0xff << 24 | grayValue << 16 | grayValue << 8 | grayValue);

Running through all that gets me the following Moloka’i map.

Molokai shaded relief grayscale


5) HYPSOMETRIC TINTING

The gray relief map wowed me enough the first time because I was amazed at how simple it was to produce (despite my lengthy explanation above). But it’s not a nice map without colors! Hypsometric tinting, in which colors are mapped to elevation ranges, is a common and aesthetically pleasing technique on relief maps. To get hypsometric tints on my relief map, another BitmapData is drawn based on the grayscale source image, but this time the math is much simpler. All that needs to be done is to match the gray value of each source pixel to a color class or a point in a continuous gradient of colors. I settled on the following gradient for elevation colors. It stays mostly in the greens to reflect the tropical setting, and moves a bit toward warmer colors at just the highest elevations.

Colors for hypsometric tints

I won’t explain all the code in detail, but basically I drew a vector gradient 256 pixels wide, then sampled the color at each x location and saved it in an array. Thus each 0-255 gray value can be matched to an index in that array. Then, in the loop where all that reflectance math is taking place, I added a couple of lines to draw that color to another BitmapData.

Note there is some additional funny business involving the first couple of values in the array. This is to deal with the ocean in my map. I had trouble getting the blues to look good in the gradient, so after creating the array I set those blues using brute force. I also didn’t want the ocean to show up as gray in the relief image, so if a pixel has a value at that low end that indicates it as water, it isn’t drawn to the relief map.

var gradientColors : Array = [0x0000ff,0x004000,0x679167, 0x81b279, 0xbfdfa8, 0xd0b8aa];
var gradientRatios : Array = [1,3,1*255/4,2*255/4,3*255/4,4*255/4];
var gradientAlphas : Array = [0,1,1,1,1,1];
var matrix : Matrix = new Matrix();
matrix.createGradientBox(255,5);
var gradientShape : Shape = new Shape();
gradientShape.graphics.beginGradientFill("linear",gradientColors,gradientAlphas,gradientRatios,matrix);
gradientShape.graphics.drawRect(0,0,255,5);
gradientShape.graphics.endFill();
 
gradBmp = new BitmapData(255,5);
gradBmp.draw(gradientShape);
 
var colors : Array = [];
 
for ( var n : int = 0; n < 256; n ++ ) {
	colors.push( gradBmp.getPixel32(n,2) );
}
// ARGB again
colors[0] = 0xffccccff;
colors[1] = 0xffccccff;
 
[...]
 
for ( var i : int = 0; i < sourceBitmap.width; i ++ ) {
	for ( var j : int = 0; j < sourceBitmap.height; j ++ ) {
		[...]
		var centerValue : Number = sourceBitmap.getPixel( i,j ) & 0xff;
		tintBitmap.setPixel32( i,j,colors[int(centerValue)] );
		[...]
		if ( centerValue >= 3 ) reliefBitmap.setPixel32(i,j,0xff << 24 | grayValue << 16 | grayValue << 8 | grayValue);
	}	
}

The hypsometric tinting by itself looks like this:

Molokai hypsometric tints


6) PUTTING IT ALL TOGETHER

The final map is a matter of combining the relief and tint images. I found that overlaying the relief map at 30% opacity looked pretty good. A blur filter on both the tint and relief maps helps to smooth out the pixelated coastline and some noise, respectively. And a glow filter on the relief map provides a nice coastal effect, even if I haven’t quite perfected it here.

tintBitmap.applyFilter( tintBitmap, new Rectangle(0,0,tintBitmap.width,tintBitmap.height), new Point(), new BlurFilter(1.2,1.2) );
var tintMap : Bitmap = new Bitmap(tintBitmap);
addChild(tintMap);
 
reliefBitmap.applyFilter( reliefBitmap, new Rectangle(0,0,reliefBitmap.width,reliefBitmap.height), new Point(), new BlurFilter(2,2) );
reliefBitmap.applyFilter( reliefBitmap, new Rectangle(0,0,reliefBitmap.width,reliefBitmap.height), new Point(), new GlowFilter(0xffffff,.75,45,45,3,2) );
var reliefMap : Bitmap = new Bitmap( reliefBitmap );
addChild( reliefMap );
reliefMap.alpha = .3;

Pow. And pau.

Shaded relief map of Molokai generated in AS3


7) ADDITIONAL RESOURCES

The formulas in the Slocum et al. textbook are derived from:
Eyton, J.R. (1991) “Rate-of-change-maps.” Cartography and Geographic Information Systems 18, no. 2:87–103.
As for the textbook itself, I’m looking at the second edition, pages 300-301.

It’s been a while since I actually sought out DEM data, but at least for the United States the best source remains the USGS. You can download data for anywhere in the country at http://seamless.usgs.gov/. The USGS also has some worldwide data with GTOPO30. A bit of Googling will likely lead you to additional sources of DEMs.

DEM data may not always come in a format suitable for this particular tutorial. (But remember that this is not the way to go most of the time anyway.) It may come in the actual DEM file format, or perhaps a GeoTIFF that you may have trouble opening in ordinary image editing software. Unfortunately I’m not up on the best ways to convert DEMs to grayscale images. I grew to like MicroDEM, a Windows program from several years ago, for all kinds of DEM manipulation. I’m a Mac user but have fired up Windows on my machine more than once to use this program. Anyway, again the best I can suggest is to use your Google-fu.

Finally, there is no avoiding mentioning the man who is barely short of a deity of relief shading: Tom Patterson. Mr. Patterson, who is known to and revered by all cartographers, does mind-blowing work for the National Park Service, and he has been kind enough to impart some of his wisdom through shadedrelief.com. Do check out that site for a lot of good information and tutorials. He also has created the wonderful Natural Earth I, II, and III maps and has co-spearheaded the Natural Earth vector effort.

Tagged , , ,

7 Comments

  1. Alternatively, to find the slopes / central difference, you can use a pixelbender filter or a combination of 2 convolution filters (one for x, one for y). This enables you to do this in realtime.

    Ralph Hauwert
    25 March 2010 @ 5:58am

  2. You might want to play around with blending modes for different effects when combining the maps.

    For example, try applying the relief map over the top with blend mode set to multiply:

    reliefMap.blendMode = BlendMode.MULTIPLY;

    :)

    Tim Oxley
    25 March 2010 @ 7:14am

  3. Ralph: thanks for that tip! I’m sort of aware of those things but don’t know much about them. I’ll be sure to check it out when possible.

    Tim: indeed, thank you for pointing that out. Working with this kind of thing in the past in Photoshop, I was always most pleased with just a regular ol’ semi-transparent overlay instead of different blending modes, but that’s definitely worth mentioning. Different color schemes will probably work best with different blending modes.

    Andy Woodruff
    25 March 2010 @ 10:05am

  4. In addition to Tom’s site I would add http://www.reliefshading.com which has a lot of detail on the theory and methodology behind effective shading. And it’s worth mentioning the book “Cartographic Relief Presentation” by Eduard Imhof, the man who made relief shading an academic subject.

    David Medeiros
    25 March 2010 @ 12:48pm

  5. Looks nice! You might find BitmapData.paletteMap() useful for hypsometric tinting. You can define a 256-element array to map the gray values into some other color, and Flash will do the loop and lookup for you.

    Amit Patel
    19 June 2010 @ 9:08pm

  6. Wow, how did I not know about that paletteMap method? Thanks!

    Andy Woodruff
    20 June 2010 @ 12:23pm

  7. That grayscale elevation map on section 1) is exactly what I have been needing. I’ve been trying to find ones of the other islands too but can and when I look up DEMs all I find is one that do the shadowing but I need it with out shadow and just grayscale elevation. Can someone help me find that?

    kaleo Puaa
    16 January 2011 @ 12:56pm