## The routing code

``````

In [1]:

import rasterio
import numpy as np
from affine import Affine

def pixelToCoordFn(raster):
""" Converts from pixel indices to actual coordinates.

:param raster: file for which the coordinates have to be calculated
:param x: x coordinate
:param y: y coordinate
:returns transFunc: function which does the transformation

"""

with rasterio.open(raster, 'r') as rasterFile:
T0 = rasterFile.affine

#By default, transformation is for pixel edge. convert to centre
T1 = T0#*Affine.translation(0.5,0.5)

transFunc = lambda (r, c): (c, r) * T1

return transFunc

def coordToPixelFn(raster):
""" Converts from actual coordinates to pixel indices.

:param raster: file for which the coordinates have to be calculated
:param x: x coordinate
:param y: y coordinate
:returns transFunc: function which does the transformation

"""

with rasterio.open(raster, 'r') as rasterFile:
T0 = rasterFile.affine

a1 = 1./T0.a
b1 = T0.b
c1 = -T0.c/T0.a

d1 = T0.d
e1 = 1/T0.e
f1 = -T0.f/T0.e

#print T0.a, T0.b, T0.c, T0.d, T0.e, T0.f
T0 = Affine(a1,b1,c1,d1,e1,f1)
#print T0.a, T0.b, T0.c, T0.d, T0.e, T0.f
#print T0*(69.99958357546711, 20.000417005657994)

#By default, transformation is for pixel edge. convert to centre
T1 = T0#*Affine.translation(0.5,0.5)

transFunc = lambda (r, c): tuple(np.array((r, c)*T1)[::-1])

return transFunc

``````
``````

In [2]:

from skimage.graph import MCP_Geometric

class returnCost(object):
"""
Simple class to store the graph and its
cumulative cost
"""
def __init__(self, graph, cumulativeCost, name, px2Coord, coord2Px):

self.graph = graph
self.cumCost = cumulativeCost
self.name = name
self.px2Coord = px2Coord
self.coord2Px = coord2Px
return self

def returnRouteTo(destination):

index = self.coord2Px(destination)
route = self.graph.traceback(index)
coords = [self.px2Coord(point) for point in route]
return (coords, self.cumCost[index])

def createPathFrom(startPlaceName,costSurfaceFile,startCoord):
"""Create a least cost path from startCoord to stopCoord given a cost surface.

:param costSurfaceFile: name of file containing costs
:param startCoord: (x,y) values in world coordinates to start from
:param startPlaceName: Name of the place we are starting from

:returns: an object which gives the costs and route given an end point
"""

with rasterio.open(costSurfaceFile,'r') as costFile:
# creates array from cost surface raster

#
px2Coord = pixelToCoordFn(costSurfaceFile)
coord2Px = coordToPixelFn(costSurfaceFile)

# coordinates to array index
startIndices = coord2Px(startCoord)

print 'Starting from: ', startIndices
print 'Initialising'
# create path
graph = MCP_Geometric(costSurfaceArray, fully_connected=True)
print 'done creating graph, now calculating all routes'
cumCost, traceback = graph.find_costs([startIndices])
print done

costObject = returnCost(graph, cumCost, startPlaceName, px2Coord, coord2Px)

return costObject

``````
``````

In [3]:

import rasterio
import fiona

"""Read all the points and their properties from a vector shape file

:param shapeFile: name of a file to read from
:param feature: either 'Point', 'Linestring' or 'Polygon' for now
:returns: a list of feature coordinates and their names

"""

with fiona.open(shapeFile) as shpFile:
allFeatures = list(shpFile)

coordinates = []
names = []

print 'Number of features: ', len(allFeatures)
for p in allFeatures:
if p['geometry']['type'] == feature:

if feature is 'Polygon':
coords = p['geometry']['coordinates'][0]
if len(coords) == 3:
# We want only the first two coordinates
coords = coords[0:-1]
coordinates.append(coords)
elif feature is 'MultiPolygon':
coords = p['geometry']['coordinates'][0][0]
if len(coords) == 3:
# We want only the first two coordinates
coords = coords[0:-1]
coordinates.append(coords)
else:
coords = p['geometry']['coordinates']
if len(coords) == 3:
# We want only the first two coordinates
coords = coords[0:-1]
coordinates.append(coords)

if 'Name' in p['properties'].keys():
names.append(p['properties']['Name'])

return dict(zip(names, coordinates))

def writeRasterFile(outputArray, outputFilename, refRaster):

with rasterio.open(refRaster) as src:
kwargs = src.meta

with rasterio.open(outputFilename, 'w', **kwargs) as dst:
dst.write_band(1, outputArray.astype(kwargs['dtype']))

``````
``````

In [4]:

``````
``````

Number of features:  28
Number of features:  58

``````
``````

In [5]:

deccanpoints

``````
``````

Out[5]:

u'Alluru': (80.44099004685778, 16.76914862020891),
u'Amaravathi': (80.35748150779486, 16.57265804076962),
u'Belgaum': (74.49954574080886, 15.84901894017051),
u'Bellary': (76.92117873163025, 15.13823924528715),
u'Bidar': (77.519216811633, 17.91842182127456),
u'Bijapur': (75.7093407630056, 16.82743100851483),
u'Bodhan': (77.90009660906304, 18.66974366553383),
u'Chandrapur': (79.29934030661626, 19.9474431901183),
u'Gajulabanda': (79.44493733632086, 17.44215256110829),
u'Goli': (79.5434675017246, 16.59450020186393),
u'Gulbarga': (76.83396651534765, 17.32551520806241),
u'Guntakal': (77.37276063520018, 15.16365983163976),
u'Guntupalli': (81.13092964834577, 17.01980841825057),
u'Hanamkonda': (79.55700126542254, 18.00451828047292),
u'Hinganghat': (78.84351576485443, 20.54616862239718),
u'Jaggayyapeta': (80.09765382462204, 16.89172811164058),
u'Jangaon': (79.17975332696246, 17.71930952888685),
u'Jeypore': (82.57296241810019, 18.86498408451859),
u'Karimnagar': (79.12776108314662, 18.43201177520683),
u'Koilkuntla': (78.31660254474086, 15.23213338730874),
u'Kolhapur': (74.24471429038874, 16.69068836133482),
u'Kondapur': (78.70924827361775, 18.5078100896339),
u'Koraput': (82.70984909855451, 18.80514157617955),
u'Korukonda': (81.83168629707318, 17.17148761059733),
u'Mancherial': (79.46376933502866, 18.8673716603848),
u'Nagarjunakonda': (79.24199678951538, 16.52183413460003),
u'Nagpur': (79.08789817754764, 21.14436177623766),
u'Nalgonda': (79.26249043624492, 17.05606774654516),
u'Nanded': (77.29983754081528, 19.14899325321492),
u'Nandori': (78.97262228635823, 20.51549885475078),
u'Nandyal': (78.47532948984444, 15.46354714616056),
u'Nashik': (73.78958686068472, 19.99540231835371),
u'Nellakondapalli': (80.06844451703626, 17.11177535790434),
u'Paithan/Pratisthana': (75.37984949409807, 19.47967189425886),
u'Parbhani': (76.77932842138735, 19.26734494918001),
u'Pauni': (79.63662767566476, 20.79293732145555),
u'Phanigiri': (79.47240516099104, 17.42685330178541),
u'Pulivendula': (78.2263479018568, 14.42097483173875),
u'Pune': (73.85642271698597, 18.51794103099387),
u'Raichur': (77.360229, 16.256867),
u'Raipur': (81.62096789818713, 21.21497295599469),
u'Sambalpur': (83.96255252011352, 21.43249096082672),
u'Satara': (74.00081569295679, 17.6905105618208),
u'Vinukonda': (79.74046144589191, 16.04918314969372),
u'Yerrampalem': (81.87762278372801, 17.22616295399498)}

``````
``````

In [6]:

coastpoints

``````
``````

Out[6]:

{u'Anakapalle': (83.00875379544561, 17.68576243311793),
u'Arugolanu': (80.95597753800922, 16.5624674260166),
u'Bhattiprolu': (80.77990331869378, 16.09976381993026),
u'Buddhani': (80.55216526074642, 15.9471966617142),
u'Chinnaganjam': (80.23677804109577, 15.69547066681743),
u'Dharapalem': (82.930650153634, 17.5360386476163),
u'Ganjam': (85.05139343552105, 19.38692022064951),
u'Ghantasala/Kantakasaila': (80.94419671522414, 16.16882457923828),
u'Guntur': (80.43636906697354, 16.3047782626829),
u'Kanuparthi': (80.21399871732065, 15.58640997565835),
u'Kodavali': (82.2752559344329, 17.22965596702113),
u'Koduru': (81.0339529191514, 16.0083183783955),
u'Korukonda': (81.83168629707318, 17.17148761059733),
u'Machilipatnam': (81.12942577890752, 16.16784580293908),
u'Nellore': (79.98633385330119, 14.44214673434968),
u'Pedavegi': (81.10654445718636, 16.81006050148935),
u'Peddaganjam': (80.22737229545297, 15.64196473414722),
u'Pithapuram': (82.2529052307822, 17.11207269213028),
u'Ramateertham': (83.48542029622485, 18.17601006760968),
u'Ramatirtham': (80.14070000135608, 14.64445603145569),
u'Salihundam': (84.04993714992355, 18.33315124578951),
u'Sankaram': (83.02134742034013, 17.68568964751766),
u'Srikakulam': (83.89642756047344, 18.29596163470776),
u'Srikalahasti': (79.69986833292941, 13.75969203264884),
u'Uppugunduru': (80.21710145938552, 15.63952786244812),
u'Yerrampalem': (81.87762278372801, 17.22616295399498)}

``````

## Actual Routing steps follow...

``````

In [7]:

startPoint = coastpoints['Bhattiprolu']
costFile = '../Shapefiles/CompleteCost.tif'

``````
``````

In [ ]:

with rasterio.open(costFile) as src:

graph = MCP_Geometric(costArray, fully_connected=True)

``````
``````

In [ ]:

``````