9

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.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
user3367601
  • 163
  • 1
  • 7
  • 1
    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 '15 at 01: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. – radouxju Aug 20 '15 at 13:57

1 Answers1

4

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()
graph.add_edge('coord_a', 'coord_b', weight=distance)
# add more edges in the same way
start_node = find_closest_node(graph, lat_start, long_start)
end_node = find_closest_node(graph, lat_end, long_end)
ship_path = networkx.shortest_path(graph, start_node, end_node)
Community
  • 1
user55937
  • 1,293
  • 11
  • 18
  • Thanks a lot. That is just the kind of help I needed. Will try this on my developer machine next week. Will request for help if I get stuck. Thanks a lot. – user3367601 Aug 21 '15 at 16:30
  • Great! Let me know if I can help. It would be good to include any additions in the answer. Also, if you don't mind shelling out some money, here's a good resource where I first learned some of this (I'm not the author): https://www.packtpub.com/big-data-and-business-intelligence/python-geospatial-analysis-essentials. It's pretty concise, and has a lot of code samples for geospatial stuff, shortest paths, and map matching. – user55937 Aug 21 '15 at 17:09