The routing code



In [1]:
# %load transformsUtils.py
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
    """

    print 'reading data'
    with rasterio.open(costSurfaceFile,'r') as costFile:
        costSurfaceArray = costFile.read() 
        # 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]:
# %load readWriteUtils.py
import rasterio
import fiona

def readFeature(shapeFile, feature='Point'):
    """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]:
coastpoints = readFeature('../Shapefiles/CoromandelCoastPlaces.shp')
deccanpoints = readFeature('../Shapefiles/InteriorPlaces.shp')


Number of features:  28
Number of features:  58

In [5]:
deccanpoints


Out[5]:
{u'Adam': (79.44988317286726, 20.99930153541378),
 u'Adilabad': (78.52950372547053, 19.66824489882036),
 u'Alluru': (80.44099004685778, 16.76914862020891),
 u'Amaravathi': (80.35748150779486, 16.57265804076962),
 u'Aurangabad': (75.34127720572363, 19.86848230793264),
 u'Badvel': (79.05858948334388, 14.72978737897093),
 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'Daulatabad': (75.21303051920609, 19.94264550944658),
 u'Gajulabanda': (79.44493733632086, 17.44215256110829),
 u'Garikapadu': (79.7007591234751, 16.33464846208296),
 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'Hubli-Dharwad': (75.08483583551289, 15.36121334290463),
 u'Hyderabad': (78.48684875152986, 17.38030290599092),
 u'Jaggayyapeta': (80.09765382462204, 16.89172811164058),
 u'Jangaon': (79.17975332696246, 17.71930952888685),
 u'Jeypore': (82.57296241810019, 18.86498408451859),
 u'Kadiri': (78.16024257463974, 14.11123859788977),
 u'Karad': (74.20027654621495, 17.27562920444157),
 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'Nizamabad': (78.0929612153006, 18.6674997429364),
 u'Paithan/Pratisthana': (75.37984949409807, 19.47967189425886),
 u'Parbhani': (76.77932842138735, 19.26734494918001),
 u'Pauni': (79.63662767566476, 20.79293732145555),
 u'Pedamadduru': (80.39783444460723, 16.53820612323536),
 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'Tadipatri': (77.99042277711548, 14.90442253101064),
 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'Gudivada': (80.98916947371474, 16.42483255880556),
 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'Vijayaderpuram (Vijayawada)': (80.64722940849252, 16.50547588848681),
 u'Yerrampalem': (81.87762278372801, 17.22616295399498)}

Actual Routing steps follow...


In [7]:
startPoint = coastpoints['Bhattiprolu']
#endPoint = deccanpoints['Nizamabad']
costFile = '../Shapefiles/CompleteCost.tif'

In [ ]:
with rasterio.open(costFile) as src:
    costArray = src.read(1)
    
graph = MCP_Geometric(costArray, fully_connected=True)

In [ ]: