10

I am currently working in the field of isochrones and the underlying algorithms. What now causes problems is not the calculation if the isochrone itself, but the visualization of the results.
The result of my isochrone algorithm are points and edges. In fact I do have a working solution, but for 3873 edges and for 1529 nodes things seem to take forever (around 2.0 seconds on my Lenovo T440s laptop which contains a 2015 Core i7 CPU and a pretty fast SSD). Instead of seconds I want something more like msec :-).

Maybe someone can help me to reduce the calculation time needed to build the polygons which visualize the reachable areas.

But wait... first things first!
Here is a visualization of the edges that I are the calculation result of my isochrone: Isochrone calculation result (skeleton existing of linestrings) These edges are stored in a PostGIS database table and are simple linestrings.

What I want to show to the user looks like this: enter image description here Note the disconnected areas in the very south and very east of the picture. These should be drawn as separate areas (so no merging allowed here :-))

Currently I am using this query:

SELECT ST_AsGeoJson(St_Transform(ST_Multi(ST_Collect(polygons)), 4326)) AS coverage FROM (
    SELECT ST_MakePolygon(ST_ExteriorRing(ST_GeometryN(segments, generate_series(1, ST_NumGeometries(segments))))) AS polygons FROM (
        SELECT ST_Union(ST_Buffer("GEOMETRY", 20, 'quad_segs=2')) AS segments FROM my_edges AS a
    ) AS b
) AS c

I already did some experimenting and I also read a lot of documentation, but I just can't find a better solution.
In my eyes the big problem is the usage of ST_Union (as stated in the docs this function can be slow). The very interesting thing is that replacing it with ST_Collect seems to slow down the ST_Buffer calculation so that all-in-all the following query takes even longer, although it does not fill the areas between the edges (it only creates a buffer around the lines):

SELECT ST_AsGeoJson(St_Transform(ST_Multi(ST_Collect(polygons)), 4326)) AS coverage FROM (
    SELECT ST_Buffer(ST_Collect(ST_LineMerge("GEOMETRY")), 20, 'quad_segs=2') AS polygons FROM my_edges AS a
) AS b

This takes around 3.8 seconds on my system (so nearly twice the time). My first conclusion out of this little benchmark is that ST_Buffer gets unexpectedly slow when it comes to MultiLineStrings (even slower than when creating buffers for each line and merging the buffers - which in my eyes is just weird)

I also tried to use alpha-shapes (using the implementation from pgRouting), but since there is no alpha value to set (and in fact I would not really now to which value to set such a value) I just get one great polygon (so I'd lose the regions in the very south and east as separate regions which is not what I want).
Also ST_Polygonize (which was the first thing that came in my mind) did not produce any usable results, but maybe I missed something here...

Is there a better way to create the area shown in PostGIS? Maybe also by using java code (jts) or client side javascript code (jsts)? In fact I could live with loosing some detail as long as the areas shown in my result stay separated and the calculation gets (much) faster.

  • Can you not just use ST_Exteriorring(ST_Dump(ST_Union(ST_Buffer(geom,....))).geom seeing as anything buffered will be a polygon anyway, and ST_Union will connect all intersecting geometries, so no need for MakePolygon or GeometryN. You might need to test for Linestrings that sometimes result from ST_Union after a buffer, but that is easy with ST_GeometryType(geom). As far as using Java or jsts is concerned, you can, but it is unlikely to be faster, given that a large part of the Postgis (GEOS) functions are C/C++ ports of JTS in the first place. – John Powell Jul 13 '15 at 16:58
  • You are right, this works, but in fact it is not faster (takes ~3.1secs, while using GeometryN takes 2secs). Here is what I used: SELECT ST_AsGeoJson(ST_Transform(ST_Exteriorring((ST_Dump(ST_Union(ST_Buffer("GEOMETRY", 20)))).geom), 4326)) FROM my_edges; – Nikolaus Krismer Jul 13 '15 at 17:16
  • @john-barça: Oh.. I mist the quad_segs=2 part in the ST_Buffer when trying your approach... with that change the queries are even (both at around 2secs). However, this still is very slow (in my eyes), is there another way to try it? – Nikolaus Krismer Jul 13 '15 at 18:49
  • Interesting problem....do you want to share some test data? – dbaston Jul 15 '15 at 1:56
  • If it helps I am happy to share some data. All of the things I do here are open-source, so this should not be a big problem. First thing to notice: A web application for testing is located at dbis-isochrone.uibk.ac.at:8080/testing. More information about the things I work on can be found at dbis-isochrone.uibk.ac.at. In the "links" section of the website there are some further references (including some test data) – Nikolaus Krismer Jul 15 '15 at 8:32
5

Setting aside the GeoJSON serialization, the following takes about 6.3 seconds on my laptop:

SELECT
  ST_MakePolygon(
    ST_ExteriorRing(
      (ST_Dump(
        ST_Union(
          ST_Buffer(geom, 20, 2)))).geom))
FROM bz_edges

Looking at the data in OpenJUMP, I noticed quite a bit of detail in the street segments, relative to the desired level of detail in the output. It seems that even an on-the-fly simplification of these lines can produce a big speedup in PostGIS:

SELECT
  ST_MakePolygon(
    ST_ExteriorRing(
      (ST_Dump(
        ST_Union(
          ST_Buffer(ST_Simplify(geom, 10), 20, 2)))).geom))
FROM bz_edges

which brings things down to 2.3 seconds. I thought I could maybe do better by storing the generalized geometry in a separate column, instead of calculating it on the fly, but that actually provided no additional benefit.

Depending own how much code you're willing to write, you can almost certainly do better in Java, if nothing else because you can take advantage of multiple cores. (For what it's worth, JTS performs the above operation in 2.8 seconds). One approach might be to extend CascadedPolygonUnion to make some of the union operations happen in parallel. (update - here is a ParallelCascadedPolygonUnion)

I noticed in the sample data that the edges are stored with references to their start and end nodes, i.e., you have a pre-built graph. I suspect you can generate these polygons much more quickly if you work from the graph instead of using generic geometry operations. For example, I think you could so something like this:

  1. identify the connected components of the graph
  2. for each connected component, find the node with the minimum X coordinate (guaranteed to be on the outside of the component)
  3. walk the edges of the component, always turning left (or right) when possible. This will give you the exterior ring of each component.
  4. polygonize the exterior ring and buffer appropriately.
  • Thanks... simplification is a great and even "simple" improvement. It took the time needed on my laptop down to 1.5sec. It is not where I want to be, but a bit better. – Nikolaus Krismer Jul 15 '15 at 13:46
  • Regarding your suggested solution (points 1-4). Sounds also very simple and is worth a try. I thought of something similar, but I am stuck at point1 (so very early :-)). How would one identify connected components (only thing I can think of is a recursive query which might also be very slow). – Nikolaus Krismer Jul 15 '15 at 13:48
  • @NikolausKrismer I use both JGraphT and loom for tasks like this. If you write your own graph methods instead (not a bad idea for best performance), a depth-first search will find you the components. (You could find them in the upcoming PostGIS 2.2 with ST_ClusterIntersecting but I think you'll want any kind of graph processing to happen outside of the database anyway, so this probably isn't useful). – dbaston Jul 15 '15 at 14:10
  • these are some great hints. I looked at JGraphT and this could certainly help solve my problem. However, I also looked at Postgis 2.2 and the ST_ClusterIntersecting function -> it takes about 200-250msec to identify the different clusters in the case above. That's OK to me (JGraphT could certainly do better). Now I have to deal with creating the exteriorRing (ST_ExteriorRing fails, since ST_MakePolygon says my linkes are no shell) – Nikolaus Krismer Jul 16 '15 at 9:06
  • I see two complications: (a) you need not only the exterior ring, but also any segments that extend outward from that ring, and (b) it looks like your lines don't actually intersect at some intersections. You'll need to fix (b) if you're going to try to construct a geometry from the results of a graph walk. – dbaston Jul 16 '15 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.