Quickly find and render terrain above a given elevation

Question

Given an elevation map consisting of lat/lon/elevation pairs, what is the fastest way to find all points above a given elevation level (or better yet, just the the 2D concave hull)?

I'm working on a GIS app where I need to render an overlay on top of a map to visually indicate regions that are of higher elevation; it's determining this polygon/region that has me stumped (for now). I have a simple array of lat/lon/elevation pairs (more specifically, the GTOPO30 DEM files), but I'm free to transform that into any data structure that you would suggest.

We've been pointed toward Triangulated Irregular Networks (TINs), but I'm not sure how to efficiently query that data once we've generated the TIN. I wouldn't be surprised if our problem could be solved similarly to how one would generate a contour map, but I don't have any experience with it. Any suggestions would be awesome.

Sachet

Answer 1

It sounds like you're attempting to create a polygonal representation of the boundary of the high land.

If you're working with raster data (sampled on a rectangular grid), try this.

Think of your grid as an assembly of right triangles.

Let's say you have a 3x3 grid of points

• a b c
• d e f
• g h k

Your triangles are:

• abd part of the rectangle abed
• bde the other part of the rectangle abed
• bef part of the rectangle bcfe
• cef the other part of the rectangle bcfe
• dge ... and so on

Your algorithm has these steps.

1. Build a list of triangles that are above the elevation threshold.

2. Take the union of these triangles to make a polygonal area.

3. Determine the boundary of the polygon.

4. If necessary, smooth the polygon boundary to make your layer look ok when displayed.

If you're trying to generate good looking contour lines, step 4 is very hard to to right.

Step 1 is the key to this problem.

For each triangle, if all three vertices are above the threshold, include the whole triangle in your list. If all are below, forget about the triangle. If some vertices are above and others below, split your triangle into three by adding new vertices that lie precisely on the elevation line (by interpolating elevation). Include the one or two of those new triangles in your highland list.

For the rest of the steps you'll need a decent 2d geometry processing library.

If your points are not on a regular grid, start by using the Delaunay algorithm (which you can look up) to organize your pointss in into triangles. Then follow the same algorith I mentioned above. Warning. This is going to look kind of sketchy if you don't have many points.

Bost

DON'T use this algorithm for civil engineering purposes! Unless you want massive puddles in the highway you're designing! Get a trained cartographer and high-resolution surveys for that purpose. – Bost

Answer 2

Assuming you have the lat/lon/elevation data stored in an array (or three separate arrays) you should be able to use array querying techniques to select all of the points where the elevation is above a certain threshold. For example, in python with numpy you can do:

indices = where(array > value)

And the indices variable will contain the indices of all elements of array greater than the threshold value. Similar commands are available in various other languages (for example IDL has the WHERE() command, and similar things can be done in Matlab).

Once you've got this list of indices you could create a new binary array where each place where the threshold was satisfied is set to 1:

binary_array[indices] = 1

(Assuming you've created a blank array of the same size as your original lat/long/elevation and called it binary_array.

If you're working with raster data (which I would recommend for this type of work), you may find that you can simply overlay this array on a map and get a nice set of regions appearing. However, if you need to convert the areas above the elevation threshold to vector polygons then you could use one of many inbuilt GIS methods to convert raster->vector.

Brythonic

While that would definitely work, but I'm afraid it would be too inefficient - I'm working with something on the scale of 1 elevation point per square 1/2 mile, for the entire U.S. (although I will only care about the 100 or so square miles currently on the screen). I'll need to update the screen very rapidly, as we're wanting to visually warn pilots of upcoming obstacles. – Sachet

Sorry for dumb Q: do you know the aircraft travel vector? Would it be a good place to start for optimising the data processing by filtering against the lat/long for the space the craft is moving towards? – Vento

Answer 3

I would use a nested C-squares arrangement, with each square having a pre-calculated maximum ground height. This would allow me to scan at a high level, discarding any squares where the max height is not above the search height, and drilling further into those squares where parts of the ground were above the search height.

If you're working to various set levels of search height, you could precalculate the convex hull for the various predefined levels for the smallest squares that you decide to use (or all the squares, for that matter.)

Scintillant

Consider scoring the lat/long/alt triplets such that the closest ones score the highest. Then handle them in batches of relative distance (10 miles, 20 miles, 50 miles ...? How fast is the craft going to be travelling ...? This is not an answer, but perhaps a pointer in a new direction? – Vento

Answer 4

I'm not sure whether your lat/lon/alt points are on a regular grid or not, but if not, perhaps they could be interpolated to represent even 100' ft altitude increments, and uniform lat/lon divisions (bearing in mind that that does not give uniform distance divisions). But if that would work, why not precompute a three dimensional array, where the indices represent altitude, latitude, and longitude respectively. Then when the aircraft needs data about points at or above an altitude, for a specific piece of terrain, the code only needs to read out a small part of the data in this array, which is indexed to make contiguous "voxels" contiguous in the indexing scheme.

Of course, the increments in longitude would not have to be uniform: if uniform distances are required, the same scheme would work, but the indexes for longitude would point to a nonuniformly spaced set of longitudes.

I don't think there would be any faster way of searching this data.

Carhart

Answer 5

It's not clear from your question if the set of points is static and you need to find what points are above a given elevation many times, or if you only need to do the query once.

The easiest solution is to just store the points in an array, sorted by elevation. Finding all points in a certain elevation range is just binary search, and you only need to sort once.

If you only need to do the query once, just do a linear search through the array in the order you got it. Building a fancier data structure from the array is going to be O(n) anyway, so you won't get better results by complicating things.

If you have some other requirements, like say you need to efficiently list all points inside some rectangle the user is viewing, or that points can be added or deleted at runtime, then a different data structure might be better. Presumably some sort of tree or grid.

If all you care about is rendering, you can perform this very efficiently using graphics hardware, and there is no need to use a fancy data structure at all, you can just send triangles to the GPU and have it kill fragments above or below a certain elevation.

Galangal