Degenerate Shape Files

The shapefile specification is more like a pipe dream than a file format; it is a document that describes what people wish would happen in shapefiles. Examples of things one might wish for:

• No self-intersecting, non-simple, or generally hosed polygons.
• Polygon windings are CW for outer boundaries and CCW for holes.
• Whatever convention a shapefile decides to follow, the same convention is used for the entire file.

Unfortunately shapefiles don't always work that way. I discovered that some of the SWBD tiles violate the winding rule.

Fortunately times are changing: with so much open source GIS going on, often you can find a log of other people's misery and pretty quickly get a handle on the problem.

Of particular interest is the fight between correctness and performance in the trac log. How do we rapidly import our data if we don't know how good it is? Checking data quality is expensive?

Ways To Interpret a Broken Shapefile

There are a few ways to interpret a broken shapefile:

• We can go strictly by winding order, even when it is wrong, and perform a union on overlapping areas.
• We can use a "binary" approach to polygons (think of what you get when you rasterize a non-simple polygon) where we switch inside and outside every time we cross an edge.
• We can simply throw out data that violates a spec requirement.
• We can assume (for better or for worse) that the first winding is outer, the rest are inner.
• We can do a mix of this based on heuristics (e.g. what OGR does: if the windings look sane, trust them, otherwise slowly calculate ring nesting).

I have chosen the "binary" approach for the following reasons:

• It is never dependent on the order of windings in the file or the direction of windings.
• It works correctly with polygons with holes, polygons with islands in holes, and multipolygons (multiple disjoint outer rings), as long as the polygons are non-overlapping.
• When polygons do overlap, the results are reasonably acceptable.
• The implementation is really fast.

To understand that last point, we have to look at how CGAL processes polygons.

A Sweeping Improvement in Performance

CGAL implements sweep-line algorithms. It's beyond this post to explain how it works, but the main thing to note is that sweep-line processes a big pile of overlapping line segments in O(nlogn+klogn) time where n is the number of curves and k is the number of intersections. When adding/merging or inserting a large number of line segments into a planar map, this usually beats the pants off incremental modification algorithsm.

(A particularly painful case is building a very large ring in an empty map...each insert has to check for intersections with every other segment in the map with basically no spatial indexing or speed-up. Ouch!)

Of course we can go a lot faster when we know our data doesn't overlap, but go look at the title of the blog post...we're dealing with shape files here! We'll count ourselves lucky if polygons have 3 or more points.

Sweep-line is found in a number of interesting places:

• Bulk-adding line segments to an empty (or very full) map.
• Overlaying one map with another.
• Boolean operations between two sets of polygons.
• Checking whether a polygon is simple. (This is done with a stripped down sweep algorithm, but time complexity should be about the same.)

This list is interesting because approaches that use these different APIs are going to run at about the same time complexity. The number of times we use these operations might matter a lot.

Adding Up Lots Of Polygons

CGAL's general polygon set provides bulk insertion routines - given a range of polygons (or polygons with holes) you can join/difference them all at once. Looking at the code, it appears to use a divide and conquer scheme. I think this should be more optimal than a series of individual merges (which has to be slow, because by the end we are re-merging the finished data over and over), but frankly I can't tell where the various edges get consolidated (if ever).

Here's the problem: if we want to insert polygons with holes, we need to insert non-degenerate ones. (If we insert polygons, then subtract the holes, the results are not the same!) But...non-degenerate polygons, where do we get those?

Rather than an incremental approach (clean each winding, organize the windings, build sane polygons with holes, then merge them) I observe that we don't need clean input to use the "binary" rule for polygon interiors.

So my new shapefile processor simply inserts every single line segment from the entire file into an empty arrangement at once. This is about as fast as we can get - a single sweep that simultaneously deals with overlapping polygons, degenerate windings, and incorrect winding relationships.

I then do a breadth-first-search from the outside of the map to the inside, assigning properties. My arrangement tracks the original providing polygon so each time I "cross" a boundary I can track the change in ownership via a set.

The memory cost is O(F + E) where F is the number of faces in the final imported arrangement and E is the number of polygon edges. The time complexity is O(ElogE+KlogE + F) where K is the number of intersections our shape file turned out to have. (In practice, that F doesn't matter, we're almost guaranateed to have more edge than situations in a real world situation.

GeoTIFF and Off-By-One

I fear I have stumbled into GeoTIFF off-by-one hell. There is surprisingly little written about this, considering that (1) GeoTIFF is somewhat widely used and (2) the whole point of GeoTIFF is that the file tells you where the data goes. So...here is a summary of what I've stumbled into.

GeoTIFF has a flag to indicate whether pixels are "pixel is area" or "pixel is point". To summarize:

• Pixel is area means that each pixel is a rectangle, with infinitely thin coordinates bounding all four sides. The sample data represents a summary of the entire area, or maybe a sample at the center of the area.
• Pixel is point means that each pixel is the infinitely tiny intersection of two infinitely thin grid-points and represents the data as sampled at that point.

Who cares? Well, it turns out it matters. Imagine we have "3 second" data, which means that there are 1200 samples per degree, each one approximately 90 meters apart. (This is geographic projection, so the postings vary with latitude - I mention 90 meters only for the "feel" of the data.)

Let's consider what we have to do if we want data that covers exactly one degree.

If we are using area pixels, that one degree needs 1200 samples, each one 1/1200th of a degree wide. The left edge of the left most sample would be on a longitude line and the right edge of the right-most sample would be too. The next tile should have no data duplication from this tile.

If we are using point pixels, that one degree needs 1201 samples, each 1/1200th of a degree apart, but with no width on their own. The longitude line of the left edge has samples exactly on it, as does the right edge. And in fact, the right edge of our tile is a duplicate of the left-edge of the previous tile.

We could choose to exclude two of the four edges to in the pixel-center case to avoid duplication, but when we're tiling with pixel-center data, it is sort of nice to have the entire tile.

SRTM

The SRTM is pixel-center, at least in the raw data that JPL provides; we can only assume that every derived product should be that way too.

When you download a seamless SRTM tile from the USGS you get 1201 pixels, tagged as area pixels, and a bounding box that goes 1.5 arc-seconds outside the bounds you might expect. This is a little bit strange (IMO) but more or less correct.

The CGIAR SRTMs in GeoTIFF appear to be shifted by half a pixel: while the ARC ASCII files come in 6001 blocks that are clearly point-pixel (and this is what we would expect from a SRTM-derived product) the GeoTIFF version provides 6000 postings, but is listed as area-pixel with bounds on tile boundaries. Visual inspection indicates that the south and east rows have been dropped to avoid duplication; this would imply that the bounding box should be inset on the south and east edges by 1.5 arc seconds (half a pixel) and outside on the other edges.

Follow-Up 4/1/11: the v4 CGIAR SRTM files appear to have been recut at 6001x6001 for correct pixel-center formatting. The GeoTIFFs are marked as area pixels and the bounding box is expanded by 0.5 samples on all sides.

NED

NED appears to truly be area aligned (unlike SRTM); every form I have seen has been 3600 pixels centered on a lat-lon tile. In this case area pixels is an appropriate way to tag the data set.

To Sweep or Not To Sweep

I've spent the week porting the X-Plane scenery generation code back to CGAL.  I will describe the history of this code and its constant battle with robustness in another post, but for now I will just mention that credit goes (besides obviously to the authors of CGAL, who are geniuses of the Nth degree) to Andrew McGregor for a lot of the ideas on how to use CGAL with the scenery generation code - the final design is really just a refinement of his work.

The X-Plane scenery generation code (which we just call "the render farm") is the code that turns hundreds of GB of GIS data into scenery.  At its heart is a topologically integrated planar map of all vector features in a given tile - every road, river, coastline, change of terrain type, etc. is baked into one map.  A typical map might have half a million edges, covering a 1 x 1 degree area.

For X-Plane 8.20 and 9.0 we shipped using a planar map I wrote myself.  While the planar map didn't work very well, it had a few special cases that were highly optimized.  In particular, given two sub-regions of two maps, with an identical ring of half-edges identifying them (identical meaning identical coordinates) you could swap the interiors of the two regions in time proportional to only the number of elements inside the ring.

It turns out that this is a really useful trick for a few things:

• If you have a huge map you want to tile, you simply induce the cut lines, and swap each single tile with an empty map.
• If you want to "overlay" a polygon, or even a sub-map (destroying what is below) you simply induce a border ring in both polygons and swap.
• You can crop a map's interior or exterior by inducing a ring and swapping with a blank map.

This was an operation where the devil is in the details - because components are swapped "in-place" the swap doesn't have to copy, allocate, or deallocate the map, so it turns out to be fast in pracatice.  In terms of time complexity, it's very fast for very small operations on a huge map, which turns out to be an important case for us.

Now CGAL's arrangement_2 gives me some power tools* that are very useful - in particular, sweep-line algorithms.  The link shows a good explanation.  Basically the power of the link is to process a set of line segments in a lot less than O(N^2).

But there's a slight problem with sweep-line - it has to iterate over the entire set of edges in a map, including the ones already in the map.  That means that inserting a small number of edges into a complex map is surprisingly slow.

So in some cases I now have two implementations of various map operations:

• Sweep-based implementations, which are good when we're inserting a lot of data at once.
• Incremental implementations - theoretically much worse, but acceptable to keep the size of our data set down.

So far, incremental operations are proving to be a big win when "burning airports" - that is, merging or overlaying very small airport boundary areas into a full map.

* Of course, the real advantage of CGAL is that it is robust and doesn't randomly blow up during operations!