Raster Processing

def reproject(source, destination, src_transform=None, src_crs=None, 
              src_nodata=None, dst_transform=None, dst_crs=None, 
              dst_nodata=None, resampling=Resampling.nearest, num_threads=1, 
              **kwargs):
    """
    Reproject raster data from source to destination CRS.
    
    Parameters:
    - source (numpy.ndarray or DatasetReader): Source raster data
    - destination (numpy.ndarray or DatasetWriter): Destination array or dataset
    - src_transform (Affine): Source geospatial transformation
    - src_crs (CRS): Source coordinate reference system
    - src_nodata (number): Source nodata value
    - dst_transform (Affine): Destination geospatial transformation  
    - dst_crs (CRS): Destination coordinate reference system
    - dst_nodata (number): Destination nodata value
    - resampling (Resampling): Resampling algorithm
    - num_threads (int): Number of processing threads
    - **kwargs: Additional GDAL warp options
    
    Returns:
    numpy.ndarray or None: Reprojected data (if destination is array)
    """

def calculate_default_transform(src_crs, dst_crs, width, height, left, bottom, 
                               right, top, resolution=None, dst_width=None, 
                               dst_height=None):
    """
    Calculate optimal transform and dimensions for reprojection.
    
    Parameters:
    - src_crs (CRS): Source coordinate reference system
    - dst_crs (CRS): Destination coordinate reference system
    - width (int): Source width in pixels
    - height (int): Source height in pixels
    - left, bottom, right, top (float): Source bounds
    - resolution (float): Target pixel resolution
    - dst_width, dst_height (int): Target dimensions
    
    Returns:
    tuple: (transform, width, height) for destination
    """

def aligned_target(transform, width, height, resolution):
    """
    Create aligned target parameters for reprojection.
    
    Parameters:
    - transform (Affine): Source transformation
    - width (int): Source width
    - height (int): Source height
    - resolution (float): Target resolution
    
    Returns:
    tuple: (aligned_transform, aligned_width, aligned_height)
    """

from rasterio.warp import reproject, calculate_default_transform, Resampling
from rasterio.crs import CRS
import numpy as np

# Reproject between datasets
with rasterio.open('input_wgs84.tif') as src:
    src_data = src.read(1)
    
    # Calculate target parameters
    dst_crs = CRS.from_epsg(3857)  # Web Mercator
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds
    )
    
    # Create destination array
    dst_data = np.empty((height, width), dtype=src_data.dtype)
    
    # Perform reprojection
    reproject(
        source=src_data,
        destination=dst_data,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.bilinear
    )

# Write reprojected data
profile = src.profile.copy()
profile.update({
    'crs': dst_crs,
    'transform': transform,
    'width': width,
    'height': height
})

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

def transform_geom(src_crs, dst_crs, geom, antimeridian_cutting=False, 
                   antimeridian_offset=10.0, precision=-1):
    """
    Transform geometry coordinates between CRS.
    
    Parameters:
    - src_crs (CRS): Source coordinate reference system
    - dst_crs (CRS): Destination coordinate reference system  
    - geom (dict): GeoJSON-like geometry
    - antimeridian_cutting (bool): Cut geometries at antimeridian
    - antimeridian_offset (float): Offset for antimeridian cutting
    - precision (int): Coordinate precision (decimal places)
    
    Returns:
    dict: Transformed geometry
    """

def transform_bounds(src_crs, dst_crs, left, bottom, right, top, densify_pts=21):
    """
    Transform bounding box coordinates between CRS.
    
    Parameters:
    - src_crs (CRS): Source coordinate reference system
    - dst_crs (CRS): Destination coordinate reference system
    - left, bottom, right, top (float): Source bounds
    - densify_pts (int): Number of points for accurate transformation
    
    Returns:
    tuple: (left, bottom, right, top) in destination CRS  
    """

from rasterio.warp import transform_geom, transform_bounds
from rasterio.crs import CRS

# Transform geometry
src_crs = CRS.from_epsg(4326)  # WGS84
dst_crs = CRS.from_epsg(3857)  # Web Mercator

# GeoJSON polygon
polygon = {
    'type': 'Polygon',
    'coordinates': [[
        [-120.0, 35.0], [-119.0, 35.0], 
        [-119.0, 36.0], [-120.0, 36.0], 
        [-120.0, 35.0]
    ]]
}

# Transform to Web Mercator
transformed_polygon = transform_geom(src_crs, dst_crs, polygon)

# Transform bounding box
src_bounds = (-120.5, 35.2, -119.8, 35.9)
dst_bounds = transform_bounds(src_crs, dst_crs, *src_bounds)
print(f"Transformed bounds: {dst_bounds}")

def mask(dataset, shapes, all_touched=False, invert=False, nodata=None, 
         filled=True, crop=False, pad=False, pad_width=0, indexes=None):
    """
    Mask raster using vector geometries.
    
    Parameters:
    - dataset (DatasetReader): Input raster dataset
    - shapes (iterable): Geometries for masking (GeoJSON-like)
    - all_touched (bool): Include pixels touched by geometries
    - invert (bool): Invert mask (exclude geometries instead)
    - nodata (number): NoData value for masked areas
    - filled (bool): Fill masked areas with nodata
    - crop (bool): Crop result to geometry bounds
    - pad (bool): Pad cropped result
    - pad_width (int): Padding width in pixels
    - indexes (sequence): Band indexes to mask
    
    Returns:
    tuple: (masked_array, output_transform)  
    """

def raster_geometry_mask(dataset, shapes, all_touched=False, invert=False, 
                        crop=False, pad=False, pad_width=0):
    """
    Create geometry mask for raster dataset.
    
    Parameters:
    - dataset (DatasetReader): Raster dataset
    - shapes (iterable): Masking geometries
    - all_touched (bool): Include touched pixels
    - invert (bool): Invert mask
    - crop (bool): Crop to geometry bounds  
    - pad (bool): Pad result
    - pad_width (int): Padding width
    
    Returns:
    tuple: (mask_array, output_transform)
    """

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

from rasterio.mask import mask, geometry_mask
import rasterio
import geopandas as gpd

# Mask raster with shapefile
with rasterio.open('landsat_scene.tif') as dataset:
    # Load vector mask
    shapes_gdf = gpd.read_file('study_area.shp')
    shapes = [geom.__geo_interface__ for geom in shapes_gdf.geometry]
    
    # Apply mask
    masked_data, masked_transform = mask(
        dataset, shapes, crop=True, filled=True, nodata=-9999
    )
    
    # Save masked result
    profile = dataset.profile.copy()
    profile.update({
        'height': masked_data.shape[1],
        'width': masked_data.shape[2], 
        'transform': masked_transform,
        'nodata': -9999
    })
    
    with rasterio.open('masked_landsat.tif', 'w', **profile) as dst:
        dst.write(masked_data)

# Create custom mask
polygon_coords = [(-120, 35), (-119, 35), (-119, 36), (-120, 36), (-120, 35)]
polygon = {'type': 'Polygon', 'coordinates': [polygon_coords]}

# Create boolean mask
mask_array = geometry_mask(
    [polygon], 
    out_shape=(1000, 1000),
    transform=from_bounds(-121, 34, -118, 37, 1000, 1000),
    invert=True  # True inside polygon
)

def merge(datasets, bounds=None, res=None, nodata=None, dtype=None, 
          precision=7, indexes=None, output_count=None, resampling=Resampling.nearest,
          method='first', **kwargs):
    """
    Merge multiple raster datasets into a mosaic.
    
    Parameters:
    - datasets (sequence): Input datasets (DatasetReader objects)
    - bounds (tuple): Output bounds (left, bottom, right, top) 
    - res (tuple): Output resolution (x_res, y_res)
    - nodata (number): Output nodata value
    - dtype (str): Output data type
    - precision (int): Coordinate precision
    - indexes (sequence): Band indexes to merge
    - output_count (int): Number of output bands
    - resampling (Resampling): Resampling algorithm
    - method (str): Merge method ('first', 'last', 'min', 'max', 'mean')
    - **kwargs: Additional options
    
    Returns:
    tuple: (merged_array, output_transform)
    """

from rasterio.merge import merge
import glob

# Merge multiple GeoTIFF files
tiff_files = glob.glob('tiles/*.tif')
datasets = [rasterio.open(f) for f in tiff_files]

try:
    # Merge with default settings (first method)
    mosaic, output_transform = merge(datasets)
    
    # Create output profile
    profile = datasets[0].profile.copy()
    profile.update({
        'height': mosaic.shape[1],
        'width': mosaic.shape[2],
        'transform': output_transform
    })
    
    # Save merged result
    with rasterio.open('merged_mosaic.tif', 'w', **profile) as dst:
        dst.write(mosaic)
        
    # Merge with specific bounds and resolution
    bounds = (-120, 35, -119, 36)  # Custom extent  
    res = (0.001, 0.001)          # 0.001 degree pixels
    
    mosaic_custom, transform_custom = merge(
        datasets, 
        bounds=bounds,
        res=res,
        method='mean',  # Average overlapping pixels
        resampling=Resampling.bilinear
    )

finally:
    # Close all datasets
    for dataset in datasets:
        dataset.close()

class Resampling(Enum):
    """Resampling algorithms enumeration."""
    
    nearest = 'nearest'           # Nearest neighbor (fastest, categorical data)
    bilinear = 'bilinear'         # Bilinear interpolation (continuous data)
    cubic = 'cubic'               # Cubic convolution (smooth continuous data)
    cubic_spline = 'cubic_spline' # Cubic spline (very smooth)
    lanczos = 'lanczos'           # Lanczos windowed sinc (sharp details)
    average = 'average'           # Average of contributing pixels
    mode = 'mode'                 # Most common value (categorical data)
    gauss = 'gauss'               # Gaussian kernel
    max = 'max'                   # Maximum value
    min = 'min'                   # Minimum value
    med = 'med'                   # Median value
    q1 = 'q1'                     # First quartile
    q3 = 'q3'                     # Third quartile
    sum = 'sum'                   # Sum of values
    rms = 'rms'                   # Root mean square

from rasterio.enums import Resampling

# Resampling algorithm selection guide:

# Categorical data (land cover, classifications)
resampling = Resampling.nearest  # Preserves original values
resampling = Resampling.mode     # Most common value in area

# Continuous data (elevation, temperature)  
resampling = Resampling.bilinear  # Good balance of speed and quality
resampling = Resampling.cubic     # Smoother interpolation
resampling = Resampling.lanczos   # Sharp details, may create artifacts

# Statistical summaries
resampling = Resampling.average   # Mean value
resampling = Resampling.max       # Maximum value  
resampling = Resampling.min       # Minimum value
resampling = Resampling.med       # Median value

# Use in reprojection
reproject(
    source=src_data,
    destination=dst_data,
    src_transform=src_transform,
    dst_transform=dst_transform,
    src_crs=src_crs,
    dst_crs=dst_crs,
    resampling=Resampling.cubic  # High-quality interpolation
)