How to use the rasterio.features.rasterize function in rasterio

To help you get started, we’ve selected a few rasterio examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github mapbox / rasterio / tests / test_features.py View on Github external
def test_rasterize_fill_value(basic_geometry, basic_image_2x2):
    """All pixels not covered by shapes should be given fill value."""
    default_value = 2
    assert np.array_equal(
        basic_image_2x2 + 1,
        rasterize(
            [basic_geometry], out_shape=DEFAULT_SHAPE, fill=1,
            default_value=default_value
        )
github mapbox / rasterio / tests / test_features.py View on Github external
def test_rasterize_missing_out(basic_geometry):
    """If both out and out_shape are missing, should raise exception."""
    with pytest.raises(ValueError):
        rasterize([basic_geometry], out=None, out_shape=None)
github mapbox / rasterio / tests / test_features.py View on Github external
def test_rasterize_geomcollection_no_hole():
    """
    Make sure that bug reported in
    https://github.com/mapbox/rasterio/issues/1253
    does not recur.  GeometryCollections are flattened to individual parts,
    and should result in no holes where parts overlap.
    """

    geomcollection = {'type': 'GeometryCollection', 'geometries': [
        {'type': 'Polygon',
            'coordinates': (((0, 0), (0, 5), (5, 5), (5, 0), (0, 0)),)},
        {'type': 'Polygon',
            'coordinates': (((2, 2), (2, 7), (7, 7), (7, 2), (2, 2)),)}
    ]}

    expected = rasterize(geomcollection['geometries'], out_shape=DEFAULT_SHAPE)

    assert np.array_equal(
        rasterize([geomcollection], out_shape=DEFAULT_SHAPE),
        expected
    )
github mapbox / rasterio / tests / test_features.py View on Github external
def test_rasterize_geometries_symmetric():
    """Make sure that rasterize is symmetric with shapes."""
    transform = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)
    truth = np.zeros(DEFAULT_SHAPE, dtype=rasterio.ubyte)
    truth[2:5, 2:5] = 1
    s = shapes(truth, transform=transform)
    result = rasterize(s, out_shape=DEFAULT_SHAPE, transform=transform)
    assert np.array_equal(result, truth)
github geowurster / gj2ascii / gj2ascii / core.py View on Github external
# This potentially creates a large in-memory object so if processing an entire layer it is
    # best to explicitly define min/max, especially because its also faster.
    if bbox:
        x_min, y_min, x_max, y_max = bbox
    else:
        _bbox, ftrz = min_bbox(ftrz, return_iter=True)
        x_min, y_min, x_max, y_max = _bbox

    x_delta = x_max - x_min
    y_delta = y_max - y_min
    cell_size = x_delta / width
    height = int(y_delta / cell_size)
    if height is 0:
        height = 1

    output_array = rasterize(
        fill=0,
        default_value=1,
        shapes=(g for g in _geometry_extractor(ftrz)),
        out_shape=(height, width),
        transform=affine.Affine.from_gdal(*(x_min, cell_size, 0.0, y_max, 0.0, -cell_size)),
        all_touched=all_touched,
        dtype=rio.uint8
    )

    # Convert to string dtype and do character replacements
    output_array = output_array.astype(np.str_)
    if fill is not None and fill != '0':
        output_array = np.char.replace(output_array, '0', fill)
    if char is not None and fill != '1':
        output_array = np.char.replace(output_array, '1', char)
github CosmiQ / cw-geodata / cw_geodata / vector_label / mask.py View on Github external
shape = reference_im.shape
        if do_transform:
            affine_obj = reference_im.transform

    # extract geometries and pair them with burn values
    if burn_field:
        if out_type == 'int':
            feature_list = list(zip(df[geom_col],
                                    df[burn_field].astype('uint8')))
        else:
            feature_list = list(zip(df[geom_col],
                                    df[burn_field].astype('uint8')))
    else:
        feature_list = list(zip(df[geom_col], [burn_value]*len(df)))

    output_arr = features.rasterize(shapes=feature_list, out_shape=shape,
                                    transform=affine_obj)
    if out_file:
        meta = reference_im.meta.copy()
        meta.update(count=1)
        if out_type == 'int':
            meta.update(dtype='uint8')
        with rasterio.open(out_file, 'w', **meta) as dst:
            dst.write(output_arr, indexes=1)

    return output_arr
github LSDtopotools / LSDMappingTools / pythonic_rasterizer.py View on Github external
tgdf[args.lookup_field][tgdf[field] == key] = dictionary[key]  
            
        # Get the template raster
        rst = rasterio.open(rst_fn)

        # copy and update the metadata from the input raster for the output
        meta = rst.meta.copy()
        meta.update(compress='lzw')

        with rasterio.open(out_fn, 'w+', **meta) as out:
            out_arr = out.read(1)

            # this is where we create a generator of geom, value pairs to use in rasterizing
            shapes = ((geom,value) for geom, value in zip(tgdf.geometry, tgdf[args.lookup_field]))

            burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=out.transform)
            out.write_band(1, burned)

    else:
        out_fn = this_dir+args.fname_prefix+".tif"
        
        print("I am going to proceed using unique values")
        print("This will relsult in a new shapefile")
        # creating the rasterization column
        tgdf['RZ_key'] = '' 

        #Hosts the equivalences
        eqdic = {}

        # Fill the null values
        tgdf[field].fillna(-5, inplace=True)
github yghlc / Landuse_DL / sentinelScripts / get_subImages.py View on Github external
# add the center polygons to adj_polygons
    adj_polygons.extend([center_polygon])
    adj_polygons_class.extend([class_int])

    with rasterio.open(sub_image_path) as src:

        transform = src.transform

        burn_out = np.zeros((src.height, src.width))

        # rasterize the shapes
        burn_shapes = [(item_shape, item_class_int) for (item_shape, item_class_int) in
                       zip(adj_polygons, adj_polygons_class)]
        #
        out_label = rasterize(burn_shapes, out=burn_out, transform=transform,
                              fill=0, all_touched=False, dtype=rasterio.uint8)

        # test: save it to disk
        kwargs = src.meta
        kwargs.update(
            dtype=rasterio.uint8,
            count=1)
        #   width=burn_out.shape[1],
        #    height=burn_out.shape[0],#transform=transform
        # remove nodta in the output
        if 'nodata' in kwargs.keys():
            del kwargs['nodata']

        with rasterio.open(save_path, 'w', **kwargs) as dst:
            dst.write_band(1, out_label.astype(rasterio.uint8))
github GeoscienceAustralia / wagl / wagl / brdf.py View on Github external
bound_poly = ops.transform(lambda x, y: dst_geotransform * (x, y), box(0., 0., ds_width, ds_height, ccw=False))
    if not bound_poly.intersects(dst_poly):
        return BrdfTileSummary.empty()

    ocean_poly = ops.transform(lambda x, y: fid_mask.transform * (x, y), box(0., 0., fid_mask.width, fid_mask.height))
    if not ocean_poly.intersects(dst_poly):
        return BrdfTileSummary.empty()

    # read ocean mask file for correspoing tile window
    # land=1, ocean=0
    bound_poly_coords = list(bound_poly.exterior.coords)[:4]
    ocean_mask, _ = read_subset(fid_mask, *bound_poly_coords)
    ocean_mask = ocean_mask.astype(bool)

    # inside=1, outside=0
    roi_mask = rasterize([(dst_poly, 1)], fill=0, out_shape=(ds_height, ds_width), transform=dst_geotransform)
    roi_mask = roi_mask.astype(bool)

    # both ocean_mask and mask shape should be same
    if ocean_mask.shape != roi_mask.shape:
        raise ValueError('ocean mask and ROI mask do not have the same shape')
    if roi_mask.shape != ds.shape:
        raise ValueError('BRDF dataset and ROI mask do not have the same shape')

    roi_mask = roi_mask & ocean_mask

    def layer_sum(param):
        layer = ds[param][:, :]
        common_mask = roi_mask & (layer != ds.attrs['_FillValue'])
        layer = layer.astype('float32')
        layer[~common_mask] = np.nan
        layer = ds.attrs['scale_factor'] * (layer - ds.attrs['add_offset'])