When building a geospatial application it is very easy to get into problems with performance, mostly because of the volume of data - complex geometries containing many vertices and/or large datasets containing many geometries. Just think of the number of objects that Google Maps has to store to support everything we do while using the application - all the buildings in the world, all the roads, rivers, terrain information, traffic… On the other hand, users will often expect instant or close to instant responsiveness of the application while moving on the map and zooming in and out. So, how do we accomplish that?
There are a few different techniques that can dramatically improve the performance of spatial operations. Before we dig into the solution, let’s analyze the problem. Spatial queries typically involve intersection or similar operations and not the exact equality comparison that we often use with regular SQL queries. If we have a large dataset of complex geometries, it will be very computationally heavy to check which of the geometries intersect with the given geometry.
As an example, let’s analyze a problem that is faced by almost any geospatial application, such as Google Maps. When a user opens the map, it is up to the system to determine which of the objects from the database of all objects on the globe are relevant for the current AOI (area of interest - the rectangle the user is looking at). This is typically done by intersecting the AOI rectangle with polygons (or other vector objects) representing the features of the map to determine which ones should be fetched from the storage and rendered on the screen. To simplify, we can assume that we are only dealing with rivers and those rivers are stored as records of PostGIS’ geometry (Polygon,4326) type in a separate table.
CREATE TABLE river (
rid uuid PRIMARY KEY,
name text,
shape geometry(polygon,4326)
);
Then we can try to find those rivers that intersect with our AOI polygon using the query below.
SELECT *
FROM river r
WHERE ST_Intersects(r.shape, ST_MakeEnvelope(24, 47, 25, 48, 4326));
We are using ST_MakeEnvelope function to construct our AOI rectangle with the provided top-left and bottom-right latitude/longitude coordinates and 4326 as the SRID of the polygon.
As expected, running the “Explain Analyze” tool in Postgres shows that we are indeed performing a table-scan on the “river” table, which can slow things down if we have more than a few records in the table.
Spatial Indexing
Enter spatial indexing. The first and easiest optimization that we would implement is creating a spatial index on our river table. Spatial indexes are explained in more detail in the first article in the series.
CREATE INDEX river_shape_idx ON river USING gist(shape);
If we re-run the SELECT query above, the query execution plan switches from “Seq Scan” to “Index Scan” and query execution time is reduced by ~97%.
In some cases, this may be perfectly sufficient and we could stop with the optimization here. The spatial index helps to narrow down the records from the table whose bounding boxes intersect with our AOI rectangle, so PostGIS needs to perform the intersection math on only a small fraction of the polygons in the table instead of the whole table. When the polygons we are dealing with are relatively simple and small, that may be sufficiently performant. Most of the intersecting bounding boxes will in fact result in actual polygon intersections and there will be very few false positives.
As an example, we can consider Google Maps and polygons that represent the buildings (gray blocks on the screenshot below). They are relatively small and simple and there’s very little difference between the actual polygons and their bounding boxes, which minimizes the chance of a false positive index match. In other words - if our AOI polygon intersects with the bounding box of a building polygon, it’s very likely that it would also intersect with the actual building polygon. That means that spatial index will be very performant when applied to building polygons.
Polygon Subdivision
Now let’s consider the opposite case from simple building polygons - a major river that runs across the continent. The polygon that represents the river will be very large and contain many, many vertices. That means that the bounding box of that polygon will also be large, spreading across the continent, and wherever we place our rectangle, it would likely intersect with that bounding box. However, the actual area of the river will in most cases cover just a tiny fraction of the total area of its bounding box. That means that the probability of the river actually intersecting with a rectangle intersecting with its bounding box will be very low. And that is a problem. The spatial index uses bounding boxes of geometries to index them and after the crude intersection is established, using the spatial index, PostGIS will go into the detailed calculation to verify if the two geometries truly intersect. That means that in most cases, the spatial index will not be doing its job and PostGIS will still take more time to compute intersections.
Considering the example below, the river Danube runs across Eastern Europe and has a massive bounding box. It’s very easy to fool the spatial index and get a false positive match for an AOI rectangle that can be very far from the actual river polygon. If we need to find all rivers crossing the green AOI rectangle, the spatial query will consider the Danube, even though it’s hundreds of kilometers away from the AOI. PostGIS will realize that AOI does not intersect with the river only after performing the costly calculations.
The solution to this problem is to use a polygon subdivision that chops up large, complex polygons into simpler ones, store those polygon segments in a separate table, and use them to determine intersecting objects. Being simpler and physically smaller, those polygon segments will have much smaller bounding boxes, so chances of false positives will be dramatically reduced. And even in cases when a false positive happens, actual intersection calculation is much faster on a smaller, simpler polygon. That way, we are minimizing the area on the map that is part of the bounding boxes of the object but is not actually part of the object.
CREATE TABLE river_segment (
sid uuid PRIMARY KEY,
rid uuid REFERENCES river(rid),
name text,
shape geometry(polygon,4326)
);
CREATE INDEX river_segment_shape_idx ON river_segment USING gist(shape);
Now we can construct the segments and copy them to the new table. The maximum number of vertices is set to 10 as it offers a decent balance between polygon complexity and the resulting number of segments. However, further experiments with the actual use-case may reveal a slightly different number of vertices as an optimum.
INSERT INTO river_segment
SELECT uuid_generate_v4(), rid, ST_Subdivide(shape, 10)
FROM river;
Note: the step above can also be set up as a DB trigger, to automatically replicate segments as the new rivers are added.
If we now run our original query on the river_segment table instead, we would get an even better execution time, approximately a tenth of the time we got using the spatial index on the full polygons.
SELECT *
FROM river_segment r
WHERE ST_Intersects(r.shape, ST_MakeEnvelope(24, 47, 24, 48, 4326));
If we need to reconstruct the portion of the original polygon that intersects the AOI, we can run ST_Union function against the segments within the AOI.
SELECT r.rid, ST_Union(r.shape)
FROM river_segment r
WHERE ST_Intersects(r.shape, ST_MakeEnvelope(24, 47, 25, 48, 4326))
GROUP BY r.rid;
Alternatively, we can use the ST_Intersects on river_segment table to determine which of the rivers intersect with our AOI, but then run another query to get the visible segment of the river from the original table. Depending on the use-case, we might want to grab the whole original polygon of the river or just the intersection with the AOI. The query below returns both.
WITH intersected_rivers AS (
SELECT DISTINCT r.rid
FROM river_segment r
WHERE ST_Intersects(r.shape, ST_MakeEnvelope(24, 47, 25, 48, 4326))
)
SELECT r.rid,
r.shape as original_shape,
ST_Intersection(r.shape, ST_MakeEnvelope(24, 47, 25, 48, 4326)) as clipped_shape
FROM river r
INNER JOIN intersected_rivers i ON i.id = r.id;
Using the same example as above, if we were to subdivide the polygon representing the Danube into smaller segments, we would get a lot of purple bounding boxes that cover just a fraction of the original bounding box of the whole river. Given the same AOI rectangle as before, we can easily see that it will intersect no bounding boxes of the segments of the Danube. Therefore, the spatial index will be able to rule out Danube as a candidate for intersection with the AOI.
As always, there are costs to this optimization. Spatial indexes, like regular indexes, have a negative impact on the performance of INSERT operations. With the subdivision, INSERT operations will become even slower as they include the subdivision operation and storing the segments into a separate table that also has an index. Also, we are (more than) doubling the storage size for polygon objects because we still need to keep the original objects along with the segments. In our test case with major rivers of the world, the number of segments required to cover the polygons of the 1500 rivers is over 300k which requires approximately 5 times more disc space for storage.
Depending on the use case, the performance gains of intersection queries are significant and often outweigh the costs. In most cases, we will query the database much more frequently than perform insert operations into the database, so it makes sense to optimize for reads.
Index-only Queries
Probably the fastest way to get information out of a spatially indexed PostGIS table are index-only queries. Using the && operator, we are instructing PostGIS to return any geometries whose bounding boxes intersect with the bounding box of the provided geometry. This can be determined purely by examining the R-tree index without validating if the actual geometries intersect, and is inherently a very fast way of querying spatial data.
SELECT *
FROM river r
WHERE r.shape && ST_MakeEnvelope(24, 47, 24, 48, 4326);
However, index-only queries come with a warning. In the case of large and complex geometries, such as rivers or roads, the query will potentially return many more potentially heavy geometries that do not necessarily intersect with the AOI. In the example above, for the green AOI rectangle, we would get the whole polygon for the Danube river back in the results even though it’s nowhere near the AOI. This will slow down data retrieval and would unnecessarily load down the frontend components. In those cases, the downside of index-only queries is probably not worth the gain in the query speed. On the other hand, when objects are relatively simple, such as building footprints, speed increase often outweighs the extra load coming from false-positive matches. In those cases, index-only queries are a great way of improving system performance.
Polygon Simplification
Now let’s consider the following use case. We are looking at the map that shows the Balkans peninsula and our intersection query determined that the Danube should be returned to the client for rendering. The polygon that depicts the whole river of Danube is very complex, as it contains high-resolution details that show every bend and can be used to render the river when zoomed in to the maximum zoom level. That requires many thousands of vertices that capture all the details across the length of the river. However, when zoomed out, the whole river is a couple of hundreds of pixels long, so do we really need all those details? The answer is no. Every visible pixel of the river rendered zoomed out may contain many vertices which means that we would be returning much more information than is necessary to render the polygon.
The same applies to other large polygons, such as islands, major roads, etc. We can test this through Google Maps. As we zoom in and out, Google Maps reveal and hide details of complex polygons dynamically.
The solution is to leverage the polygon simplification operation that reduces the number of vertices of a polygon with a given tolerance. The output of the operation is a polygon that has no two vertices closer than the provided tolerance. The image below shows the same complex polygon representing the UK island on the left which is then simplified with increasing tolerance going to the right. The right-most polygon looks nothing like the original, but it is very light, containing only 8 vertices. Our goal is to find a good balance between object size and how well it represents the original polygon.
In order for the polygon simplification to work, we need to know the zoom level currently viewed by the client, so zoom level needs to be added to the parameters of the call. Zoom levels are typically enumerated between 0 and 23, 0 being the whole earth, 3 approximately the size of a continent, going down to zoom level 15, showing individual buildings. Extreme zoom levels over 20 show the map with a very high level of detail.
The next thing we need to know is the approximate resolution of the map in each of the zoom levels. Getting the exact number is not easy, as it depends on the latitude, but getting an approximation is good enough. With each zoom level, we are halving the resolution, meaning that the pixel size of each next zoom level is one-half of the pixel size of the previous zoom level. Knowing that the approximate pixel resolution for a given zoom level can be calculated as
pixelSize = basePixelSize / (2 ^ zoomLevel)
Where basePixelSize is the pixel size at zoom level 0 and that is approximately 78km. Depending on the SRID used, we need to convert this to the units of the SRID. In the example above we used 4326 which uses degrees as units. We can convert that 78 km to degrees approximately by dividing by 111km/px, resulting in 0.7 degrees per pixel at zoom level 0. Then our formula becomes
pixelSize = 0.7 / (2 ^ zoomLevel)
Now to tie it all together, let’s add polygon simplification to the query that does simple intersection.
SELECT *
FROM river r
WHERE ST_Simplify(
ST_Intersects(r.shape, ST_MakeEnvelope(24, 47, 25, 48, 4326)),
0.7 / (2 ^ zoomLevel),
false
);
The last parameter of ST_Simplify signals to PostGIS not to return geometries that collapse into a single point, as there’s nothing to render for them. That means that as we zoom out, we will see fewer and fewer details and at some point, the smaller objects on the map will completely disappear. Google Maps does this as well. If we navigate the map to a small island on the sea, we can observe all the details when zoomed all the way in. But as we zoom out, we will see fewer and fewer details, and eventually, the island will disappear completely.
We can take it a step further and increase the tolerance some more to further reduce the complexity of polygons without sacrificing the rendering quality too much. Suppose that we do not want to render polygon edges shorter than 5px, we can multiply the tolerance in the query above with 5, which will further reduce the complexity of the polygon. Depending on the use-case, choosing a value between 1 and 10 pixels is a good starting point.
We can apply polygon simplification in two ways:
1. On-the-fly - every time users make a request, we determine the level of detail needed for the requested zoom level and simplify the polygons before returning them to the client
2. Pre-computed - knowing the range of zoom levels, we can pre-calculate and store all simplified versions of original polygons and then just retrieve them according to the requested zoom level.
As before, the tradeoff here is between the performance of the read operations on one side and the performance of insert operations and storage space on the other. Depending on the use case, it may make more sense to take one approach over the other.
Conclusion
In this article, we got acquainted with spatial indexing, polygon subdivision, and polygon simplification techniques. When applied properly, these techniques can dramatically improve the performance of geospatial applications, dealing with large numbers of geometries or dealing with complex, computationally heavy geometries. There are many more ways to deal with spatial optimization, but hopefully, this article scratches the surface.
About the author
Branislav Stojkovic is a senior software engineer with more than 15 years of industry experience working at our engineering hub in Niš.
Branislav has extensive hands-on experience in designing, building, and maintaining large enterprise systems, mostly in the financial industry, and strong technical knowledge in back-end development using Java and .NET technologies and database modeling. The financial sector is often volatile and requires a deep understanding of the domain knowledge, and Branislav can dive deep into the details but also have a broad, high-level view of the whole system and the roadmap for the future. Agile methodologies are often employed to cope with the volatility of the financial industry. This search for optimal solutions started at an early stage of his career while working with algorithm design and competing in high school and university competitions.
