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
Verify that fwd and inverse are actually correct
In [2]:
costFile = '../Shapefiles/TotalAreaSlopeCost.tif'
totalCostFile = '../Shapefiles/CompleteCost.tif'
px2Coord = pixelToCoordFn(costFile)
coord2Px = coordToPixelFn(costFile)
coord2Px(px2Coord((1001,1011)))
px2Coord(coord2Px((72,81)))
Out[2]:
We read different kinds of geometry from a shapefile
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]
coordinates.append(coords)
elif feature is 'MultiPolygon':
coords = p['geometry']['coordinates'][0][0]
print coords
coordinates.append(coords)
if 'Name' in p['properties'].keys():
names.append(p['properties']['Name'])
return 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]:
names, coords = readFeature('../Shapefiles/SahyadriPasses.shp')
polyNames, polyCoords = readFeature('../Shapefiles/WaterBodies.shp', feature='Polygon')
riverNames, riverCoords = readFeature('../Shapefiles/DeccanRiverBuffer.shp', feature='Polygon')
In [5]:
from PIL import Image, ImageDraw
import numpy as np
import scipy as sp
def drawPolygons(polyList, rasterFile, fillValue=4, fillType='uniform'):
""" Draw water bodies (or any other polygon on an image)
:param polyList: list of polygons to draw
:param fillValue: The value to use at the filled locations
:param rasterFile: Raster from which dimensions are extracted
:param fillType: either 'uniform' or 'random'
"""
with rasterio.open(rasterFile) as raster:
width = raster.width
height = raster.height
img = Image.new('L',(height, width), 0)
coord2px = coordToPixelFn(rasterFile)
for poly in polyList:
polyInCoords = poly
polyInPixels = map(coord2Px, polyInCoords)
ImageDraw.Draw(img).polygon(polyInPixels, outline=True, fill=True)
img = np.array(img,dtype=float).transpose()
if fillType == 'uniform':
return fillValue*img
elif fillType == 'random':
print 'Filling Random values.'
indices = np.nonzero(img)
numVals = len(indices[0])
print 'in',numVals,'locations'
img[indices] = sp.random.exponential(scale=2,size=numVals)
print np.mean(img[indices])
return img
else:
print 'Fill Type not recognised'
return
In [6]:
output = drawPolygons(polyCoords, costFile, fillType='random')
outputRivers = drawPolygons(riverCoords, costFile)
In [7]:
import matplotlib.pyplot as plt
plt.imshow(output)
print 'finished one'
with rasterio.open(costFile) as src:
arr = src.read(1)
print 'finished reading ', arr.shape
#plt.figure()
#plt.imshow(arr)
#plt.figure()
#plt.imshow(outputRivers)
#print 'all done!'
#plt.show()
In [8]:
writeRasterFile(arr+outputRivers+output, totalCostFile, costFile)
In [9]:
np.nonzero(output)
Out[9]:
In [ ]: