Skip to content

Terrain based on USGS data

ricardolpd edited this page Nov 3, 2012 · 10 revisions

Summary

For a game you can use artificial terrain or - as I did - use a real terrain to save the work of designing your heightmap. You can use that for other application than games though we'll consider some aspects of a walkable terrain.

Get the data

There is a tutorial that shows how to fetch elevation data from the USGS at http://www.gandraxa.com/digital_elevation_models.aspx. Look for an area that suits your needs and download the data in signed short (BIL) format. I took roughly one to one and a half degree which resulted in a 7.2MB compressed download. The BIL file is 37MB and the data is organized in 5400 columns and 3601 rows. Even outside of the US I got a one-second resolution (may be interpolated).

The download pack consists of several files but what we really need is just the BIL and the HDR file. The BIL consists the height data in signed short and the HDR tells you some details about the data you fetched.

This is the HDR file I got:

BYTEORDER      I
LAYOUT         BIL
NROWS          3601
NCOLS          5400
NBANDS         1
NBITS          16
BANDROWBYTES   10800
TOTALROWBYTES  10800
PIXELTYPE      SIGNEDINT
ULXMAP         -110.499861110207
ULYMAP         23.6668055539431
XDIM           0.000277777777799941
YDIM           0.000277777777799943
NODATA         32767

Considerations

If you keep the world scale of the DEM data you will end up with a resolution around 30m (one second) or more. If you will never put the camera near the terrain this may be sufficient but if you need a walkable terrain it should be around one meter.

For a fine grained heightmap of one meter a resolution of 30 meters can be considered the LOD level 5 data as a 32:1 reduction. Higher LOD levels can be computed as in the Ardor3D terrain examples which just skip values i.e. for one level higher we just take every second height point of the previous level. The lower levels have to be interpolated. I used a simple linear interpolation to start with.

A DEM data source can get quite big and it would be a waste of memory to keep it in RAM. I decided to use a mapped buffer which takes advantage of the underlying O/S memory mapping.

Implementation

I implemented a HeightmapPyramid as I don't intend to keep several heightmaps in memory. The logical layers (LOD levels) are served by this map pyramid.

The functions getHeight() and getSize() tell from which layer the data is requested in the first argument and we can use that for our calculations.

Loading a BIL file is very simple. It consists of signed short values for the height in one band i.e. all values from left to right (west to east) of the first row followed by the values of the following rows from top to down (north to south). The HDR file specifies what the byte order is.

The loader just supports one format to keep things simple. This is the constructor that reads the HDR, opens the file and initializes a small noise array:

public DEMHeightmapPyramid(File hdrfile, File mapfile) throws IOException {
		loadHDR(hdrfile);
		if (!("BIL".equals(layout)))
			throw new IOException("unsupported layout: " + layout);
		if (nbits != 16)
			throw new IOException("unsupported bit depth: " + nbits);
		raf = new RandomAccessFile(mapfile, "r");
		final FileChannel channel = raf.getChannel();
		mappedFile = channel.map(FileChannel.MapMode.READ_ONLY, 0, channel.size()).load();
		mappedFile.order(byteOrder);
		channel.close();
		interpolation = Interpolation.Linear;
		noiseFactor = 0.002f;
		final Random rnd = new Random(42);
		for (int i = 0; i < 256; i++)
			noise[i] = rnd.nextFloat();
	}

First we analyze the HDR file:

private void loadHDR(File f) throws IOException {
		BufferedReader r = null;
		try {
			r = new BufferedReader(new FileReader(f));
			while (r.ready()) {
				final String l = r.readLine();
				if (l.startsWith("NROWS")) {
					rows = Integer.parseInt(l.substring(6).trim());
				} else if (l.startsWith("NCOLS")) {
					columns = Integer.parseInt(l.substring(6).trim());
				} else if (l.startsWith("XDIM")) {
					xdim = Double.parseDouble(l.substring(6).trim());
				} else if (l.startsWith("YDIM")) {
					ydim = Double.parseDouble(l.substring(6).trim());
				} else if (l.startsWith("ULXMAP")) {
					xmap = Double.parseDouble(l.substring(6).trim());
				} else if (l.startsWith("ULYMAP")) {
					ymap = Double.parseDouble(l.substring(6).trim());
				} else if (l.startsWith("NBITS")) {
					nbits = Integer.parseInt(l.substring(7).trim());
				} else if (l.startsWith("BYTEORDER")) {
					byteOrder = "I".equals(l.substring(10).trim()) ? ByteOrder.LITTLE_ENDIAN : ByteOrder.BIG_ENDIAN;
				} else if (l.startsWith("LAYOUT")) {
					layout = l.substring(7).trim();
				}
			}
			// size = smallestPowerOf2Containing(Math.max(rows, columns));
			size = Math.max(rows, columns);
			// final double xmean = xmap + xdim * columns / 2.0;
			final double ymean = ymap + ydim * rows / 2.0;
			// 1 degree = 111.32 km, ydim is resolution in degree
			yres = 111320.0 * ydim;
			// one minute = 1852 m * cos(a)
			xres = 1852.0 / (xdim * 60 * columns * Math.cos(ymean));
			final int lvlx = nearestLevel(xres);
			final int lvly = nearestLevel(yres);
			if (lvlx != lvly)
				throw new IOException("x and y resolution are too different to build one LOD level");
			suggestedLevel = lvlx;
		} catch (final IOException ex) {
			throw ex;
		} finally {
			try {
				r.close();
			} catch (final Exception ex) {
			}
		}
	}

This function extracts some parameters and calculates the level from the resolution. In my case the x and y resolution was near 30 so the suggested LOD level fo the data is 5 (2^5 = 32). The function to get the level is:

private int nearestLevel(double res) {
		int level = 0;
		double previousFault = Double.MAX_VALUE;
		for (;;) {
			final double fault = Math.abs((1 << level) - res) / res;
			if (fault > previousFault)
				break;
			previousFault = fault;
			level++;
		}
		logger.log(Level.INFO, "level {0,number} for {1,number,#.###}, relative fault {2,number,#.##}%", new Object[] {
				(level - 1), res, previousFault * 100 });
		return level - 1;
	}

Now let's take care of two methods of the HeightmapPyramid interface.

@Override
	public int getHeightmapCount() {
		return 20; // suggestedLevel + 4;
	}

	@Override
	public int getSize(int level) {
		if (level >= suggestedLevel)
			return size / (1 << (level - suggestedLevel));
		else
			return size * (1 << (suggestedLevel - level));
	}

It doesn't hurt to keep the heightmap count high as the calculation is cheap. It could be dependent on the native level or just a fixed value. The size of the heightmap is calculated using the requested level.

Next we need a function to get a height value from the mapped memory. The byte order is already set so that is easy as well. We just have to keep in mind that the value from the BIL is in meters above the sea level while we want to get float values between 0 and 1 or if we really want that we could allow -1 to 1.

private float getDirectHeight(int x, int y) {
		final int index = (y * columns + x) * 2;
		if (index < 0 || index >= mappedFile.limit())
			// throw new IndexOutOfBoundsException("index out of range: " + x +
			// "/" + y);
			return defaultHeight;
		// scale to -1..1
		return mappedFile.getShort(index) / heightScale;
	}

The above function can be used for the native level directly but we need some more logic to serve all levels:

@Override
	public float getHeight(int level, int x, int y) {
		// short cut: if the level is our "native" level
		if (level == suggestedLevel)
			return getDirectHeight(x, y);
		// higher level: just skip values
		if (level > suggestedLevel) {
			final int scaled = 1 << (level - suggestedLevel);
			return getDirectHeight(x * scaled, y * scaled);
		}
		// lower level: need to interpolate
		final int scaleToNative = 1 << (suggestedLevel - level);
		return (float) getInterpolatedHeight((float) x / scaleToNative, (float) y / scaleToNative);
	}

As mentioned before we just skip values for higher levels and interpolate for lower. The interpolation can be done in different ways and here are three to play with:

public static enum Interpolation {
		None, Linear, Noise
	};

	private Interpolation interpolation;

The Interpolation.None will give stairs so you won't need it. The Noise is very simple and could probably be replaced with a more sophisticated function that uses turbulence. The Linear just does what it says: it interpolates linearly.

protected double getInterpolatedHeight(final float x, final float z) {
		final int col1 = (int) Math.floor(x);
		final int row1 = (int) Math.floor(z);

		final float intOnX = x - col1, intOnZ = z - row1;

		float height = defaultHeight;
		if (interpolation == Interpolation.None)
			return getDirectHeight(col1, row1);
		final int col2 = col1 + 1;
		final int row2 = row1 + 1;
		final float topLeft = getDirectHeight(col1, row1);
		final float topRight = getDirectHeight(col2, row1);
		final float bottomLeft = getDirectHeight(col1, row2);
		final float bottomRight = getDirectHeight(col2, row2);
		height = MathUtils.lerp(intOnZ, MathUtils.lerp(intOnX, topLeft, topRight), MathUtils
				.lerp(intOnX, bottomLeft, bottomRight));
		if (interpolation == Interpolation.Noise)
			height += (noise[Float.floatToIntBits(x + z) & 0xff] - 0.5f) * noiseFactor;
		return height;
	}

That's all for the heightmap pyramid. You can get the full source from SVN at: http://code.google.com/p/worldofmystery/source/browse/trunk/worldofmystery/src/de/wom/a3dclient/terrain/DEMHeightmapPyramid.java

Usage of the heightmap

Now let's build a clipmap terrain that uses the heightmap pyramid. For the time being I use a copy of the procedural texture streamer. That uses a function so let's just wrap the getter in a Function3D:

public class DEMTerrain {

	private final TexturedGeometryClipmapTerrain geometryClipmapTerrain;
	private final DEMHeightmapPyramid heightmapPyramid;
	private final int clipSideSize = 127;
	private final int textureSliceSize = 256;

	public DEMTerrain(Camera camera, File hdrfile, File mapfile) {
		try {
			heightmapPyramid = new DEMHeightmapPyramid(hdrfile, mapfile);
		} catch (final IOException e) {
			throw new RuntimeException(e);
		}

		final DesignTextureStreamer streamer = new DesignTextureStreamer(65536, textureSliceSize, heightmapPyramid,
				heightmapPyramid.getHeightScale());

		final Function3D func = new Function3D() {

			@Override
			public double eval(double x, double y, double z) {
				return heightmapPyramid.getHeight(0, (int) x, (int) y);
			}
		};
		final HeightBasedTextureStreamer streamer2 = new HeightBasedTextureStreamer(65536, textureSliceSize, func,
				heightmapPyramid.getHeightScale());

		final int validLevels = streamer2.getValidLevels();
		final TextureClipmap textureClipmap = new TextureClipmap(streamer2, textureSliceSize, validLevels, 64f);

The above mentioned texture streamer looks very familiar if you check the Ardor3D examples. You can find it at http://code.google.com/p/worldofmystery/source/browse/trunk/worldofmystery/src/de/wom/a3dclient/terrain/HeightBasedTextureStreamer.java

Note: This should be replaced by some better implementation.

Finally we can build a TexturedGeometryClipmapTerrain:

	geometryClipmapTerrain = new TexturedGeometryClipmapTerrain(textureClipmap, camera, heightmapPyramid, clipSideSize,
			heightmapPyramid.getHeightScale());
	geometryClipmapTerrain.setHeightRange(-1f, 1f);

The complete class that builds the terrain is at http://code.google.com/p/worldofmystery/source/browse/trunk/worldofmystery/src/de/wom/a3dclient/terrain/DEMTerrain.java.