New vector isochrones generation

[Laurent GRÉGOIRE]

Hi,

I'm in the process of implementing a new vector isochrone polygon
generation tool and here are a description of an early version for
people that would be interested.

The current vector isochrone generator rely on building a convex hull
out of a delaunay triangulation of some vertex points from the
shortest-path tree graph. However, this implementation lack a few
features. I will be trying to address the following points:

a) support disconnected areas (the output of the isochrone should be a
set of polygons). This is important especially for transit isochrones
where the result is a highly fragmented area with small islands around
stations;
b) support holes. This is less important than a) but for eg. rivers or
large unaccessible areas, some polygons can have many holes;
c) no hard limit on the covered area. Real accessible area should be
covering an area nearby roads, with an accessibility function (off-road
walk speed, off-road max distance...);
d) output is "street-accessible area". There can be two definitions of
an isochrone (only considering area accessible on street, or accessible
even if you are still on the transit vehicle), for standard usage I
think the first definition make more sense;
e) be able to build several isochrones per request;
f) have a configurable and constant precision;
g) polygons within each isochrone should not intersect (also holes);
h) isochrones from different z values should not intersect.

Since some points would be rather hard or impossible to solve with the
current isochrone implementation, I made another one based on a rather
different algorithm, which try to cover all this points.

Broadly speaking, it uses the current analyst sample factory to get a
z=f(x,y) time (or whatever) function, and use an iterative grid sampling
with interpolation to determine the set of isochrone rings. It then
build a set of polygons (some with holes) out of it. See attached teaser
isochrone_1.jpg to see the results. The method in details:

1) Configure the isoline builder with a center point and a grid size
(unit) determining the resulting precision. Try were made with a 100m x
100m grid.

2) Seed the covered area by determining the set of 4x4 units "big-tiles"
that covers all the relevant point from the SPT. The covered area at
this stage can be non contiguous. See debug_1.jpg to see the seed points
and the resulting seed tiles set. This size (4x4) determine the
precision in small islands/holes detection.

3) Recursively split each cutting edges from the big tiles (a cutting
edge AB is an edge intersecting the z0 plane where z0 = isochrone time,
that is zA<z0<=zB or zA>=z0>zB), until the edge has length 1 unit.

4) Expand small cutting edges: for each cutting edge, explore the 2
adjacent tiles by looking at 6 other edges. Use the painter algorithm.
See debug_2.jpg to see the resulting small grid (initial edges and
cutting edges in red, isoline in blue).

5) Build polylines by walking around cutting eges. Handle tiles with 4
intersecting edges (saddle point) by an heuristics. Walk by keeping
always Z>z0 the same side. Determine exact position by interpolating
between edge end-point based on the z values. To better handle infinite
z values (which appear quite often), we further divide the edge to get a
better approximation. See grid_1.jpg.

6) Build polygons out of rings. Shells are CW polylines, holes are CCW
polylines. See isochrone_2.jpg for more details to see islands and
holes. Isochrones from different z values should never intersect.

Here is the first commit for those interested in the details:
https://github.com/laurentg/OpenTripPlanner-plannerstack/commit/922421ab41145eb52da71bb6eaebc7adf907eecc

I split the code in three sections (webservice+json conversion, SPT
computation, isoline computation). The isoline builder is in the
geometry package, as it does not rely on any OTP concept (it's a generic
z=f(x,y) isoline builder).

For an intermediate size city on a rather old laptop, the SPT takes
around 300ms to build and the 3 isolines around 500ms (result seen on
isochrone_1.jpg). Time is of course highly dependent on the covered area
and tile size, and a bit on the complexity of the isochrones.

Most of the time is took by the z sampling (by using a fast Z function
we decrease the time to ~20ms -- see linear_zfunc.jpg). The algorithm
try hard to minimize the number of sampling (using a coarse grid,
walking only when necessary, keeping sample values in cache...). With
properly defined and well-behaved z functions, the resulting output
looks quite good even by using a coarse grid (see again linear_zfunc).

The remaining issues are:
1) properly handling infinite z values, which are not always easy to
work with as we rely on linear interpolation.
2) removing the current z sampling method. The current mechanism is not
optimal, as we rely on a spatial index to determine the vertex. This
takes time, and we have all we need at hand in the first place.

I will try to come up with a solution to this two problems, by using
directly z values from the SPT and it's geometry to seed a grid of
sampling points, handling off-road values with a decay function (dz,w) =
f(d) and a proximity "brush" mechanism, adding a w "probability/weight"
factor alongside z. Probably more info soon.

Does anybody would be interested in that implementation? Anyway feel
free if you have any comments.

--Laurent

[Andrew Byrd]

On 10/08/2013 11:06 AM, Laurent GR�GOIRE wrote:
> I'm in the process of implementing a new vector isochrone polygon
> generation tool and here are a description of an early version for
> people that would be interested.

Excellent work Laurent, and as always thanks for your detailed
explanations. This is exactly the kind of thing we need to get OTP
spatial analysis development moving again.

For a much less sophisticated alternate vector isochrone implementation,
also see my branch simpleIsochrone which I just merged into master. It
is intended for long distance travel where OSM detail is not adequate
(it was used for making isochrones from an unofficial China National
Railways GTFS). Decay is a linear function of the distance to the
nearest transit stop, and rather than using interpolation on regular
grids it works "outward" from the vertices.

> Since some points would be rather hard or impossible to solve with the
> current isochrone implementation, I made another one based on a rather
> different algorithm, which try to cover all this points.

I agree that all the features you listed would be important additions to
the vector isochrone service, and justify a new line of development.

> I split the code in three sections (webservice+json conversion, SPT
> computation, isoline computation). The isoline builder is in the
> geometry package, as it does not rely on any OTP concept (it's a generic
> z=f(x,y) isoline builder).

Glad to see you're keeping this modular, which should enable
experimentation by other developers.

> For an intermediate size city on a rather old laptop, the SPT takes
> around 300ms to build and the 3 isolines around 500ms (result seen on
> isochrone_1.jpg). Time is of course highly dependent on the covered area
> and tile size, and a bit on the complexity of the isochrones.

Pruning the search at a certain time (say t0 + 2 hours) often
drastically reduces SPT construction time. At least in the case of
raster isochrones we are not displaying points beyond a certain time, so
the time spent exploring those far reaches of the graph is wasted. For
testing purposes you may also find the EarliestArrivalSPTService useful.
It prunes the search as described and does not handle Pareto co-dominant
states so it can be much faster.

> 2) removing the current z sampling method. The current mechanism is not
> optimal, as we rely on a spatial index to determine the vertex. This
> takes time, and we have all we need at hand in the first place.

Note that the intermediate results of the admittedly inefficient spatial
index lookups (the relationships between grid points and road segments)
can be cached as they are in the raster isochrone endpoints, so the
spatial index only needs to be hit once if you are constantly re-using
the same grid and have plenty of memory.

> I will try to come up with a solution to this two problems, by using
> directly z values from the SPT and it's geometry to seed a grid of
> sampling points, handling off-road values with a decay function (dz,w) =
> f(d) and a proximity "brush" mechanism, adding a w "probability/weight"
> factor alongside z. Probably more info soon.

It would indeed be interesting to try out the reverse approach and work
outward from the street geometries in a "brush like" fashion to fill in
the grid.

> Does anybody would be interested in that implementation? Anyway feel
> free if you have any comments.

I am certainly interested in seeing vector isochrone support continue to
evolve. Please do continue down this path!

-Andrew

[Andrew Byrd]

On 10/08/2013 11:06 AM, Laurent GR�GOIRE wrote:

> The remaining issues are:
> 1) properly handling infinite z values, which are not always easy to
> work with as we rely on linear interpolation.

When we first began work on OTP isochrones, before I decided to adopt a
raster approach, I was considering Delaunay triangulation on a set of
points taken from the network itself -- the vertices, whose z values are
known from the SPT.

When working through the problem I realized that linear interpolation is
simply not correct when known points are selectively located on streets.
The Z surface is defined within M meters of any road segment as the time
to travel the perpendicular distance to the nearest road segment plus
the distance along the road segment to the endpoint vertex with the
smallest travel time, plus the time to reach that vertex from the search
origin. Taking a Y-shaped road as an example, with the search origin at
the base of the Y, the linearly interpolated value half-way between the
tips of the Y's branches will be equal to the value at the tips of the
Y's branches. This clearly does not match the expected Z value; an
isochrone created this way could include a straight horizontal line
between the two branches.

Of course this is the problem you are hoping to solve by generating seed
values based on distance from the street.

Another totally different approach would be to simply directly construct
isochrones around each relevant street segment and union them all. These
would be trapezoids with (approximately) circular endcaps, or the
intersection of two such trapezoids coming from opposite ends of the
segment. Most street segments could be skipped if they are not within a
certain threshold of an isochrone boundary.

-Andrew

[Andrew Byrd]

On 10/08/2013 12:14 PM, Andrew Byrd wrote:
> Another totally different approach would be to simply directly construct
> isochrones around each relevant street segment and union them all. These
> would be trapezoids with (approximately) circular endcaps, or the
> intersection of two such trapezoids coming from opposite ends of the
> segment. Most street segments could be skipped if they are not within a
> certain threshold of an isochrone boundary.

On second thought, it would be difficult to skip some segments and make
assumptions about what's inside or outside an isochrone... but many of
the street-level access polygons would be simple rectangles.

-Andrew

[Laurent GRÉGOIRE]

On 08/10/2013 11:51, Andrew Byrd wrote:
> Pruning the search at a certain time (say t0 + 2 hours) often
> drastically reduces SPT construction time. At least in the case of
> raster isochrones we are not displaying points beyond a certain time, so
> the time spent exploring those far reaches of the graph is wasted. For
> testing purposes you may also find the EarliestArrivalSPTService useful.
> It prunes the search as described and does not handle Pareto co-dominant
> states so it can be much faster.

Thanks for this tip. Indeed I'm already pruning the tree by setting
worstTime. By the way I'm not sure what a "safe" value would be (is
setting max(z0) safe?). Anyway, I'm currently focusing on optimizing the
isochrone construction itself. SPT building is quite independant, I'm
timing both separately.

> Note that the intermediate results of the admittedly inefficient spatial
> index lookups (the relationships between grid points and road segments)
> can be cached as they are in the raster isochrone endpoints, so the
> spatial index only needs to be hit once if you are constantly re-using
> the same grid and have plenty of memory.

I will try not to rely on this caching. The problem with caching is that
you keep values in memory between requests, so if one make one request
per week you end up consuming lots of memory for a very limited scope.
Also the amount of samples can be quite large. Anyway, for the isochrone
computation we generate everything during one request, thus I do not
have the same constraint as for tile generation, which need to
reuse/cache data across different requests.

--Laurent

[Laurent GRÉGOIRE]

On 08/10/2013 12:20, Andrew Byrd wrote:
> On 10/08/2013 12:14 PM, Andrew Byrd wrote:
>> Another totally different approach would be to simply directly construct
>> isochrones around each relevant street segment and union them all. These
>> would be trapezoids with (approximately) circular endcaps, or the
>> intersection of two such trapezoids coming from opposite ends of the
>> segment. Most street segments could be skipped if they are not within a
>> certain threshold of an isochrone boundary.

This would be an interesting path to follow, probably much more
efficient for large and sparse networks. It has to rely on
ultra-optimized polygon union routines. Also that need to handle
crossing or u-shaped street geometries. JTS has a Geometry::buffer()
operation, I don't know if that work for line string and if that remove
intersecting results.

> On second thought, it would be difficult to skip some segments and make
> assumptions about what's inside or outside an isochrone... but many of
> the street-level access polygons would be simple rectangles.

I do not see what would be difficult? Determining the street geometry
clipping corresponding to the inside is just linear interpolation based
on geometry distance, walk-speed and end-point times?

--Laurent

[Stefan Steiniger]

Hi Laurent,

interesting ideas - and looking fairly well too.
I actually may say that my first implementation was based on someone
else' ideas (is now at Conveyal ;). And my method worked, more or less,
for me as it was developed having a pedestrian in mind, a less so
transit and biking. That is "detail" was important, and also it is kind
of assumed that everything is reachable when leaving a street.
Which actually brings me to the hole problem. I am surprised that in
your model/example, the river is not covered. I would have thought that
you can exclude such areas, like lakes or so, only by introducing
additional information. Because how would you know that you can not walk
towards a direction? I.e. it would mean to analyze additional land-use
info in OSM files, and perhaps, afterwards cookie-cut it from the final
access polygon.

And a thought on interpolation... well, if you want something more
advanced and smoother than what you get with a z-plane, you would
certainly end up with a kernel density function (often used for what is
called heat-maps). Then you could try finding a min-value to draw the
isoline (in animal home range analysis/estimation they often take the
(arbitrary) z-value that would encompass p=95% Volume : of the heat
mountains). However kernel functions require your to define at least the
interpolation-window size (bandwidth h) and then the user is left to
figure what is best ;) (I assume a kernel function calculation is also a
bit more costly?)

and then... see below:

Am 08.10.13 13:00, schrieb Laurent GRÉGOIRE:

>
> This would be an interesting path to follow, probably much more
> efficient for large and sparse networks. It has to rely on
> ultra-optimized polygon union routines. Also that need to handle
> crossing or u-shaped street geometries. JTS has a Geometry::buffer()
> operation, I don't know if that work for line string and if that remove
> intersecting results.

well, there is still the (JTS) danger that an union of buffers leads to
topology exceptions (due to nearly parallel lines that need to be
intersected, e.g. parallel roads that are buffered and the buffers may
barely touch).

Anyway, keep on the good work!

cheers,
stefan

[Laurent GRÉGOIRE]

Hi Stefan,

On 08/10/2013 22:08, Stefan Steiniger wrote:
> I actually may say that my first implementation was based on someone
> else' ideas (is now at Conveyal ;). And my method worked, more or less,
> for me as it was developed having a pedestrian in mind, a less so
> transit and biking. That is "detail" was important, and also it is kind
> of assumed that everything is reachable when leaving a street.

That make sense. In walk shed computation (or car) you probably can't
have islands (at least until teleportation is discovered), the problem
arise with transit where islands are numerous. Also for walk/car
isochrones holes are much less of a nuisance, usually the metric will be
smoother. Since OTP is mainly transit focused having to handle both
cases could be helpful.

> Which actually brings me to the hole problem. I am surprised that in
> your model/example, the river is not covered. I would have thought that
> you can exclude such areas, like lakes or so, only by introducing
> additional information. Because how would you know that you can not walk
> towards a direction? I.e. it would mean to analyze additional land-use
> info in OSM files, and perhaps, afterwards cookie-cut it from the final
> access polygon.

The trick is that the z function properly handle this by itself because
it will find the shortest edge geometry to the considered point, and if
no geometry is found within a specified distance it will bail-out: in
this case we consider the z value to be infinite. This brings the
problem of interpolating with infinite values, but the new algorithm
(soon another mail) solve this by adding the off-road distance to the mix.

> And a thought on interpolation... well, if you want something more
> advanced and smoother than what you get with a z-plane, you would
> certainly end up with a kernel density function (often used for what is
> called heat-maps). Then you could try finding a min-value to draw the
> isoline (in animal home range analysis/estimation they often take the
> (arbitrary) z-value that would encompass p=95% Volume : of the heat
> mountains). However kernel functions require your to define at least the
> interpolation-window size (bandwidth h) and then the user is left to
> figure what is best ;) (I assume a kernel function calculation is also a
> bit more costly?)

More info soon, but indeed the new algorithm make this kind of kernel
density computation, in a discrete fashion, using the off-road distance.
It's hard to say what exactly would be the corresponding continuous
kernel density function as everything is discrete, this is left as an
exercice to the reader :)

> well, there is still the (JTS) danger that an union of buffers leads to
> topology exceptions (due to nearly parallel lines that need to be
> intersected, e.g. parallel roads that are buffered and the buffers may
> barely touch).

A prototype would be nice to make, but I'm quite afraid that the
computation cost would be way too expensive. In a medium-sized city you
easily end-up with 100k edges (lots with simple geometry, but some with
more complex one such as round-abouts modeled with 100 points...) The
buffering only by itself could take too much time, and I'm not talking
about making the union of 100k polygons, with lots of them
intersecting... I would be curious about the performances though.

> Anyway, keep on the good work!

Thanks! I hope all this can be of some help to the community.

--Laurent