Giter Club home page Giter Club logo

pyvisgraph's Introduction

Pyvisgraph - Python Visibility Graph

MIT License PyPI

Given a set of simple obstacle polygons, build a visibility graph and find the shortest path between two points.

Figure 1

Pyvisgraph is a MIT-licensed Python package for building visibility graphs from a list of simple obstacle polygons. The visibility graph algorithm (D.T. Lee) runs in O(n^2 log n) time. The shortest path is found using Djikstra's algorithm.

To see how visibility graphs work interactively, take a look at the Visibility Graph Simulator built with Pyvisgraph.

Installing Pyvisgraph

$ pip install pyvisgraph

Pyvisgraph supports Python 2 and 3.

Usage

Here is an example of building a visibility graph given a list of simple polygons:

>>> import pyvisgraph as vg
>>> polys = [[vg.Point(0.0,1.0), vg.Point(3.0,1.0), vg.Point(1.5,4.0)],
>>>          [vg.Point(4.0,4.0), vg.Point(7.0,4.0), vg.Point(5.5,8.0)]]
>>> g = vg.VisGraph()
>>> g.build(polys)
>>> shortest = g.shortest_path(vg.Point(1.5,0.0), vg.Point(4.0, 6.0))
>>> print shortest
[Point(1.50, 0.00), Point(3.00, 1.00), Point(4.00, 6.00)]

Once the visibility graph is built, it can be saved and subsequently loaded. This is useful for large graphs where build time is long. pickle is used for saving and loading.

>>> g.save('graph.pk1')
>>> g2 = VisGraph()
>>> g2.load('graph.pk1')

For obstacles with a large number of points, Pyvisgraph can take advantage of processors with multiple cores using the multiprocessing module. Simply add the number of workers (processes) to the build method:

>>> g.build(polys, workers=4)

Pyvisgraph also has some useful helper functions:

  • g.update([list of Points]): Updates the visibility graph by checking visibility of each Point in the list.
  • g.point_in_polygon(Point): Check if Point is in the interior of any of the obstacle polygons. Returns the polygon_id of said polygon, -1 if not inside any polygon.
  • g.closest_point(Point, polygon_id): Return the closest point outside polygon with polygon_id from Point.

For further examples, please look at the scripts provided in the examples folder.

Example & performance

This example uses a shapefile representing world shorelines as obstacles. Two vessels were picked randomly and their current location found using AIS. Red lines are laden voyage legs (carrying cargo) and dotted blue lines are ballast legs (no cargo, moving to load destination). Pyvisgraph has the following performance on a Microsoft Surface Pro 3 (Intel i7-4650U @ 1.7Ghz, 8GB DDR3 RAM), where time is in seconds:

Shoreline obstacle graph [points: 4335 edges: 4335]
Using 4 worker processes...
Time to create visibility graph: 554.683238029
Visibility graph edges: 118532
Time to update visgraph & find shortest path: 1.09287905693
Shorest path nodes: 19
Time to find shortest path between existing points: 0.508340835571

For one vessel the origin and destination Points were not part of the built visibility graph and had to first be computed. For the second vessel, these Points were first added to the visibility graph using update, then finding the shortest path is faster. Using Matplotlib basemap to visualize the routes: Figure 2

For more information about the implementation, see these series of articles:

pyvisgraph's People

Contributors

deboc avatar jianantian avatar taipanrex avatar tjansson60 avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

pyvisgraph's Issues

Detecting points on the edges of a polygon

Would it be possible to alter the visibile_vertices function to detect points on the edge of a polygon as well?

For example in the following example:

p1 = vg.Point(0.20,0.1)
p2 = vg.Point(0.50,0.5)
p3 = vg.Point(0.49,0.5)

polys = [[(0,0),(1,0),(0,1)]]
G = vg.Graph(polys)

vg.visible_vertices.visible_vertices(p1,G1,destination= p2) will return [Point(0.00, 1.00), Point(0.00, 0.00), Point(1.00, 0.00)]
vg.visible_vertices.visible_vertices(p1,G1,destination= p3) will return [Point(0.49, 0.50), Point(0.00, 1.00), Point(0.00, 0.00), Point(1.00, 0.00)]

Would it be possible to maybe alter the visibile_vertices function to get [Point(0.00, 1.00), Point(0.00, 0.00), Point(1.00, 0.00), Point(0.50, 0.50)] from vg.visible_vertices.visible_vertices(p1,G1,destination= p2) as well? Or is there another way to achieve this result?

open_edges: consider when initializing checking against bounding box first

When open_edges is initialized, the scanline checks intersection against all Edges in graph. It might be faster on graphs with many obstacle polygons to check if the scanline crosses the bounding box of each polygon first. For the boxes it does not cross, all edges in that polygon can be ignored.

Extending to 3D possible?

I am using this to model a city with the buildings and no-fly zones(NFZ) which will have 3d representations, for navigating an UAV through these buildings and NFZs. Has there been any attempt to extend this to 3D? Or are there any leads in this direction? If not I would love to make an attempt at extending this to 3D

Implement AVL Tree for open_edges

An AVL Tree has O(log n) insert, delete and search, which is better than Python lists' sort O(n log n) and delete O(n). Critical for the performance of Lee's algorithm.

Rounding error in ccw function causes mistakes in edge intersection evaluation

In the ccw function defined in visible_vertices.py, the term area = (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x) is evaluated. Rounding error caused by Python number representation when computing the term can mean that the function evaluates colinear points as noncolinear. This can cause lots of problems later on (evaluated many edges incorrectly in my particular application's visibility graph).

I am using Python 2.7 so I cannot make use of math.isclose(), so I defined an isclose method in the file and then used that to check if the area as computed by ccw was close to zero with an absolute tolerance of 0.0000000000000001.

Progressbar output not functional for workers>1

As I have mentioned before I really enjoy your project and I have an idea for an improvement.

One worker
When building a graph from polygons derived from shorelines such as GSHHS_c_L1 containing 742 shapes using 1 worker the output works nicely even though there is no indication in time when it will be done:

[Done][Rem.][Avg t]                                                                                                     
[  79][1608][0.051]                                                                                                     

This is sufficient for small graphs, but complex sets of polygons such as the intermediate shoreline set as GSHHS_i_L1 this can take a long time with one worker and with a time estimate it is even harder to plan around and perhaps do over the weekend.

More than one worker
Using more than one worker the different workers overwrite the progress output on the same line and the output is garbled and it is even harder to find out when the process will be done. Using say 16 workers it is almost impossible to decipher the progress.

 [Done][Rem.][Avg t] [Done][Rem.][Avg t] [Done][Rem.][Avg t] [Done][Rem.][Avg t] 
[63.0C[ 281][ 140][0.059]

Proposed solution
I was thinking of rewriting the current progress bar with a modern alternative such as tqdm or progressbar2, but wanted to make sure that this is inline with your thoughts on the progress bar before I started.

License file is not packaged in the PyPi tarball

It is generally best to include the license file in the distributed source tarball. It is helpful for downstream packaging (and generally expected). Can you add that to the next release on PyPi?

shortest_path when two solutions, issue with Python 3 and heapq

If there are two points of equal distance, the priority in the priority_dict will be same. heapq will then revert to the next value (which are the Points themselves) and try and compare. Point does not have a lt (or a cmp for that matter, so not sure why this isnt a problem in Python2) and it raises a ValueError as Points are not comparable.

A quick fix is to add a lt to Point class, and just compare their hash. This way its kind of random which solution is chosen.

Performance hit & js port

Hi there @TaipanRex

I just wanted to say thanks for this awesome library and the associated blog writeups! I've been porting it across to js here

In doing so I was tinkering with the performance and found that the int() operation this line made a significant difference to my performance.
On data with 1400 points calculating the graph takes either 4 seconds (using the existing float) vs 30 seconds (converting on the integer). I don't know if it's just the js parseInt function that is super slow but just thought I'd give you a heads up.

Anyway thanks again :)
Rowan

UnboundLocalError: __lt__() in vis_graph.py

same_point is not defined and should be. Initiate to None, then throw an exception if still None after the two if statements. There are instances where a point of a polygon is on the edge of another polygon. That means their distances are equal, but have no common points. same_point therefore does not exist...

Outdated link to GSHHS shapefile for 1_build_graph_from_shapefiles.py: add this update and clarification

In the examples folder and the 1_build_graph_from_shapefiles.py file, change the description from line 27:31 with the following:

# In this example a visibility graph will be created from the GSHHS shorelines
# available at: http://www.soest.hawaii.edu/pwessel/gshhg/
# 1. Download the shorelines in the format of shapefiles and select the wanted resolution 
#    and level (in this example we use GSHHS (type), Crude (quality), L1 (region)). 
# 2. Copy all the related files (shx, shp, prj, dbf) from the following subfolder (GSHHS -> c) 
#    and only the *C1 relevant files to the folder that contains this script.

Support of polygons with interior holes

Request for support of polygons with interior holes. I think the core of the code will work as is if there is some way to represent them.

For now my hack is this: when creating the list of polygons for the graph creation, specify the exterior polygon ring first, then put a blank polygon (empty list) in front of every matching interior ring. Then in the init method for the graph I catch those empty lists and rewind the pid by one. The resulting behavior for the visibility graph seems fine so far.

Add Pyvisgraph documentation

I have started developing this on the documentation branch. I will host these on readthedocs.io, using Sphinx.

There will be two parts to the documentation 1) from the user perspective 2) from the developer perspective.

Start or end point on polygon boundary

When the start or end point is on the boundary of a polygon, VisGraph.shortest_path() throws an error, as demonstrated in this test case:

import pyvisgraph as vg
polys = [[vg.Point(3.0,2.0), vg.Point(6.0,2.0), vg.Point(6.0,7.0), vg.Point(3.0,7.0)]]
g = vg.VisGraph()
g.build(polys)
s = vg.Point(3.0,5.0)
e = vg.Point(8.0,5.0)
path = g.shortest_path(s, e)
print path
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-64-23f53555023b> in <module>()
      5 s  = vg.Point(3.0,5.0)
      6 e = vg.Point(8.0,5.0)
----> 7 path = g.shortest_path(s, e)
      8 print path

pyvisgraph\vis_graph.py in shortest_path(self, origin, destination)
    128             for v in visible_vertices(destination, self.graph, origin=orgn):
    129                 add_to_visg.add_edge(Edge(destination, v))
--> 130         return shortest_path(self.visgraph, origin, destination, add_to_visg)
    131 
    132     def point_in_polygon(self, point):

pyvisgraph\shortest_path.py in shortest_path(graph, origin, destination, add_to_visgraph)
     68         path.append(destination)
     69         if destination == origin: break
---> 70         destination = P[destination]
     71     path.reverse()
     72     return path

KeyError: Point(8.00, 5.00)

Not sure if this is the intended behavior.

It seems unusual that pyvisgraph will run the shortest path along a polygon boundary, e.g. as in the illustration below, yet not allow the start or end point to be on the boundary.

image
(source code for this illustration is here)

In this application of a pixelated map with integer coordinates, it's not a trivial issue: the start and end points must be 1-pixel outside the polygon, yet the interior points of the path can be on the polygon boundary.

My clunky workaround is to use Shapely to move the start and/or end points outward along the edge normal by 0.1 units. This avoids the error, and gives the correct result when the path coordinates are rounded back to integers.

Bug

Help please.
I didn't use shape to specify the obstacles. I just manually input the polygons like:
polys = [[vg.Point(0.0,1.0), vg.Point(3.0,1.0), vg.Point(1.5,4.0)],[vg.Point(4.0,4.0), vg.Point(7.0,4.0), vg.Point(5.5,8.0)]]
Is it wrong? The boundary is specified by the convex hull vertices of all points in a polygon.
2

Requesting support for smaller distances

I found that this doesnt work for very small shapes. Because in the result, coordinates are rounded to 2 decimal places.

So i kindly request you to add a parameter for decimal places which will support feet or meter accuracy.

Crash when vertices overlay

Let's say I have a first shape with points [(0, 0), (1, 0), (1, 1)] and a second shape with points [(0, 0), (1, 0), (1, -1)]
This crashes the program with error :

File "pyvisgraph/vis_graph.py", line 82, in build
    for edge in _vis_graph(self.graph, batch):
  File "pyvisgraph/vis_graph.py", line 156, in _vis_graph
    for p2 in visible_vertices(p1, graph, scan='half'):
  File "pyvisgraph/visible_vertices.py", line 110, in visible_vertices
    insort(open_edges, k)
  File "pyvisgraph/visible_vertices.py", line 346, in insort
    if x < a[mid]: hi = mid
  File "pyvisgraph/visible_vertices.py", line 388, in __lt__
    aslf = angle2(self.p1, self.p2, self.edge.get_adjacent(same_point))
UnboundLocalError: local variable 'same_point' referenced before assignment

Exceptions: polygons

If try to add a incorrect polygon, should raise an exception. Should be simple polygons, i.e. min. 3 points, no intersecting lines etc.

Visibility graphs have a huge number of polygon crossings

I am trying to create a visibility graph from a shapefile, following the examples.
The example programs run fine, and I am able to get a correct path.
For my own shapefile, the visibility graph seem to never avoid crossing the polygons. So solution paths always cross a polygon if a polygon exists between them.

Here is my original map, plotting the polygons with matplotlitb.
sample_region

For faster debugging, I made a reduced shapefile with only two polygons.
Here, I plotted all the edges in graph.visgraph.get_edges(). You can also see the solution path passing through the polygon.

sample_region_path

Any idea what might cause this? Happy to share code if it helps.
The code starts by converting a raster to polygons with fiona, saving them as a shapefile, then following example 1. Not too experienced with geospatial data, so I could have made an error here. But the polygons load fine in QGIS.

Handling voyages crossing the antimeridian

First of all, thanks for the very nice library. I really enjoy it. I think I have found a bug, which I have detailed below.

I am using the polygon shoreline set GSHHS_c_L1 and I am trying to calculate the distance between the the two points:

p1 = (179, 0)
p2 = (-179, 0)

The true distance between the points is probably something close to the geopy.distance.greater_circle distance of 120 nm as the points are very close to each other with no polygons in the way.

The pyvisgraph distance is 21595 nm and the calculated graph is shown below using folium. It is quite clear that the problem is that pyvisgraph does not take into account the crossing of the antimeridian. Is there any way this could be enabled or is there something I have overlooked?

antimederian

VG generate wrongly due to polygon_crossing function not accurate

  below is a convex polygon, each point should only have 2 visible points, but point_3 got 3.
  polygon_crossing function need review, as it return false when checking mid point (208.5, 227)

    self.point_1 = Point(298, 163.0)
    self.point_2 = Point(263, 168.0)
    self.point_3 = Point(199, 194.0)
    self.point_4 = Point(192, 227.0)
    self.point_5 = Point(218, 260.0)
    
    g = vg.VisGraph()
    g.build([[self.point_1, self.point_2, self.point_3, self.point_4, self.point_5]])

Killing Pyvisgraph with CTRL-C when running multiple processes does not work

If you are running build() with multiple workers, CTRL-C will not kill them and you have to end the processes through the task manager/kill in separate terminal. The below solution should work:

def do_work(i):
    try:
        print 'Work Started: %d %d' % (os.getpid(), i)
        time.sleep(2)
        return 'Success'
    except KeyboardInterrupt, e:
        pass

def main():
    pool = multiprocessing.Pool(3)
    p = pool.map_async(do_work, range(6))
    try:
        results = p.get(0xFFFF)
    except KeyboardInterrupt:
        print 'parent received control-c'
        return

IndexError

Hi! This is an amazing library and I have really enjoyed working with it. I have come across an issue I haven't been able to root cause. Some of the shapefiles I have used (such as: https://geodata.lib.berkeley.edu/catalog/ark28722-s7d02x) will have an IndexError during if self._open_edges[index] == edge, where the index is -1 (calculated from index = self._index(p1, p2, edge) - 1 in the delete function of class OpenEdges. I have been unable to figure out why some shapefiles would lead to this. I have examined and cleaned the polygons closely in QGIS. Any ideas on how to overcome this?

VisGraph.shortest_path() updates self.visgraph

See test test_shortest_path. I had first made a new visible_edges list for origin and destination, then make a shallow copy of visgraph and combine them. The problem is that visgraph is a dict of lists and a Python shallow copy will only go 1 level deep, so when adding edges to the copy dict, its adding to visgraph as well. I am considering reverting the new edges by tracking a list of new edges and using set operations to remove them.

shortest path crossing the obstacle

_Hi
I was trying to use the algorithm with shapely polygons.
To pass the points to the graph I'm using polygon.exterior
When I run the algorithm I have the shortest path that is crossing the polygon. Below are the link to the plot and zoom of the area where the path crosses the polygon.
All the land is represented by one polygon, no multipolygons or overlaping polygons.

https://www.dropbox.com/s/rlyxf8e7sw6yn0t/land%20crossing.png?dl=0 (polygon plot with shortest path)
https://www.dropbox.com/s/i6b5ziwplauxqj9/land%20crossing%20zoom.png?dl=0 (zoom)

Any clue how to solve this problem and what am I doing wrong ?_

SOLVED:
I checked than when you call polygon.exterior you might have points with a huge number of decimal points, so just round the points coordinates to a limited number of decimal points for the coordinates. This solved my problem

Optimize removing edges that are inside polygons

Currently implemented by checking if the mid-point of each visible edge is on the interior of a polygon using point_in_polygon(). This runs in O(n^2), is there a faster way?

One optimization is to create a polygon class and track if it is convex or concave. Point_in_polygon is only needed for concave polygons.

Visibility graph bug

I tried to build a visibility graph with the following input:
[[vg.Point(x=119.55624999999999, y=16.48333333333367), vg.Point(x=203.7872976182326, y=81.64246915103985), vg.Point(x=203.78733213720346, y=112.326984303335), vg.Point(x=189.70101758191677, y=112.326984303335), vg.Point(x=184.99113388465582, y=75.10000000000025), vg.Point(x=184.99113388465582, y=103.01666666666668), vg.Point(x=133.69454408285583, y=103.01666666666668), vg.Point(x=133.69454408285583, y=75.10000000000025), vg.Point(x=131.22787741618916, y=75.10000000000025), vg.Point(x=131.22787741618916, y=103.0166666666666), vg.Point(x=119.55624999999999, y=103.01666666666759)]].

This was the result:
image

I solved this by rounding the points before giving them as input.
I suggest fixing this bug.

Suggestion: examples folder

I was wondering if might be a good idea to have an examples folder in the library showing how to use the package.

For myself at least I needed write some extra code to use the package to read shapefiles, calculate shortest distance and plot the voyages on maps, so I was thinking that a small examples folder could contain some minimal working examples to do these things just to get users of for a flying start?

I have examples ready if you are interested I could make a pull-request with these for the project?

Kind regards
Thomas Jansson

closest_point: not actually outside the polygon

When you run closest_point, the point it returns is actually on a edge of the polygon. This does not seem to have an effect when using shortest_path or update, but must be fixed. Need to extend the line just slightly to get the point outside of the polygon. See test_closest_point.

Support method to combine overlapping polygons in a Graph

overlaps = True
while overlaps:
overlaps = False
for each polygon n:
for each polygon m:
# https://techoverflow.net/2017/02/23/computing-bounding-box-for-a-list-of-coordinates-in-python/
# checking bounding boxes is a quick way to see if two polygons potentially overlap.
if the bounding box of n overlaps with the bounding box of m:
# combine the polygons or return both if they don't actually overlap
https://stackoverflow.com/questions/2667748/how-do-i-combine-complex-polygons
overlaps = True

Worst case this runs in O(n^3), i.e. all polygons overlap.

polygon_crossing misses a few cases

I narrowed down a mistake in evaluating visible vertices to a failure to detect a point inside a polygon. In the function polygon_crossing there is a section intended to evaluate edges where one or both edge points sit on the ray from the point to +INF in the X direction. It fails in some cases.

In the picture below, the red line indicates two points mistakenly identified as visible to each other. In the process, a point inside the T-shaped polygon along the red line shown is mistakenly identified as outside the polygon.
bugb

I never found the bug exactly, but I came up with some slightly reworked logic that seems functional so far. The reworked logic is as follows:

  • When both edge points lie on the ray, do nothing
  • When just one edge point is on the ray, pick the other edge point and figure out if it lies above or below the ray (in Y dimension). Count total above-ray cases and below-ray cases in two separate counters.
    -After every edge was evaluated, see if both counters are odd (I think it is actually sufficient to see if either is odd, based on assumptions on the polygon being properly formed). If so, this counts as a crossing overall, so increment the crossing counter by 1.

is haversine distance in use?

Since this package seems to be intended for large scale navigation on earth I wonder if the input polygons are given in (latitude, longitude) coordinates. In that case the distance the actual distance between two points would not be euclidian but the haversine distance right? I am curious how this issue has been incorporated.
I ask because my packages https://github.com/MrMinimal64/extremitypathfinder and https://github.com/MrMinimal64/timezonefinder are related to this issue. Thanks

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.