# Find shortest path between points avoiding shapefile?

I am an intermediate python developer with very limited GIS experience. I have a list of latitude/longitudes of a ship at sea. I need to join them into a polyline using the shortest distance between them to trace the path followed by the ship. However, the path cannot pass over land (for which I have a shapefile). How can I do this? I use Python 2.7.

I am presently not using any GIS software and hope to do the task using the python modules like matplotlib and scipy.

• Could you edit the question to include the software and version (if relevant)? That way any answers can be tailored to your needs. Welcome to GIS.SE!
– Paul
Aug 20, 2015 at 1:36
• I would consider a weighted cost distance analysis in raster (with huge cost on land pixels). Following the coast when your shortest path intersect the land will not give the the shortest path. Aug 20, 2015 at 13:57

This sounds like a fun problem. I would consider:

1. Read the land boundary shp file using one of these recipes: How to install Fiona to read Shapefile attributes with OSGeo4W?. It would probably be good to read it into a Shapely geometry. I like using ogr to read shapefiles:
```import ogr
import shapely
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open('landboundary.shp', 0)
layer = dataSource.GetLayer()
feature = layer[0] # first boundary
geom = feature.GetGeometryRef()
shapely.geometry.base.geom_from_wkt(geom.ExportToWkt())
```
1. Use NetworkX to build a network of all your lat long ship coordinates, but don't connect coordinates whose straight line connection intersects the land boundary. Here's an example to solve for the intersection of geometries: Shapely LineString and Polygon intersect?.

2. Then compute the least cost path using NetworkX with distance as the cost. The LCP won't cross over land because those nodes are not connected.

Also, just fyi, shapely geometries convert to WKT (well known text which is a text-based format) like this, `myShapelyGeometry.wkt`, which is very easy and intuitive to parse and use how you want.

Edit: Just some extra info: To compute the distance between two lat long coords, you can use pyproj:

``````import pyproj
geod = pyproj.Geod(ellps='WGS84')
_,_,distance = geod.inv(long_start,lat_start,long_end,lat_end)
``````

Also, you can start a networkx graph like:

``````import networkx
graph = networkx.Graph()