From 2b9ce6dcf0ed5c77ff44171ec2295e7e7077442e Mon Sep 17 00:00:00 2001
From: Tsunghao Huang <j22700126@gmail.com>
Date: Tue, 16 Jan 2018 10:53:00 +0100
Subject: [PATCH] Add files via upload

---
 DistanceCalculator.py | 74 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 74 insertions(+)
 create mode 100644 DistanceCalculator.py

diff --git a/DistanceCalculator.py b/DistanceCalculator.py
new file mode 100644
index 0000000..7a385eb
--- /dev/null
+++ b/DistanceCalculator.py
@@ -0,0 +1,74 @@
+import gdal, osr
+from skimage.graph import route_through_array
+import numpy as np
+from geopy.distance import vincenty
+import matplotlib.pyplot as plt
+import pandas as pd
+
+#Transforming a raster map to array datatype.
+raster = gdal.Open('~/map.tif')
+band = raster.GetRasterBand(1)
+mapArray = band.ReadAsArray()
+
+#Get geotransform information and declare some variables for later use
+geotransform = raster.GetGeoTransform()
+originX = geotransform[0]
+originY = geotransform[3] 
+pixelWidth = geotransform[1] 
+pixelHeight = geotransform[5]
+
+ports = pd.read_csv('~/port.csv')
+
+#Visualize the map base on the array, if you want, not neccessary.
+plt.imshow(mapArray)
+plt.gray()
+plt.show()
+
+#transform the coordinates to the exact position in the array.
+def coord2pixelOffset(x,y):
+    
+    xOffset = int((x - originX)/pixelWidth)
+    yOffset = int((y - originY)/pixelHeight)
+    return xOffset,yOffset
+
+#create a path which travels through the cost map.
+def createPath(costSurfaceArray,startCoord,stopCoord):   
+
+    # coordinates to array index
+    startCoordX = startCoord[0]
+    startCoordY = startCoord[1]
+    startIndexX,startIndexY = coord2pixelOffset(startCoordX,startCoordY)
+
+    stopCoordX = stopCoord[0]
+    stopCoordY = stopCoord[1]
+    stopIndexX,stopIndexY = coord2pixelOffset(stopCoordX,stopCoordY)
+
+    # create path
+    indices, weight = route_through_array(costSurfaceArray, (startIndexY,startIndexX), (stopIndexY,stopIndexX),geometric=True,fully_connected=True)
+    indices = np.array(indices).T
+    indices = indices.astype(float)
+    indices[1] = indices[1]*pixelWidth + originX
+    indices[0] = indices[0]*pixelHeight + originY
+    return indices
+
+#Calculate the vincenty distance starts from the first pair of points to the last.
+def calculateDistance(pathIndices):
+    distance = 0
+    for i in range(0,(len(pathIndices[0])-1)):
+        distance += vincenty((pathIndices[1,i], pathIndices[0,i]), (pathIndices[1,i+1], pathIndices[0,i+1])).miles*0.868976
+    return distance
+
+
+def distanceCalculator(startCoord, stopCoord):
+    pathIndices = createPath(mapArray,startCoord,stopCoord)
+    distance = calculateDistance(pathIndices)
+    print(distance)
+    return distance
+
+
+#Test, in this case it's from CATANIA to MANTA
+ports.Port[1000]
+startCoord = (ports.longitude[2205], ports.latitude[2205])
+stopCoord = (ports.longitude[1000], ports.latitude[100])
+
+distanceCalculator(startCoord, stopCoord)