[streets] Streets and intersections, first experiments

[polarkernel]

In order to create realistic streets and intersections, many of their features have to be invented by the software, as they are not provided by the OSM data. In the following posts, first problems and findings, experienced with experimental code, are discussed.

[polarkernel]

Collisions with objects

Using default widths for the roads, it can happen that they collide with other objects (grassland, buildings, etc.). For the Karl-Marx-Allee for instance, using a lane width of 3.5 m, we get collisions with the grassland of the median strip, as depicted at the arrows in the left image below (blue: way-sections, green: grassland). In reality, the three lanes have a width of 8.15 m (measured using Google Earth), so that one lane has a width of 2.72 m. Setting the lane width to 2.7 m, the collisions remain, because the center lines drawn by the operator is too near to the median strip (right image).

Some other examples of collisions are depicted below:

Two main problems will have to be solved:

• How to detect collisions efficiently? It is difficult to find them manually.
• How to resolve them consistently?

I tried also to extract some statistics about the street widths from the OSM width tag, but almost none of the ways have this tag. We will have to create our own rules from comparisons with satellite data, I think

[polarkernel]

Simple Intersections

For simple intersections, it is relatively easy to construct the intersection area. The centerlines of the outgoing way-sections are expanded to the left (green in left image) and to the right (red) by half of the way width to get the borders of the carriage way. This can be done using pyGEOS. Then the intersections between a left border with the right border of the following way-section in counter-clockwise order is computed (dots in left image). There, circular arcs are fitted to the angles built by the way borders (middle image). The distance to the farthest endpoint on the left or right way-border of every outgoing way is then used to construct the perpendicular connector line of the way. This finalizes the polygon of the intersection area (right image).

In many cases, this leads to pleasing results (different scales have been used for these images, the radius of the circular arc is constant for these examples):

There are in detail some issues to solve. Sometimes, the circular arc can't be fit to the intersecting way-segments of the way-section, because they are too short so that the tangent points of the arc are outside the segment. (left image below). A first solution, trying to intersect more distant way-segments and fit the circular arc to them, failed in many cases, because the way-section were too short. In the current solution, the radius of the arc is decreased iteratively until it fits to the original segments. But then the radius of the arc depends on how long the intersecting segments have been drawn. The radii are then different for the same intersection (right image) and can become very small. So far, a better idea is missing.

Another issue arises, when there is no intersection between left border and right border (on top of the left image below). Then just a straight line may be used, as shown in the right image, or an S-shaped curve has to be constructed.

We will also have to find a rule for the radius of the circular arcs. One formula I have found is derived from the turning radius estimation based on a centrifugal force formula. The radius r becomes there r = V_max / (µ_s g), where V_max is given by the max speed tag in OSM, g is the gravitational constant and µ_s the coefficient of friction between the wheel and the surface of the road (~0.7).

[polarkernel]

Further Problems

Street width cannot be determined solely by the number of lanes. In the example below, the Karl-Marx-Allee widens from three to four lanes. In the satellite image on the right, however, you can see that the additional lane is a turning lane, so the road only widens on one side of the OSM centerline. We will have to introduce some concept of lanes and their orientation relative to the OSM centerline. This will also be useful for the realistic representation of the roads.

We will need some intelligent algorithms to construct realistic way-ends:

Finally, the problem of complex clustered intersections is still completely unresolved:

Converting OSM tags to relevant information is complicated. Perhaps we can share some knowledge with A/B Street, where the development team is still highly active currently (see osm2lanes on GitHub).

[vvoovv]

Sometimes, the circular arc can't be fit to the intersecting way-segments of the way-section, because they are too short so that the tangent points of the arc are outside the segment.

In which OSM file can I find this case? An where exactly is it located?

[polarkernel]

In which OSM file can I find this case

The depicted case is from _berlin_karl_marxallee.osm, at the left top intersection of the roundabout. The intersecting OSM way ids are 932638437 and 147646845.

I committed some code where I have highlighted all the locations in green, where currently the radius is reduced due to short intersecting segments. Use the arguments --highways --roadIntersections --buildings to run the code. Note that this is just experimental code producing a lot of crap, just useful for me currently.

[polarkernel]

Unexpectedly, a new problem appeared. The code of the function _paralleloffset() in pyGEOS sometimes returns multiline output. The multiline output appears to be missing a large gap of the expected offset curve. The bug is known since 5 years and also occurs for example in shapely or other libraries that are based on GEOS. It is still unresolved. I tried unsuccessfully to trace the code line by line to find out where the bug occurs.

Now I decided to write my own code for both, this function and, related to it, the expansion of a line to a polygon (the function buffer()). Since roads never have self-intersections, a lot of code can be omitted. As I have found, pyGEOS also executes an incredible number of unnecessary function calls, which take a lot of execution time. I am convinced that own code will be faster.

[vvoovv]

I am convinced that own code will be faster.

Sounds good.

[vvoovv]

Just an idea.

I suggest so start from an approximate solution of the problem of an area partition in accordance with OSM data.

We know all graph-cycles. We know which objects are located inside each graph-cycle.

From all graph-cycle we select "narrow" ones. How to find the "narrow" graph-cycles is to be developed.

We offset the way-segments for a given graph-cycle. If all objects located within the graph-cycle are OSM landuses (for example, grass). we assume that they occupy the remaining space in the graph-cycle.

I think this approach can solve the configurations like Karl-Marx-Allee and its sidewalks and cycleways.

[vvoovv]

The depicted case is from _berlin_karl_marxallee.osm, at the left top intersection of the roundabout. The intersecting OSM way ids are 932638437 and 147646845.

This intersection has sharp angles in the reality:

[vvoovv]

The radii are then different for the same intersection (right image) and can become very small. So far, a better idea is missing.

Can't we use the minimum radius of all corners of the intersection?

[polarkernel]

Can't we use the minimum radius of all corners of the intersection?

This could possibly be a viable idea. Eventually, the angle could also play a role, as for the example of the flat intersection of the Karl-Marx-Allee above.

However, I don't like the fact that these radii shall depend on the operator who drew the roads. Maybe some restrictions can be applied there. Possibly something like a straight angle removal or a spline fit. There exist also regulations (at least in the EU and the USA) on how intersections should be built. Eventually an idea can be found from there.

I am also currently looking for an idea on how to easily overlay a drawing of a result onto a satellite image, so that experiments can be verified in a simple way. Using Photoshop, this is very tedious.

[vvoovv]

I am also currently looking for an idea on how to easily overlay a drawing of a result onto a satellite image, so that experiments can be verified in a simple way. Using Photoshop, this is very tedious.

Do you think using Blender for that would be convenient? I can develop a renderer for that.

[vvoovv]

It seems to be possible to draw over an image using matplotlib. https://stackoverflow.com/questions/34458251/plot-over-an-image-background-in-python https://stackoverflow.com/questions/15160123/adding-a-background-image-to-a-plot

Downloading and stitching the satellite imagery together using the addon in the command line mode requires setting:

bpy.context.scene.blosm.commandLineMode = True

[polarkernel]

It seems to be possible to draw over an image using matplotlib.

This would be very convenient! Is there also a method available to create the georeference between pixels and the scene data? Where has this setting to be added?

[vvoovv]

I will write the code to load an image stitched from tiles in matplotlib. This will be done in the coming days.

I will probably need your help for the matplotlib integration.

The target branch is streets_intersections. Right?

[polarkernel]

I will write the code to load an image stitched from tiles in matplotlib.

Thanks!

The target branch is streets_intersections. Right?

Yes.

[vvoovv]

I will write the code to load an image stitched from tiles in matplotlib. This will be done in the coming days.

Created #42 to track the progress.

[polarkernel]

I am convinced that own code will be faster.

The own code for the creation of offset lines and buffered lines is written and committed. I used a simplified version of the algorithm described in this paper. In contrast to the general algorithm used there and also in pyGEOS, I introduced some additional preconditions to simplify the solution:

• The way-section must not intersect itself.
• The (offset) right border of a way-section must not intersect its left border.

The first one is satisfied by the removal of way-intersections before the way-network gets built. So far, I have never found a violation of the second condition, nor can I imagine how this situation could arise. However, I will keep an eye on it.

Input and output vertices use the mathutils.Vector class, so that no type conversion to Coordinate is required. There was a big surprise when I compared the new solution with the method of pyGEOS. For example, when I tested them with the following way-section from _beijingcenter.osm as input (black: way-section, red: offset line)

and using the following code

distance = 10
testline = [Vector((1069.54638671875, -533.5795288085938)),  ...... ]
polyline = PolyLine(testline)
offsetGen = WayOffsetGenerator()

geosF = GeometryFactory()
coords = [Coordinate(v[0],v[1]) for v in testline]
geosLine = geosF.createLineString(coords)

print('BLOSM',timeit.timeit('offsetGen.parallel_offset(polyline,distance)', number=100, globals=globals()) )
print('GEOS ',timeit.timeit('geosLine.parallel_offset(distance)', number=100, globals=globals()) )

I got the result

BLOSM 0.1132895
GEOS  6.2422631

I never expected such an improvement of speed!

[vvoovv]

I got the result

BLOSM 0.1132895
GEOS  6.2422631

I never expected such an improvement of speed!

What an incredible performance boost!

[vvoovv]

I tried to run the script with the arguments

--highways --roadIntersections

I am getting the error:

ModuleNotFoundError: No module named 'action.plotUtilities'

[polarkernel]

ModuleNotFoundError: No module named 'action.plotUtilities'

Sorry, plotUtilities is a module that I use since a while only for debugging purposes. It will not stay in the branch. For the moment I committed it together with a new special case of way-sections. These are a way-sections that have both ends at the same intersection, topologically in a network, they are called loops, a case that was not yet treated correctly before.

Now, the script using the arguments --highways --roadIntersections should run correctly.

[polarkernel]

I committed a new version, where WaySection, Intersection and PolyLine have been completely refactored. There were a lot of bugs before, some of them very subtle. The module plotUtilities is now removed. Most of the intersection areas are now computed correctly, as long as they do not overlap. The street widths and the radii of the intersection fillets are still not yet detailed. Another innovation is the expansion (buffering) of the way center lines to polygons. An own method has been written, with similar gain in speed as already reported. Until now, pyGEOS is no more used in these modules.

As a next step, I will write a lot of code comments, so that I will still remember my thoughts in two months. Then I will care on street widths and the radii of the intersection fillets. The satellite images will be very helpful for this.

[vvoovv]

How can tram tracks (railway=tram) be added?

[polarkernel]

How can tram tracks (railway=tram) be added?

I did not yet care about them. But aren't they very similar to streets? I would assume that they could be feed into the same processing pipeline, just with their own rules for widths and radii at crossing areas. But this still needs to be investigated.

[vvoovv]

Can we define the types of ways that took part in the calculation of intersections in defs.way?

For example,

roadwayIntersectionCategories = (
    ...
)

railwayIntersectionCategories = (
    "tram",
    "light_rail",
    "funicular"
)

wayIntersectionCategories =  roadwayIntersectionCategories + railwayIntersectionCategories

[polarkernel]

Can we define the types of ways that took part in the calculation of intersections in defs.way?

Yes. I will additionally also introduce some assignments of default widths and intersection radii for them. I will put these also to defs.way. In a first (short and fast) search, I found that mainly these widths seem to be different in US than in Europe. Could you imagine that we introduce something like internationalization and localization? Eventually by PML?

To introduce railwayIntersectionCategories, will we additionally need the class RailwayManager or can railways be provided by the WayManager?

[vvoovv]

Could you imagine that we introduce something like internationalization and localization? Eventually by PML?

Yes, we will need to introduce a kind of "locale" in PML later.

To introduce railwayIntersectionCategories, will we additionally need the class RailwayManager or can railways be provided by the WayManager?

Both railways and pavement ways are treated by WayManager in the dev branch: line 60 line 142

[polarkernel]

Both railways and pavement ways are treated by WayManager in the dev branch:

I don't get any member of allRailwayCategories by wayManager.getAllWays(), neither in the branch _streetsclusters nor in _streetsintersections. The file setup/init.py, where you point on lines 60 and 142 in dev, is empty in these branches. I am sorry, but I am not able to reconstruct the changes that you did between these branches. Can you help?

[vvoovv]

I'll care about this problem today.

[vvoovv]

My bad! It works in its current state. It is simply required to add the command line argument --railways:

[polarkernel]

My bad! It works in its current state. It is simply required to add the command line argument --railways:

I hope we can solve many more problems this way :)

[vvoovv]

Just an idea.

We could try to cluster adjacent narrow graph-cycles. Geometry can be created for clusters that form carriage ways. Then geometry can be created for the clusters formed by cycleways and footways and adjacent to the carriage ways clusters.

A robust method of detecting a narrow graph-cycle will be required. One could you use the ratio of the squared half-perimeter to the area of the graph-cycle as a measure of its narrowness: factor = half_perimeter * half_perimeter / area The factor should be large enough for a graph-cycle to be narrow.

[polarkernel]

In the recent commit, I introduced more way properties. The following changes were made:

• In defs/way.py, the tuples wayIntersectionCategories, roadwayIntersectionCategories, and railwayIntersectionCategories have been introduced.
• In way/manager.py, a new method getAllIntersectionWays() has been introduced, which delivers all ways in wayIntersectionCategories.
• In a new file _way/wayproperties.py, values and methods have been introduced that describe street widths depending on their category and lane information. Additionally, radii for fillet arcs have been defined for different way categories. These values have been evaluated using several of our data files. However, they have to been improved by increasing experience.
• In _way/wayintersections.py, The radii of fillet arcs are reduced continuously when the angle between the roads gets smaller.

[vvoovv]

Just an idea. We could try to cluster adjacent narrow graph-cycles.

Do you think this idea makes sense?

[polarkernel]

Do you think this idea makes sense?

I did not yet think a lot about, at a first glance, yes. The whole task is quite complicated and needs to consider a lot of details. I first try to find the intersections and the lanes somewhat properly. This gives me a deeper idea of the complexity of the problem. It will also depend on what details need to be distinguished within the clusters (turning lanes, etc., for example). I would like to postpone this issue for now, if you agree.

[vvoovv]

I would like to postpone this issue for now, if you agree.

Ok.

[polarkernel]

In my further experiments, I studied transitions. Transitions are second-order "crossings", i.e. points between two way-sections that are of a different category or have a different number of lanes. These must already be detected when the way-section network is created and get assigned a special node in the network. In the branch _wayintersections, they are replaced by a transition area, which adapts the different way widths, if any. The connecting borders are constructed by Catmull-Rom splines that keep the direction of the border endings of the existing ways and may be curved. These fit the ways quite smooth, as can be seen for the green areas in the following pictures:

To create these areas, the way-sections connected there have to be trimmed a bit and the question is, what length this area shall have. When the width difference between the ways is large, a short length results in an unnatural transition, as shown in the left picture below. Therefore, I made this length dependent of the width difference, given limits between 2 m and 20 m. The result is shown in the right picture:

A special case are way-sections that change their width within the section. This is the case when a turn lane is to be created within the section. The transition is then within the way-section and changes its shape. The detection and distinction of these cases is by far not yet complete, and the inspirations of the operators of OSM will keep us busy for a long time to come. The shape of the turn lane is currently half of a parabola, although I found sketches that propose to use two halves:

Turn lanes may be very different, but there is no data in OSM to distinguish them. For some of them, an additional transition area is required to fit the way ends, as depicted in the middle image below. Sometimes, there are traffic islands, as shown in the rightmost image below, that are not marked in OSM.

As mentioned above there is still a lot to do, but for now I will stop development at this point, so we don't get lost in details before we have an overview of the whole project. I'm now going to spend some time looking at clusters of intersections and streets.

[vvoovv]

A turn lane can be generated or removed using a kind of interactive editor right in the Blender's 3D View.

[polarkernel]

A robust method of detecting a narrow graph-cycle will be required. One could you use the ratio of the squared half-perimeter to the area of the graph-cycle as a measure of its narrowness: factor = half_perimeter * half_perimeter / area

For the following experiments, I used the reciprocal of your suggestion, as it is more suitable for visualization. Then it is known as the Polsby–Popper test. The equation

PP = 4π A/P²

where A is the area and P is the perimeter of the graph-cycle, describes its compactness in the range [0…1]. For a circle, PP is 1 and for an elongated structure, it is smaller, which is what we are looking for. For the visualization, I used the term 1-PP to get a high value for parallel streets. The result is similar to your proposition, but your equation delivers somehow unbounded values, which makes it difficult to visualize them.

For some of our datasets, the results look very encouraging:

But for others they are completely unusable:

The reason is that complicated structures like the one below also have a large perimeter and a small area, they are not compact. But they are not elongated, as streets are expected to be. Do you know a measure for elongation? I don't even know a correct term for that in English.

[vvoovv]

The reason is that complicated structures like the one below also have a large perimeter and a small area, they are not compact. But they are not elongated, as streets are expected to be. Do you know a measure for elongation? I don't even know a correct term for that in English.

I hope tiling should help for graph-cycles like those.

[polarkernel]

I hope tiling should help for graph-cycles like those.

Unfortunately, tiling does also split the long graph-cycles between the roads of a road cluster. The reason is obvious, long graph-cycles would require a large search radius for the KD-tree, which would include too many objects to subtract.

[vvoovv]

Tiling of an elongated graph-cycle into 2 smaller graph-cycles increases PP for both parts approximately by 2.

I'd expect an increase by a larger factor for at least one tile for a graph cycles like here.

[polarkernel]

Tiling of an elongated graph-cycle into 2 smaller graph-cycles increases PP for both parts approximately by 2. I'd expect an increase by a larger factor for at least one tile for a graph cycles like here.

You are right! I am surprised, but this seems to be a viable approach. Below, you find the results for the same files when tiling is applied. I will try to elaborate this method and already have many ideas. I will also need to create a data structure that allows easily to find the neighbors of a graph-cycle.

[vvoovv]

Below, you find the results for the same files when tiling is applied.

Thank you for the encouraging result

I really enjoy looking at these mosaic pictures! :)

[polarkernel]

I really enjoy looking at these mosaic pictures! :)

Then, here are some more results. I tried, as another indicator of lengthiness, using a minimum bounding rectangle (I supplied the code some time ago in this repo). The ratio height/width (where height < width) of such a rectangle gets smaller, the more eccentric the graph-cycle is. Compared to the Polsby–Popper (PP, left) test, the bounding rectangle (BR, right) generates less false positives (I displayed 1-height/width to make the plots comparable):

As an additional advantage, width and height are physical dimensions and can be limited. The height of the bounding rectangle corresponds to the distance between the roads in a cluster. This distance is limited, for instance to 15 m, as an example, which again decreases the number of false positives (the fine white lines are all the ways, to ease the navigation):

I also tried to apply automatic classification methods I know for the binarization of images, unfortunately without success. The histograms nor for PP neither for BR show classes and vary significantly from scene to scene.

Now I want to try additionally another completely different approach. Clusters of ways always end at clusters of intersections. It should be relatively easy to find clusters of intersections using the results of generating polygons of intersections. Way clusters are then the shortest connections between these intersection clusters in the way-sections network. There exist suitable and efficient graph algorithms for this purpose.

[vvoovv]

Clusters of ways always end at clusters of intersections.

In our example file _streets/berlin_highways_withoutintersections.osm clusters of ways do not end at clusters of intersections in the left and right sides of northern part of the scene.

[polarkernel]

In our example file _streets/berlin_highways_withoutintersections.osm clusters of ways do not end at clusters of intersections in the left and right sides of northern part of the scene.

This is not a problem. To limit cycles at the scene border, I had to introduce the scene border segments as virtual way-sections of the category 'scene_border'. The intersections of the real way-sections with these virtual sections produce useful clusters:

However, playing around with clusters, I encountered a lot of subtle difficulties. Currently, a cluster is defined as a group of intersections, that are connected by way-sections of important categories and not longer than a given distance (15 m in the example above). This distance is critical. For instance, on the left end, the motorway lines and the railroad tracks are separated, but not on the right end. Or the blue cluster on the top right ends to the left in two clusters (green and orange) although this should be a sngle way-cluster.

I am not yet sure whether these clusters are a good idea. Maybe I will go back to the previous method using lengthy cycles.

[vvoovv]

The approach with lengthy cycles seems to give more reliable results than the one with the clusters of intersections.

[polarkernel]

The approach with lengthy cycles seems to give more reliable results than the one with the clusters of intersections.

I agree, currently I stopped the investigations with the approach that uses intersection clusters.

Currently, I am investigating whether tiles with smaller values can be merged to the cluster if they are neighbors of tiles with larger values. But due to reasons that are difficult to explain in a few words, the sliced graph-cycles are each independent polygons, that often have no common vertices with their neighbors. Let me illustrate that with an example. The polygons 1 and 2 in the left plot seem to have common vertices, pointed to by the red arrows. But when I shrink the tiles a bit to make them visible individually, one can see that this is not the case (right image). They have no common vertices.

I am now seeking for a method to find the neighbors of a polygon. Naturally, these are not always quadrangles. ArcGIS provides a Polygon Neighbors (Analysis) tool, but there is nowhere a description on how it works. I tried for instance to use the centers of the polygons and to construct a Delaunay triangulation to get the neighbors. The blue lines in the image below show the result for the same section as above. But this does seemingly not work.

Maybe I could use the directions given by the minimum bounding rectangles to point to neighbors.

Should you have any idea on how to find the neighbors of a polygon, please tell me.

[vvoovv]

I'd like to mention the approach with the classes Vector and Edge that we discussed in #40 starting from here. I use this approach to find the neighbors of a building or a building part.