In [23]:
import rasterio
import numpy
import math
In [24]:
# Ripped from gdalrasterize.cpp
class Rasterizer:
"""Rasterizer rasterizes a geometry and yields for every pixel or scanline rasterized"""
def __init__(self, xsize, ysize):
"""yields for point( xpix, ypix, value)
yields for scanline ( ypix, fromxpix, toxpix, value)"""
self._xsize = xsize
self._ysize = ysize
def rasterize_line(self, points, value = None):
"""Rasterizes an array of LineString points in pixel coords.
If value is not present third ordinate is interpolated and used as value.
yields ( xpix, ypix, value)"""
# ported from GDALdllImageLine in llrasterize.cpp
if points is None or len(points) == 0:
return
if value is None and len(points[0]) < 3:
raise Exception("No value specified and coordinates are two-dimensional")
usevalue = not value is None
for ix in xrange( 1, len(points) ):
x = int( math.floor( points[ ix - 1][0]) )
y = int( math.floor( points[ ix - 1][1]) )
x1 = int( math.floor( points[ ix ][0]) )
y1 = int( math.floor( points[ ix ][1]) )
if usevalue:
v = value
v1 = value
else:
v = points[ ix - 1 ][2]
v1 = points[ ix ][2]
dx = abs( x - x1 )
dy = abs( y - y1 )
xstep = -1 if x > x1 else 1
ystep = -1 if y > y1 else 1
if dx >= dy:
xerr = dy << 1
yerr = xerr - (dx << 1)
err = xerr - dx
dv = 0 if (dx == 0 or usevalue) else (v1 - v) / float(dx)
while dx >= 0 :
if 0 <= x and x < self._xsize and 0 <= y and y < self._ysize:
yield ( x, y, v )
v += dv
x += xstep
if err > 0:
y += ystep
err += yerr
else:
err += xerr
dx -= 1
else:
# dy > dx
xerr = dx << 1
yerr = xerr - (dy << 1)
err = xerr - dy
dv = 0 if (dy == 0 or usevalue) else (v1 - v) / float(dy)
while dy >= 0:
if 0 <= x and x < self._xsize and 0 <= y and y < self._ysize:
yield ( x, y, v)
v += dv
y += ystep
if err > 0:
x += xstep
err += yerr
else:
err += xerr
dy -= 1
def rasterize_polygon( self, rings, value):
"""( ypix, fromxpix, toxpix, value)"""
if value is None:
raise Exception("Cannot rasterize polygon using Z")
if rings is None or len(rings) <= 0:
return
n = 0
miny = rings[0][0][1]
maxy = miny
for ring in rings:
n += len(ring)
for point in ring:
if miny > point[1]:
miny = point[1]
if maxy < point[1]:
maxy = point[1]
miny = int( miny )
maxy = int( maxy )
miny = max( miny, 0 )
maxy = min( maxy, self._ysize - 1 )
minx = 0
maxx = self._xsize - 1
for y in xrange(miny, maxy + 1):
dy = y + 0.5 # Center height of line
part = 0
polyInts = []
for ring in rings:
for ind2 in xrange(1, len(ring)):
ind1 = ind2 - 1
dy1 = ring[ind1][1]
dy2 = ring[ind2][1]
if (dy1 < dy and dy2 < dy) or (dy1 > dy and dy2 > dy):
continue
if dy1 < dy2:
dx1 = ring[ind1][0]
dx2 = ring[ind2][0]
elif dy1 > dy2:
dy2 = ring[ind1][1]
dy1 = ring[ind2][1]
dx2 = ring[ind1][0]
dx1 = ring[ind2][0]
else:
# DO NOT skip bottom horizontal segments
# -Fill them separately-
# They are not taken into account twice.*/
if ring[ind1][0] > ring[ind2][0]:
horizontal_x1 = int( math.floor( ring[ind2][0] + 0.5 ) )
horizontal_x2 = int( math.floor( ring[ind1][0] + 0.5 ) )
if horizontal_x1 > maxx or horizontal_x2 <= minx:
continue
yield (y, horizontal_x1, horizontal_x2 - 1, value)
else:
# skip top horizontal segments (they are already filled in the regular loop)
continue
if dy < dy2 and dy >= dy1:
intersect = ( dy - dy1 ) * ( dx2 - dx1 ) / ( dy2 - dy1 ) + dx1
polyInts.append( int( math.floor( intersect ) ) )
polyInts.sort()
for i in xrange( 0, len(polyInts), 2):
if polyInts[i] <= maxx and polyInts[i + 1] > minx:
xstart = polyInts[i]
xend = polyInts[i+1]
if xstart > maxx or xend < 0:
continue
if xstart < 0:
xstart = 0
if xend > maxx:
xend = maxx
yield ( y, xstart, xend, value)
In [25]:
poly = { "type": "Polygon",
"coordinates": [
[ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0] ],
[ [100.2, 0.2], [100.8, 0.2], [100.8, 0.8], [100.2, 0.8], [100.2, 0.2] ]
]
}
In [26]:
raster = rasterio.open('./dummy.asc')
In [30]:
poly_pix = [[raster.index(*coord) for coord in ring] for ring in poly['coordinates']]
poly_pix
Out[30]:
In [28]:
rizer = Rasterizer(1000,1000)
In [29]:
list(rizer.rasterize_polygon(poly_pix, 1))
Out[29]:
In [ ]: