GeoSpiel

Once apon a time, a long long time ago, when MySQL first came out with spatial support, I looked at it with Dave Blasby and wrote an e-mail, which become the first-and-last word in terms of comparisons. (Update: Actually the last and best word is this real-world experience documented by RedFin.)

Since I was learning a bit about JMeter and JDBC testing last week, returning to the question of MySQL and PostGIS performance seemed like a reasonable thing to do.

First, many overall impressions remain the same as five years ago: the amount of functionality in MySQL Spatial remains very very small. You can do simple store-and-retrieve operations. Many of the spatial operations that are standard in full spatial databases don’t exist or are (confusingly, as some users have commented) stubbed out against bounding box tests instead. The documentation notes:

Currently, MySQL does not implement these functions according to the specification. Those that are implemented return the same result as the corresponding MBR-based functions. This includes functions in the following list other than Distance() and Related().

I encountered none of the stability or consistency issues that we experienced five years ago, MySQL worked as documented throughout.  So, what is left to test? The only thing MySQL users care about: raw query performance.

For testing I used an extract of TIGER data, 3.5M roads in the state of Texas. I loaded this data from a Shape file into MySQL and PostGIS using the ogr2ogr data conversion tool.

Once the data was loaded, I built a spatial index on the geometry columns.  For MySQL, this took about 4 minutes, and for PostGIS about 2 minutes.

Then I got to the meat of the testing, a concurrent load test. The test measures spatial query speed — given a bounding box, how long does it take to return the result set? In order to ensure that both databases got the same bounding boxes, I built a single file of bounding boxes that varied in size and location and all returned at least 20 features, but returned no more than 100. I ran this file through JMeter CSV Dataset Converters to create the load.

In addition applying a raw load, I also measured while varying two other parameters: the number of simultaneous threads working against the database; and, the percentage of queries during the load that were insert queries. The more insert queries, the more the load is a read/write load.

The results with 0% insert queries (no write load at all) look like this:

Both databases ramp up their throughput as the threads increase. PostGIS does somewhat better, but not orders of magnitude.

And what about when we add in a write load? How does query performance change? Using the load at 20 threads as the baseline, here’s what the numbers look as we vary the write load.

Surprisingly neither MySQL or PostGIS blink as the write load moves up to a fairly high proportion of the queries. I was expecting MySQL to fall down as writes increased, but it did not. The reason (I think) is that the reports of bad MySQL concurrency I had heard were around the InnoDB table type, and the spatial support is in MyISAM tables. So it’s possible that some inserts were being lost or ignored during selects (because MyISAM is not ACID) but the throughput was unaffected as writes went up.

I also tried some more complex spatial SQL, like spatial table joins, but the results did not reflect well on MySQL.  The MySQL query planner does not seem to have a good grasp of how to optimize a spatial query.

Recently a few of us developers were reviewing an OpenLayers map mixing raster and vector layers and we found that tiles for the vector layer were around 20kb on average, while those in the raster layer were more like 200kb.  Here’s a similar example based on the Blue marble data from sigma.openplans.org:

A tile from the Blue Marble layer in PNG format, about 132kB in size

The same tile, in JPEG format, about 11kB in size.

There’s a factor of about 10 in difference between the file size of these tiles, and that equates to a factor of 10 difference in speed for tile loading.  As it turned out, the map was just hard-coded to use PNG for all layers when JPEG would have made more sense in this case.

To explain, we need to examine the workings of the different image formats that GeoServer supports.  To simplify things a bit, I will focus just on PNG and JPEG, as these are reasonable representatives of two common families of image file formats. The PNG file format uses a lossless compression algorithm, while JPEG uses a lossy compression algorithm designed to reduce the size of images by removing information that isn’t noticed by the human eye anyway. But first, I’ll start off with a brief discussion of BMP, as it is one of the simplest image file formats there is.

In computer graphics, colors are generally represented as tuples of numbers, such as values for the Red, Green, and Blue channel in a color. (The exact meaning of these colors is called a color model, and, while RGB(A) is used on the web, there are plenty of others to check out if you’re interested).  Raster graphics—images that are defined in terms of pixels—are simply 2-dimensional matrices of these tuples.  The “resolution” or range of possible values for each channel determines how well a color can match the intended value.  For example, an HTML-style color specification like #8CC63F (aka GeoServer Green) includes 2 hexadecimal digits for each channel, for a total of 3 bytes or 24-bits per pixel.  In a BMP (BitMaP) image, these pixel values are stored, uncompressed, directly in the file. So, as an example, a 256-pixel-square BMP image using 24-bit colors will take 256 * 256 * 3 = 196608 bytes, plus some extra header information, totaling around 200kb for each tile.

By comparison, the PNG (Portable Network Graphics) image format uses a form of lossless compression based on differences in pixel values.  Because of the way this is implemented, solid fields of a single color compress very well, but file sizes increase dramatically when color varies a lot within the image.  In 24-bit mode, this is the only compression applied, but an 8-bit mode PNG also uses a color index.  A color index limits the space required to specify a color by predetermining a list of colors used and assigning each one a number.  With an 8-bit image, this restricts the image to only using 256 distinct colors.  The key points to take away are that PNG images tend to be very efficient when the image contains few colors and large solid fields, both typical of rendered vector layers.  You can see the imperfections introduced more easily by creating a difference image, subtracting all the individual color values in one picture from the corresponding values in another.  The darker the resulting image, the better; black (0 for each channel) implies that the color was reproduced exactly in the compressed image.  I created this image by loading both images into the GIMP and setting the layer mode to ‘Difference.’  You can see that there is very little difference between the two, mostly at boundaries between states where anti-aliasing has introduced extra color variation to this mostly flat image.  The differences here are pretty small, so you will need to click on the image to see it full size to see them.

The difference (right) between 24-bit PNG (left) and 8-bit indexed PNG (center)

On the other hand, JPEG (named after the Joint Photographic Experts Group that designed it) uses a more complicated compression scheme taking advantage of details of human color perception. For example, variations in brightness are more prominent than variations in hue.  Because JPEG encoding “throws out” information that wouldn’t be noticed by humans anyway, it is a lossy compression format, meaning that just repeatedly loading and re-encoding the same image in JPEG format can degrade its quality.  Since JPEG takes advantage of variation in color, it can efficiently store images that would result in large PNG images.  However, storing line art in JPEG usually introduces color variations, causing noticeable artifacts.  Here’s a difference image for PNG vs. JPEG for a raster layer.  Even scaled down, you can see the differences, but they are not really noticable in the JPEG image itself.

The difference (right) between the 24-bit PNG (left) and JPEG (right) versions of a Blue Marble tile.

So to sum up: vector layers should generally use PNG, and raster layers should generally use JPEG. If you are overlaying vector and raster data in the same WMS request, you can try GeoServer’s OpenLayers preview to easily switch between formats for comparison of visual artifacts versus network traffic.  Make sure to give PNG8 a try; it is a pseudoformat exposed by GeoServer that switches it to 8-bit indexed mode instead of the default 24-bit PNGs.  Often, this produces much smaller tiles with little apparent change in the rendering.

A figure-eight is a good pattern for ice-skating, but a bad pattern for a polygon. It is, in PostGIS terminology, “invalid”, it is a ring that self-intersects. It is very easy to draw figure-eights and other forms of invalid geometry in most GIS desktops, and so any data loaded into a spatial database comes with an in-built uncertainty — is the data valid or not?

Figure Eight

Validity is important because geometry algorithms depend on some basic structure in order to return correct results. Run a standard area calculation on a figure-eight polygon and you will not get the area that, to your human eye, appears to be bounded by the two lobes of the figure-eight.

Here is the text representation of a figure-eight shape, each lobe being a square 100 units on a side:

POLYGON ((100 100, 100 200, 300 200, 300 300, 200 300, 200 100, 100 100))

Perceived area: two lobes at 10,000 square units a piece = 20,000 units.

Calculated area: 0 units.

Why? Because the lobes of the figure-eight run in opposite directions, the calculated area of one cancels out the calculated area of the other. The data is invalid, or rather, what we visually perceive the data to be (two bounded areas touching at a point) is not the same as the way the computer interprets it. The algorithm used to fill in the pixels on a screen (scan-line fill) is not as sensitive to ring semantics as the algorithm for calculating area.

So, what if you are building a production database, and you don’t know what kind of data people will be entering, and you don’t want to check your data loads by hand. What is a data manager to do?

Create a RULE!

The PostgreSQL “rules system” is perfect for solving this problem.

Create two tables, one for your working data and one for data determined to be invalid.

CREATE TABLE mytable (id INTEGER, geom GEOMETRY);
CREATE TABLE mytable_invalid(id INTEGER, geom GEOMETRY);

Create a rule, that test inserts for validity and shunts the data into the invalid table where appropriate.

CREATE RULE mytable_valid_insert AS ON INSERT TO mytable WHERE NOT ST_IsValid(NEW.geom) DO INSTEAD INSERT INTO mytable_invalid VALUES ( NEW.id, NEW.geom );

Voila! Now you can be guaranteed at all times that the data in your table is valid, and you have an easy place to check for invalid inputs to repair and add to the production table.

In the GeoExt community, we had a long ongoing discussion about whether we want to use inheritance or factory methods. Why bother at all? Let’s say we want to have a data store with records representing the layers of an OpenLayers map. This store should always be in sync with the map’s layers.

Obviously, we will want to register some event listeners on the map to handle adding and removing of layers. And also listeners on the store to handle adding and removing of records.

Let’s think about inheritance for a minute: we just create a subclass of Ext.data.Store with a setMap method for making the store aware of a map and registering the listeners. And while we’re at it, also a unsetMap method that unregisters the listeners when we don’t want to sync with the map any more. Good. Now what if we want to use a different flavor of a store, say a GroupingStore? Well, mh, another subclass of Ext.data.GroupingStore, but with the same functionality. So would we provide a subclass for all subclasses of Est.data.Store? Sounds like a lot of overhead. And what if someone creates a custom subclass? Yet another subclass? There must be easier ways…

What would the factory option look like? We would just create a utility method that creates our store, assigns it a map and registers the required listeners. Optionally, when calling the method, we can pass a store class or instance that we want to use. Very well. The only remaining question seems to be which namespace to use for this method, since we will probably have lots of factory methods in our library. But what if we want to ask the store later which map it belongs to? Well, there’s no map property. We would have to dig into a record and ask a layer what its map is. And what if we want to stop syncing with the map? Well, no luck…

Now an alternative way to do this would be to create something between the map and the store that takes care of the syncing. But then our users come either from ExtJS or from OpenLayers and are used to say map.addLayer() or store.loadData(). They would have to learn about this new “inBetween” object, and we would probably want add and remove functions there that can handle both layers and records. Still a clean solution, but maybe not so easy to understand for developers.

Back to thinking. The “factory” function could modify the object to look like it has a different interface. But this is not good practice, because we are not properly advertising our modifications. Well we could document this in full text in the api docs for the factory method, but would someone look there when debugging unknown code? Unlikely.

More thinking. Extending classes in ExtJS is really easy:

var NewProto = Ext.extend(ExistingProto, overrides);

What is “overrides”? Just a javascript object, a mixin. For extending Ext.data.Store, our overrides object could look like this:

GeoExt.data.LayerStoreMixin = {
    map: null,
    constructor: function(config) {
        arguments.callee.superclass.constructor.apply(this, arguments);
        // our customized constructor code here
    }
    setMap: function(){}
};

Now this is a structure that every apidoc processor for JavaScript recognizes, so we can document the map property and the setMap method inside this object.

Another interesting part here is the first line of the constructor, especially the arguments.callee.superclass thing. Usually when we extend prototypes in ExtJS, we do something like this to call the superclass’ constructor:

NewProto = Ext.extend(OldProto, {
    constructor: function(config) {
        NewProto.superclass.constructor.apply(this, arguments);
        // our constructor code here
    }
});

Obviously this makes our overrides object totally dependent of the new prototype we create with it. And this is where arguments.callee comes in handy. This property always returns the currently executing function.

We can now just add this overrides object to the library, and users can override any store prototype directly in their code in one pass with instantiation:

var myInstance = new (Ext.extend(
    Ext.data.Store,
    GeoExt.data.LayerStoreMixin
))(config)

Of course we could also provide the above as a default implementation, create a file called LayerStore.js, put the documented mixin and the default implementation into that, and we have a nice convention for handling such mixins:

File: CustomClass.js

CustomClassMixin = {
    // mixin definition here
}
// default implementation
CustomClass = Ext.extend(Class, CustomClassMixin);

What’s the take home message of this? We can easily create a new constructor by adding our mixin to an existing prototype. The result is that we keep a sensible prototype chain and always work with instances that have the unaltered interface of their prototype. This approach has made it into the GeoExt LayerStore.

I’m a fanboi, I can’t help it, it’s true, I don’t deny it.

I found this performance test run at tweakers.net really great. Since graphs are worth 1000 words, they don’t waste a lot of time chattering.

First, MySQL throughput under increasing concurrency, for multiple hardware platforms.

MySQL Throughput vs Concurrency

Then, PostgreSQL in the same conditions. Note how PostgreSQL not only achieves higher peak throughput, but maintains that throughput as concurrency goes up.

Now that’s a database I can trust to bring home the bacon.

If you are lucky enough to use OS/X for development and have a somewhat recent version of the operating system and XCode, you have access to one of the nicest profilers around, the Shark.  Unlike the old GNU tool “gprof”, the Shark does not require that you compile your program with special flags in order to profile. This is especially nice for programs with lots of dependencies and functionality in linked libraries — it was very hard using gprof to figure out where the effort was going when cycles spent in libraries were not visible in the profile output.

Unlike gprof, the Shark operates by sampling at the operating system level, so 100% of the activity on your computer is visible to it, and processes don’t need special compilation flags to be available for profiling.

As a result, the user interface for the Shark is dead simple. When you start up the Shark (/Developer/Applications/Performance Tools/Shark) you get this small window.

The default options (Time Profile, Everything) are sufficient for profiling the performance of C/C++/ObjC applications. When you press “Start”, Shark starts sampling. It will sample for about 20 seconds or until you press “Stop”, whichever comes first.

Because the Shark is sampling all processes, you want to build up as much time on your process-of-interest as possible. For Mapserver, that means running as FastCGI, so that multiple map requests can feed the data associated with the mapserv process id. If you run Mapserver in standard CGI mode, you’ll get lots of entries for individual mapserv runs, with relatively few samples taken on each.

When the sample run is finished, Shark spends a few seconds analyzing the data, then pops up a window to browse the results.

At the bottom of the window is a “Process” selector. The process that used the most CPU during the sample run (probably the one you are interested in) is opened by default, but you can look at the profile of any process that existed during sampling by selecting it here.

The default view of results is “Heavy (Bottom Up)” which lists the functions that used the most CPU, with a collapsed tree underneath each row showing the call stack(s) down into that function.  The opposite view is also available, “Tree (Top Down)”, which starts from main() and allows you to explore downwards, seeing which parts of the execution used the most CPU time.  The “Heavy and Tree” view shows both views and once, and slaves them together so you can click on a function in “Heavy” and see where it falls in “Tree” and vice versa.

The columns of the result table are: “Self”, how much time was spent within a function; “Total”, how much time was spent within a function and all functions called by that function; “Library”, the process name or library name the function resides in; and, “Symbol”, the function name itself.

If you tried profiling with gprof but gave up because it was too hard to read the results or too hard to interpret them, it’s time to forget about that experience, mount up and ride the Shark to a faster tomorrow!

I’ve been remiss in not announcing this widely sooner, but now we are sure it’s happening, and hey, better late than never.

We are proud to announce February 28 and March 1st, Columbia Center For New Media Teaching and Learning and OpenGeo are hosting a GIS sprint in New York City.

If you happen to be in the greater NYC area next weekend and are interested in doing some coding and carousing, Columbia University has been kind enough to open one of their spatial labs to have a sprint on Saturday.  The OpenGeo team (who will be coming in from the far lands of Austria, Spain, Canada, Italy, Montana and elswhere) will be hosting at the TOPP offices on Sunday.

For us, it’s an opportunity for us (especially the remote OpenGeo folks) to connect with all the great GIS stuff going on around Columbia and New York.  If you are interested in the opensource projects we are working on, like GeoServer, GeoExt, OpenLayers, etc, or just opensource GIS in general, it’s a great opportunity to come down and hang out and talk shop.

Likewhise, if you’ve got an opensource project in the spatial realm, please c’mon out and code with us!

If you think you might come, join the mailing list and drop us a line so we can arrange for access and ammenities.

-whit

Eric Raymond says that open source software comes from developers “scratching their itch” — an unpleasant bodily metaphor for “solving their own problems” — a side effect of which is software being released and usable by everyone.

My current itch is attacking the “well-known text” (WKT) representation of spatial reference systems. In theory, an Open Geospatial Consortium standard spatial database like PostGIS stores spatial reference information in the SPATIAL_REF_SYS table, and the actual information about reference system parameters is serialized in a “well-known text” string, stored in the SRTEXT column of that table. In practice, what PostGIS really uses for coordinate transformations is a PROJ4TEXT string in a spare column of the table, and the SRTEXT is just window dressing — we carry it, but we don’t use it.

Using the WKT representation directly is attractive, because it drops needless duplication of information and allows more direct interoperation with things like ESRI “prj” files, which are themselves just WKT serializations.  Unfortunately WKT is not as “well-known” as the name would have us believe. Every vendor has used slightly different naming for things like projection operations, parameters, datum names, and so on.

So my itch is multi-fold: I want to be able to parse WKT, I want to learn the technologies necessary to parse WKT (bison and flex), I want to be able to standardize WKT (to strip out the vendor-specific bits) and I want to be able to turn my parsed form into PROJ4 projection objects, because I’ll still be using the PROJ4 engine for transformations at the end of the day. I’m an itchy guy.

So far, I have achieved the parsing and learning-how-to-parse goals, and placed my results in a spike in the PostGIS SVN repository. Next up is standardization, and finally creating PROJ4 objects. Then I’ll try to hook the whole thing into PostGIS.

I’ve been doing a bit of work on Drake recently, improving the layout, adding features, you name it. I got pretty excited about all the sweet stuff I had gotten working (in no small part due to Seb’s CSS-fu) and so I turned to one of my less techie colleagues and said “Hey, check this out!”

The first words from her mouth? “Why does Montana have two labels right next to each other?”

Tiles rendered one at a time, with one tile outlined

The answer’s pretty simple: when grabbing the rendered images from a WMS server, OpenLayers uses 256×256 pixel tiles. The advantage is that if the user pans a little bit, OpenLayers can reuse the tiles that were already on the screen (using tiles also enables us to use caches that implement the WMS Tiling Recommendation from OGC, but that’s another post). There is a drawback, though: GeoServer only takes into account the part of the map that is actually being rendered when deciding where labels should go. That means that if a feature straddles a tile border, it’s very likely that it will end up having duplicate labels close to each other.

Fortunately, there is a nice trick to help with this, called ‘metatiling.’ Instead of rendering the 256×256 pixel graphic that’s requested, the server renders a cluster of tiles centered on the requested tile, doing the layout logic as normal. Then, everything aside from the requested pixels is cropped off, so you get back the tile you wanted, but with the renderer aware of the surrounding space. The result? Montana only has one label

GeoServer rendering tiles in context

You may, at this point, be wondering, if metatiling is so great for layout, why doesn’t GeoServer use it by default?  There are a few reasons why it may not be the best option in all cases:

• Referencing. The main reason that GeoServer doesn’t metatile by default is that, given only a tile’s bounding box, GeoServer is oblivious as to where in the metatile it is.  GeoServer requires tile requests to give it a reference point to use for the origin of the tile grid for disambiguation purposes.
• Performance: Since GeoServer is rendering multiple tiles for each request, the memory usage is multiplied as well.  GeoServer uses a 3×3 metatile, meaning each request takes roughly 9 times as much memory to handle.  GeoServer does cache the entire metatile for a short time, but in general, larger metatiles mean better rendering at a fairly high memory cost.
• Lack of guaranteed results: Depending on your data, you may still have large polygons that sit on the bounds of metatiles.  Metatiling does help with label layout, but it’s still subject to the same fundamental problem: only a part of the map is taken into consideration for label positioning.  A careful analysis of where your features lie with respect to tile bounds may help, but I’m not aware of any tools to help with such analysis.

Despite all these issues, metatiling is usually a good visual improvement.  To enable it in GeoServer, you just need to append ‘&tiled=true&tilesorigin=[x,y]‘ to your requests.  In OpenLayers, you do that by something along the lines of:

layer = new OpenLayers.Layer.WMS(name, this.wmsUrl, {
layers: "theworld",
format: "image/png",
tiled: "true",
tilesOrigin: "-180,-90"
}, {
isBaselayer: false
});

Metatiling can also be done with no support from the WMS server at all by WMS tile caches such as Tilecache or GeoWebCache.

One pinned and one detached popup coexisting

The Ext.Window component provides a nice way to put some widgets in a frame that the user can move around your application. Moving that window from one parent component to another, or from a component to a ‘plain’ HTML element, proves a bit less straightforward than one might hope.

I found this out recently when Sebastian Benthall and I were working on some improvements to his popup code for GeoExt.  The main feature we really wanted to get working was detaching popups from the map.  I’ll elaborate a bit: the popups in question show up when the user clicks on a particular feature in the map.  Initially these popups are pinned to the map, meaning they have a nice tail pointing to the feature they came from, and follow that feature around if the user navigates around the map.  However, if you want to, say, have an editable map, it would be useful to be able to detach those popups from the feature so you can move vertices around while editing the feature data.  We initially tried simply removing the popup’s Ext component (an Ext.Window) from the map container and adding it to the body, but, as it turns out, not all layouts support dynamically removing and adding components after the layout has been rendered once.  The body element isn’t even an Ext Component (and therefore lacks an layout object), so we had to go with a less straightforward solution.

We thought of just making a new Window and copying the contents of the original to the new one so that Ext’s layout code would have a fresh component to layout, but there doesn’t seem to be a good way to get all the contents (there is Component.cloneConfig(), which initializes a new Component with the same parameters originally passed to the constructor, but that doesn’t take care of changes in state on the contained widgets, so it doesn’t really do what we needed.)  So, we grabbed our headlamps and pickaxes, and started digging beneath the documented ExtJS API.

Once an Ext.Window has been laid out, it keeps a reference named el to the Ext.Element that it draws itself to.  It also maintains a reference named container to the Ext.Element that is the parent of that element.  container is used for things like building the proxy object that moves with the mouse while dragging the window, and is a private property specific to the Window object. The container is set when the window is rendered for the first time, which is normally triggered by calling doLayout on it or a component that contains it. However, the rendered HTML element is cached afterward, so subsequent layouts don’t call the window’s render method.  Therefore, the layout code that handles things like:

panel.add(window);
panel.doLayout();

fails to update the container, resulting in some fairly funky behavior: the window is rendered in the new panel, but continues to interpret drags relative to its original container.  This results in all moves being offset by the difference between the top left of the original container and that of the new container.  We can fix this by updating ‘container’ manually:

window.container = panel.el;

Things are a little different when we want to detach a window from any component at all.  In this case, we need to move the window to the body of the html document, which isn’t a component at all!  We have to do all the layout work manually:

panel.remove(window, false); // remove() destroys the component unless false is passed in
panel.doLayout();
window.el.appendTo(Ext.getBody());
window.container = Ext.getBody();

In either case, changing the container causes a couple of visible glitches for the user.  One is that the window moves visually (since its position is specified relative to its containing component); the other is that the window’s shadow is focused remains in the original container.  Conveniently, once the container is hooked up properly, the shadow can fix itself automatically when the Window’s hide() and show() methods are called.  There is also a pretty simple trick to fixing the window’s absolute position; just use Window.getPagePosition() and Window.setPagePosition() to ensure the page-relative coordinates are the same before and after the switch.  So, to move a Window to a new component:

var pos = window.getPagePosition();
panel.add(window);
panel.doLayout();
window.container = panel.el;
window.setPagePosition(pos.x, pos.y);
window.hide();
window.show();

or to the page body:

var pos = window.getPagePosition();
panel.remove(window, false);
panel.doLayout();
window.el.appendTo(Ext.getBody());
window.container = Ext.getBody()
window.setPagePosition(pos.x, pos.y);
window.hide();
window.show();

Edit: Fixed a broken link and corrected some mistakes in the code examples. Thanks, Claude!