Feature Request: Create Numpy Mask from Geometry about rasterio

[sgillies]

@brendan-ward Yes, rasterizing would be a good feature.

I think we could simplify it quite a bit by having the same constraint as gdal_rasterize, that the vectors and raster have the same CRS, no reprojection. Also, I think GeoJSON-ish Python objects (the kind of things you could pass to shapely.geometry.shape()) are a better fit for rasterio than Shapely objects.

What would you think of adding a rasterio.features.rasterize() function that is symmetric with features.shapes()? Usage something like:

from rasterio.features import rasterize, shapes

# get the shapes and values from "source" and rasterize them into "destination".
destination = rasterize(shapes(source, transform=None), source.shape, transform=None)

the above should produce a destination array equal to the source. A round trip. Rasterio is all about symmetry :)

Of course, the shapes above would be in image coordinates. Passing not-None transforms in would let us use coordinates appropriate to a CRS.

[sgillies]

@brendan-ward Most importantly: pull request would be very welcome! Feel free to grab code from rasterio/_features.pyx.

[sgillies]

Sorry, all. I'm going to keep these issues on track by deleting tangential comments.

[brendan-ward]

Do we want to broaden this beyond simple binary masks, and allow parameters on this function to allow creation of either a mask, or rasterize values of a given property on the features?

I was also planning to add this in _features.pyx, since this seems highly related, rather than creating a new file. Sound good?

[brendan-ward]

Progress! I do think there has to be a better way to do some of these things - like looping over the features.

def _rasterize(features, size, transform=None):
    """
    :param features: fiona style feature iterator
    :param size: (rows, cols)
    :param transform: GDAL style transform.  If provided, will be set on output
    """

    cdef int retval, rows, cols
    cdef size_t i, num_features
    cdef void *memdriver
    cdef void *out_ds
    cdef void *out_band
    cdef double geotransform[6]
    cdef double pixel_values[1]
    cdef int dst_bands[1]

    pixel_values[0] = 1 #TODO: Allow other pixel values

    rows, cols = size

    #Do the boilerplate required to create a band
    memdriver = _gdal.GDALGetDriverByName("MEM")
    if memdriver == NULL:
        raise ValueError("NULL driver for 'MEM'")
    out_ds = _gdal.GDALCreate(memdriver, "output", cols, rows, 1, <_gdal.GDALDataType>1, NULL) #TODO: revisit data type
    if out_ds == NULL:
        raise ValueError("NULL output datasource")

    if transform:
        for i in range(6):
            geotransform[i] = transform[i]
        err = _gdal.GDALSetGeoTransform(out_ds, geotransform)
        if err:
            raise ValueError("transform not set: %s" % transform)
    out_band = _gdal.GDALGetRasterBand(out_ds, 1)
    if out_band == NULL:
        raise ValueError("NULL output band")
    dst_bands[0] = 1

    #TODO: figure out a cleaner way to do this
    features_json = []
    for feature in features:  #have to loop over features, since it may be yielded from a generator
        features_json.append(json.dumps(feature['geometry']))
    num_features = len(features_json)

    cdef void **ogr_geoms = <void **>_gdal.CPLMalloc(num_features * sizeof(void*))
    try:
        for i in range(num_features):
            ogr_geoms[i] = _ogr.OGR_G_CreateGeometryFromJson(features_json[i])

        #TODO: add options
        retval = _gdal.GDALRasterizeGeometries(out_ds, 1, dst_bands, 1, ogr_geoms, NULL, geotransform, pixel_values, NULL, NULL, NULL)
        out = np.zeros(size, np.uint8) #FIXME - data type
        retval = io_ubyte(out_band, 0, 0, 0, cols, rows, out) #FIXME - data type

    finally:
        _gdal.CPLFree(ogr_geoms)

    if out_ds != NULL:
        _gdal.GDALClose(out_ds)

    return out

[sgillies]

Exciting!