Feature Operations

def shapes(image, mask=None, connectivity=4, transform=None):
    """
    Extract polygon shapes from raster data.
    
    Parameters:
    - image (numpy.ndarray): Input raster array
    - mask (numpy.ndarray): Mask array (same shape as image)
    - connectivity (int): Pixel connectivity (4 or 8)
    - transform (Affine): Geospatial transformation for coordinates
    
    Yields:
    tuple: (geometry, value) pairs where geometry is GeoJSON-like dict
    """

from rasterio.features import shapes
import rasterio
import json

# Extract shapes from classified raster
with rasterio.open('landcover.tif') as dataset:
    # Read classification data
    landcover = dataset.read(1)
    
    # Extract all shapes
    shape_gen = shapes(
        landcover, 
        transform=dataset.transform,
        connectivity=8  # 8-connected pixels
    )
    
    # Convert to GeoJSON features
    features = []
    for geometry, value in shape_gen:
        feature = {
            'type': 'Feature',
            'geometry': geometry,
            'properties': {'class_id': int(value)}
        }
        features.append(feature)
    
    # Save as GeoJSON
    geojson = {
        'type': 'FeatureCollection',
        'features': features
    }
    
    with open('landcover_polygons.geojson', 'w') as f:
        json.dump(geojson, f)

# Extract shapes with custom mask
with rasterio.open('elevation.tif') as dataset:
    elevation = dataset.read(1)
    
    # Create mask for areas above threshold
    high_elevation_mask = elevation > 2000  # meters
    
    # Extract high elevation areas as shapes
    high_areas = list(shapes(
        elevation,
        mask=high_elevation_mask,
        transform=dataset.transform
    ))
    
    print(f"Found {len(high_areas)} high elevation areas")

def rasterize(shapes, out_shape=None, fill=0, out=None, transform=None, 
              all_touched=False, merge_alg=MergeAlg.replace, default_value=1, 
              dtype=None):
    """
    Burn vector geometries into raster array.
    
    Parameters:
    - shapes (iterable): Geometries with optional values [(geom, value), ...]
    - out_shape (tuple): Output array shape (height, width)
    - fill (number): Fill value for background pixels
    - out (numpy.ndarray): Pre-allocated output array
    - transform (Affine): Geospatial transformation
    - all_touched (bool): Burn all pixels touched by geometry
    - merge_alg (MergeAlg): Algorithm for merging values
    - default_value (number): Default value for geometries without values
    - dtype (numpy.dtype): Output data type
    
    Returns:
    numpy.ndarray: Rasterized array
    """

from rasterio.features import rasterize
from rasterio.enums import MergeAlg
import geopandas as gpd
import numpy as np

# Rasterize vector polygons
# Load vector data
gdf = gpd.read_file('study_areas.shp')

# Define raster parameters  
height, width = 1000, 1000
transform = from_bounds(*gdf.total_bounds, width, height)

# Prepare geometries with values
shapes_with_values = [
    (geom, value) for geom, value in 
    zip(gdf.geometry, gdf['area_id'])
]

# Rasterize polygons
rasterized = rasterize(
    shapes_with_values,
    out_shape=(height, width),
    transform=transform,
    fill=0,  # Background value
    all_touched=False,  # Only pixels with center inside geometry
    dtype='uint16'
)

# Save rasterized result
profile = {
    'driver': 'GTiff',
    'dtype': 'uint16',
    'nodata': 0,
    'width': width,
    'height': height,
    'count': 1,
    'crs': gdf.crs,
    'transform': transform
}

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

# Rasterize with different merge algorithms
# Create overlapping geometries
overlapping_shapes = [(geom, 1) for geom in gdf.geometry]

# Replace (default) - later values overwrite earlier ones
raster_replace = rasterize(overlapping_shapes, (height, width), transform=transform)

# Add - sum overlapping values
raster_add = rasterize(
    overlapping_shapes, 
    (height, width), 
    transform=transform,
    merge_alg=MergeAlg.add  
)

def geometry_mask(geometries, out_shape, transform, all_touched=False, 
                  invert=False):
    """
    Create boolean mask from geometries.
    
    Parameters:
    - geometries (iterable): Input geometries (GeoJSON-like)
    - out_shape (tuple): Output array shape (height, width)
    - transform (Affine): Geospatial transformation
    - all_touched (bool): Include pixels touched by geometries  
    - invert (bool): Invert mask (True inside geometries)
    
    Returns:
    numpy.ndarray: Boolean mask array (True = masked by default)
    """

from rasterio.features import geometry_mask

# Create mask from study area boundary
study_area = {
    'type': 'Polygon',
    'coordinates': [[
        [-120.5, 35.2], [-119.8, 35.2],
        [-119.8, 35.9], [-120.5, 35.9],
        [-120.5, 35.2]
    ]]
}

# Define raster grid
height, width = 500, 700
transform = from_bounds(-121, 35, -119, 36, width, height)

# Create mask (True outside polygon, False inside)
mask = geometry_mask(
    [study_area],
    out_shape=(height, width),
    transform=transform,
    invert=False  # False = not masked (inside geometry)
)

# Create inverted mask (True inside polygon)  
inside_mask = geometry_mask(
    [study_area],
    out_shape=(height, width),
    transform=transform,
    invert=True  # True = inside geometry
)

# Apply mask to raster data
with rasterio.open('satellite_image.tif') as dataset:
    # Read data
    data = dataset.read(1)
    
    # Apply study area mask
    masked_data = np.where(inside_mask, data, -9999)  # -9999 outside study area
    
    # Or use numpy masked arrays
    import numpy.ma as ma
    masked_array = ma.array(data, mask=~inside_mask)  # ~ inverts boolean

def sieve(image, size, out=None, mask=None, connectivity=4):
    """
    Remove small polygons from raster data.
    
    Parameters:
    - image (numpy.ndarray): Input raster array
    - size (int): Minimum polygon size (pixels)
    - out (numpy.ndarray): Output array
    - mask (numpy.ndarray): Mask array
    - connectivity (int): Pixel connectivity (4 or 8)
    
    Returns:
    numpy.ndarray: Sieved raster array
    """

def bounds(geometry):
    """
    Calculate bounding box of geometry.
    
    Parameters:
    - geometry (dict): GeoJSON-like geometry
    
    Returns:
    tuple: (left, bottom, right, top) bounds
    """

from rasterio.features import sieve, bounds

# Remove small polygons from classification
with rasterio.open('classification.tif') as dataset:
    classified = dataset.read(1)
    
    # Remove polygons smaller than 100 pixels
    cleaned = sieve(
        classified, 
        size=100,
        connectivity=8  # 8-connected neighborhoods
    )
    
    # Save cleaned classification
    profile = dataset.profile.copy()
    with rasterio.open('cleaned_classification.tif', 'w', **profile) as dst:
        dst.write(cleaned, 1)

# Calculate geometry bounds
polygon = {
    'type': 'Polygon', 
    'coordinates': [[
        [-120.5, 35.2], [-119.8, 35.2],
        [-119.8, 35.9], [-120.5, 35.9], 
        [-120.5, 35.2]
    ]]
}

geom_bounds = bounds(polygon)
print(f"Geometry bounds: {geom_bounds}")  # (-120.5, 35.2, -119.8, 35.9)

class MergeAlg(Enum):
    """Algorithms for merging rasterized values."""
    
    replace = 'replace'  # Later values overwrite earlier (default)
    add = 'add'          # Sum overlapping values

from rasterio.enums import MergeAlg

# Different merge strategies for overlapping features
shapes_with_weights = [(geometry, weight) for geometry, weight in data]

# Sum weights in overlapping areas
density_raster = rasterize(
    shapes_with_weights,
    out_shape=(height, width),
    transform=transform,
    merge_alg=MergeAlg.add,  # Sum overlapping values
    dtype='float32'
)

# Last feature wins in overlapping areas  
priority_raster = rasterize(
    shapes_with_weights,
    out_shape=(height, width), 
    transform=transform,
    merge_alg=MergeAlg.replace,  # Default behavior
    dtype='uint16'
)

# Integration with GeoPandas
import geopandas as gpd

# Read vector data
gdf = gpd.read_file('polygons.shp')

# Convert to rasterio format
geometries = gdf.geometry.apply(lambda x: x.__geo_interface__)
values = gdf['field_value']
shapes_list = list(zip(geometries, values))

# Rasterize
raster = rasterize(shapes_list, out_shape=(1000, 1000), transform=transform)

# Integration with Shapely
from shapely.geometry import Polygon, Point

# Create Shapely geometries
polygon = Polygon([(-120, 35), (-119, 35), (-119, 36), (-120, 36)])
point = Point(-119.5, 35.5)

# Convert to GeoJSON for rasterio
polygon_geojson = polygon.__geo_interface__
point_geojson = point.__geo_interface__

# Use in rasterio operations
mask = geometry_mask([polygon_geojson], (100, 100), transform)