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')
In [5]:
deccanpoints
Out[5]:
In [6]:
coastpoints
Out[6]:
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 [ ]: