Rasterizer



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]:
[[(5, 100), (5, 101), (4, 101), (4, 100), (5, 100)],
 [(5, 100), (5, 101), (4, 101), (4, 100), (5, 100)]]

In [28]:
rizer = Rasterizer(1000,1000)

In [29]:
list(rizer.rasterize_polygon(poly_pix, 1))


Out[29]:
[(100, 4, 4, 1), (100, 5, 5, 1)]

In [ ]: