partial pixel rasterization

_partial_pixel_rasterization.md

This is a proof-of-concept for a numpy/rasterio/shapely based implementation of partial coverage rasterization. It works, barely.

The current GDAL algorithms allow for two methods, both binary: the default centroid method and the all-touched method.

This is a third alternative which provides the percentage of coverage of each cell from 0 to 100 which can be thought of as pixel weights for many spatial analyses.

See discussion at rasterio/rasterio#232

alltouched.png

benchmark_results.txt

Performance is not great but at least it doesn't increase exponentially with array size.
It's probably dependent on the number of exterior cells as that is the slow python loop. 

Scale = 10 , shape = (50, 50)
rasterize           0.658988953 ms
rasterize exterior  0.497102737 ms
rasterize pctcover  39.159059525 ms (59x slower)

Scale = 100 , shape = (500, 500)
rasterize           2.088069916 ms
rasterize exterior  1.469135284 ms
rasterize pctcover  363.667011261 ms (174x slower)

Scale = 1000 , shape = (5000, 5000)
rasterize           62.785863876 ms
rasterize exterior  43.768882751 ms
rasterize pctcover  3739.351987839 ms (59x slower)

Scale = 2500 , shape = (12500, 12500)
rasterize           349.674940109 ms
rasterize exterior  345.547199249 ms
rasterize pctcover  10565.802097321 ms (30x slower)

create_data.py

import numpy as np
import json
import rasterio
from affine import Affine

data = np.random.rand(5, 5) * 1000
profile = {'tiled': False,
           'driver': 'GTiff',
           'affine': Affine(1.0, 0.0, 0.0, 0.0, -1.0, 5.0),
           'height': 5,
           'width': 5,
           'dtype': 'float64',
           'count': 1}

with rasterio.open('test.tif', 'w', **profile) as dst:
    dst.write(data, 1)

fc = {
    'type': 'FeatureCollection',
    'features': [
        {'type': 'Feature',
         'properties': {},
         'geometry': {
            "type": "Polygon",
            "coordinates": [
                    [[2.0, 0.0],
                     [2.4, 0.0],
                     [2.4, 0.6],
                     [2.6, 0.6],
                     [2.6, 0.0],
                     [4.0, 0.0],
                     [4.0, 5.0],
                     [0.0, 4.2],
                     [0.0, 2.0],
                     [2.0, 0.0]]]}}]}

with open('test.geojson', 'w') as fh:
    fh.write(json.dumps(fc))

default.png

diff.png

exterior.png

pctcover.png

rasterize.py

import rasterio
import fiona
import numpy as np
from rasterio import features
from affine import Affine
from shapely.geometry import shape, box


def _rasterize_geom(geom, shape, affinetrans, all_touched):
    indata = [(geom, 1)]
    rv_array = features.rasterize(
        indata,
        out_shape=shape,
        transform=affinetrans,
        fill=0,
        all_touched=all_touched)
    return rv_array


def rasterize_pctcover(geom, atrans, shape):
    alltouched = _rasterize_geom(geom, shape, atrans, all_touched=True)
    exterior = _rasterize_geom(geom.exterior, shape, atrans, all_touched=True)

    # Create percent cover grid as the difference between them
    # at this point all cells are known 100% coverage,
    # we'll update this array for exterior points
    pctcover = (alltouched - exterior) * 100

    # loop through indicies of all exterior cells
    for r, c in zip(*np.where(exterior == 1)):

        # Find cell bounds, from rasterio DatasetReader.window_bounds
        window = ((r, r+1), (c, c+1))
        ((row_min, row_max), (col_min, col_max)) = window
        x_min, y_min = (col_min, row_max) * atrans
        x_max, y_max = (col_max, row_min) * atrans
        bounds = (x_min, y_min, x_max, y_max)

        # Construct shapely geometry of cell
        cell = box(*bounds)

        # Intersect with original shape
        cell_overlap = cell.intersection(geom)

        # update pctcover with percentage based on area proportion
        coverage = cell_overlap.area / cell.area
        pctcover[r, c] = int(coverage * 100)

    return pctcover


if __name__ == "__main__":
    with fiona.open("test.geojson") as src:
        geom = shape(next(src)['geometry'])  # take first feature's shapely geometry

    scale = 20
    atrans = Affine(1.0 / scale, 0.0, 0.0, 0.0, -1.0 / scale, 5.0)
    shape = (5 * scale, 5 * scale)
    profile = {
        'affine': atrans,
        'height': shape[0],
        'width': shape[1],
        'count': 1,
        'crs': {},
        'driver': 'GTiff',
        'dtype': 'uint8',
        'nodata': None,
        'tiled': False}

    pctcover = rasterize_pctcover(geom, atrans=atrans, shape=shape)
    with rasterio.open('percent_cover.tif', 'w', **profile) as dst:
        dst.write(pctcover, 1)

test.geojson