rasterio

# Rasterio — Raster Geospatial Processing Rasterio is the standard Python library for reading and writing georeferenced raster data. It wraps GDAL but exposes a clean Pythonic API. Every raster has two components: **pixel values** (a numpy array) and **geospatial metadata** (CRS, transform, bounds) that maps pixels to real-world coordinates. ## Raster vs Vector — When to Use What ``` RASTER (Rasterio) VECTOR (GeoPandas / Shapely) ──────────────────────────────── ──────────────────────────────── Grid of pixels Points, lines, polygons Satellite imagery Administrative boundaries Elevation models (DEM) Road networks Land cover maps Locations / GPS tracks Temperature grids Census regions Spectral data Feature geometries KEY RULE: When your data IS a grid → Rasterio. When your data IS shapes/points → GeoPandas. When you need BOTH → Rasterio clips rasters TO vector shapes. ``` ## Reference Documentation **Rasterio docs**: https://rasterio.readthedocs.io/en/latest/ **GDAL formats**: https://gdal.org/drivers/raster/index.html **GitHub**: https://github.com/rasterio/rasterio **Search patterns**: `rasterio.open`, `dataset.read`, `dataset.transform`, `rasterio.mask`, `show` ## Core Principles ### The Dataset Object `rasterio.open()` returns a dataset. A dataset has **bands** (layers of pixel data), a **transform** (maps pixel coordinates to geographic coordinates), a **CRS** (coordinate reference system), and **bounds** (geographic extent). Always use context manager (`with` statement) — it handles file handles correctly. ### The Transform An affine transform maps pixel (col, row) to geographic (x, y). `dataset.transform` gives this mapping. `dataset.index(x, y)` gives the reverse: geographic → pixel. This is how you go from "latitude/longitude" to "which pixel?" ### Bands Most satellite imagery is multi-band: Band 1 = Red, Band 2 = Green, Band 3 = Blue, Band 4 = NIR, etc. Band numbering starts at **1** (not 0). `dataset.read(1)` reads band 1 as a 2D numpy array. ### CRS — Coordinate Reference Systems Every raster is projected into some CRS. `EPSG:4326` = WGS84 lat/lon (GPS coordinates). `EPSG:32633` = UTM zone 33N (meters, good for Europe). Operations between rasters in different CRS require reprojection first. ### NoData Pixels outside valid coverage are marked with a `nodata` value (e.g., -9999, 0, or NaN). Always mask these before computation — including them corrupts statistics. ## Quick Reference ### Installation ```bash pip install rasterio numpy matplotlib # For vector integration: pip install geopandas shapely fiona ``` ### Standard Imports ```python import rasterio from rasterio.transform import from_bounds, Affine from rasterio.crs import CRS import numpy as np import matplotlib.pyplot as plt ``` ### Basic Pattern — Read, Inspect, Visualize ```python import rasterio import numpy as np import matplotlib.pyplot as plt with rasterio.open('image.tif') as src: # Metadata print(f"Bands: {src.count}") print(f"Shape: {src.height} x {src.width}") print(f"CRS: {src.crs}") print(f"Bounds: {src.bounds}") # (left, bottom, right, top) print(f"Resolution: {src.res}") # (x_res, y_res) in CRS units print(f"NoData: {src.nodata}") print(f"Transform: {src.transform}") # Read all bands → shape: (bands, height, width) data = src.read() # Read single band → shape: (height, width) band1 = src.read(1) # Read with masked nodata → numpy masked array band1_masked = src.read(1, masked=True) # Visualize single band plt.imshow(band1, cmap='viridis') plt.colorbar(label='Value') plt.title('Band 1') plt.show() ``` ### Basic Pattern — Write a Raster ```python import rasterio from rasterio.transform import from_bounds import numpy as np # Create a 100x100 single-band raster covering a geographic extent height, width = 100, 100 data = np.random.rand(height, width).astype(np.float32) transform = from_bounds( west=10.0, south=50.0, east=11.0, north=51.0, width=width, height=height ) with rasterio.open( 'output.tif', 'w', driver='GTiff', height=height, width=width, count=1, # Number of bands dtype=data.dtype, crs='EPSG:4326', # WGS84 transform=transform, nodata=-9999 ) as dst: dst.write(data, 1) # Write to band 1 ``` ## Critical Rules ### ✅ DO - **Always use `with rasterio.open(...)` context manager** — Ensures file handles are closed properly. Never do `src = rasterio.open(...)` without `with`. - **Use `masked=True` when reading** — Returns a numpy masked array that automatically excludes nodata pixels. Prevents nodata from corrupting calculations. - **Match CRS before any spatial operation** — If two rasters have different CRS, reproject one before combining. Use `rasterio.warp.reproject()`. - **Use windowed reading for large files** — Files > 1GB will OOM if read entirely. Use `src.read(window=...)` to read chunks. - **Preserve geospatial metadata when writing** — Copy `transform`, `crs`, `nodata` from the source when creating derived rasters. - **Use `float32` for computed indices** — NDVI, NDWI produce float values [-1, 1]. Input bands are often `uint8` or `uint16` — cast before division to avoid integer truncation. - **Check `src.nodata` before computing** — If nodata is None, the file has no declared nodata value. Handle accordingly. - **Use `rasterio.features` for vector-raster conversion** — Don't roll your own rasterization. ### ❌ DON'T - **Don't mix up band indexing** — Rasterio bands start at **1**. `src.read(1)` = first band. numpy arrays are 0-indexed, so `data[0]` = first band after `src.read()`. - **Don't forget to cast dtypes before band math** — `uint8 / uint8` = integer division = truncated to 0. Always cast: `band.astype(np.float32)`. - **Don't assume all rasters share the same CRS** — Even "standard" datasets may differ. Always check and align. - **Don't ignore resolution differences** — Two rasters covering the same area may have different pixel sizes. Resample to match before pixel-wise operations. - **Don't hardcode nodata as 0** — Nodata values vary by file. Always read `src.nodata`. - **Don't read entire multi-GB files into memory** — Use windowed I/O or overview levels. ## Anti-Patterns (NEVER) ```python import rasterio import numpy as np # ❌ BAD: Integer division in band math — truncates to 0 with rasterio.open('sentinel.tif') as src: nir = src.read(4) # uint16 red = src.read(3) # uint16 ndvi = (nir - red) / (nir + red) # INTEGER DIVISION → all zeros or ±1! # ✅ GOOD: Cast to float BEFORE any arithmetic with rasterio.open('sentinel.tif') as src: nir = src.read(4).astype(np.float32) red = src.read(3).astype(np.float32) ndvi = (nir - red) / (nir + red + 1e-10) # +1e-10 avoids div-by-zero # ───────────────────────────────────────────────────────────── # ❌ BAD: No context manager — file handle leak src = rasterio.open('image.tif') data = src.read(1) # src never closed → resource leak, potential corruption on write # ✅ GOOD: Context manager with rasterio.open('image.tif') as src: data = src.read(1) # File closed automatically # ───────────────────────────────────────────────────────────── # ❌ BAD: Ignoring nodata in statistics with rasterio.open('dem.tif') as src: elev = src.read(1) mean_elevation = elev.mean() # Includes -9999 nodata → wrong answer! # ✅ GOOD: Mask nodata before computing with rasterio.open('dem.tif') as src: elev = src.read(1, masked=True) # Masked array mean_elevation = float(elev.mean()) # Ignores masked (nodata) pixels # OR explicitly: valid = elev[elev != src.nodata] mean_elevation = valid.mean() # ───────────────────────────────────────────────────────────── # ❌ BAD: Operating on rasters with different CRS without reprojection with rasterio.open('landcover_utm.tif') as src1: with rasterio.open('population_wgs84.tif') as src2: # src1.crs = EPSG:32633 (UTM), src2.crs = EPSG:4326 (WGS84) pop = src2.read(1) land = src1.read(1) result = pop * land # WRONG — pixels don't align geographically! # ✅ GOOD: Reproject to common CRS first (see Reprojection section) ``` ## Reading Rasters ### Full Read vs Selective Read ```python import rasterio import numpy as np with rasterio.open('multispectral.tif') as src: # All bands at once → shape: (n_bands, height, width) all_bands = src.read() # Single band → shape: (height, width) band1 = src.read(1) # Multiple specific bands → shape: (n_selected, height, width) rgb = src.read([1, 2, 3]) # With nodata masking band1_masked = src.read(1, masked=True) # numpy.ma.MaskedArray # Specific spatial subset (window) from rasterio.windows import Window window = Window(col_off=100, row_off=200, width=50, height=50) patch = src.read(1, window=window) # shape: (50, 50) # Read at lower resolution (overview) # out_shape forces resampling to smaller array small = src.read(1, out_shape=(src.height // 4, src.width // 4)) ``` ### Coordinate ↔ Pixel Mapping ```python import rasterio with rasterio.open('image.tif') as src: # Geographic coordinates → pixel (row, col) row, col = src.index(lon, lat) # e.g., src.index(10.5, 50.3) # Pixel (row, col) → geographic coordinates (x, y) x, y = src.xy(row, col) # Center of that pixel # Geographic extent of a specific pixel # Top-left corner of pixel (row, col): x_tl = src.transform.c + col * src.transform.a y_tl = src.transform.f + row * src.transform.e # Bounds of the entire raster print(f"Left={src.bounds.left}, Bottom={src.bounds.bottom}, " f"Right={src.bounds.right}, Top={src.bounds.top}") ``` ## Writing Rasters ### Write with Full Metadata ```python import rasterio from rasterio.transform import from_bounds, Affine import numpy as np def write_raster(data: np.ndarray, output_path: str, transform: Affine, crs: str = 'EPSG:4326', nodata: float = -9999, band_descriptions: list[str] = None): """ Write a numpy array as a georeferenced GeoTIFF. data shape: (bands, height, width) or (height, width) for single band. """ if data.ndim == 2: data = data[np.newaxis, :] # Add band dimension → (1, H, W) n_bands, height, width = data.shape with rasterio.open( output_path, 'w', driver='GTiff', height=height, width=width, count=n_bands, dtype=data.dtype, crs=crs, transform=transform, nodata=nodata, compress='lzw', # Compression — reduces file size significantly tiled=True, # Tiled layout — better for large file random access blockxsize=256, blockysize=256 ) as dst: dst.write(data) if band_descriptions: for i, desc in enumerate(band_descriptions, 1): dst.set_band_description(i, desc) # ─── Copy metadata pattern (most common: derive new raster from existing) ─── with rasterio.open('input.tif') as src: computed = src.read(1).astype(np.float32) * 2.0 # Some computation # Copy ALL metadata from source, override only what changed profile = src.profile.copy() profile.update(dtype=computed.dtype, count=1, nodata=-9999) with rasterio.open('output.tif', 'w', **profile) as dst: dst.write(computed, 1) ``` ## CRS and Reprojection ### Check and Compare CRS ```python import rasterio from rasterio.crs import CRS with rasterio.open('image.tif') as src: print(f"CRS: {src.crs}") print(f"EPSG: {src.crs.to_epsg()}") print(f"WKT: {src.crs.to_wkt()}") # Check if two CRS are the same target_crs = CRS.from_epsg(32633) # UTM zone 33N print(f"Same as UTM 33N? {src.crs == target_crs}") # Common EPSG codes: # EPSG:4326 — WGS84 lat/lon (GPS, most web maps) # EPSG:3857 — Web Mercator (Google Maps, OpenStreetMap tiles) # EPSG:32601–32660 — UTM zones 1–60 North (meters, good for local analysis) # EPSG:32701–32760 — UTM zones 1–60 South ``` ### Reproject a Raster ```python import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling from rasterio.crs import CRS import numpy as np def reproject_raster(input_path: str, output_path: str, target_crs: str = 'EPSG:4326'): """Reproject an entire raster file to a new CRS.""" target_crs = CRS.from_user_input(target_crs) with rasterio.open(input_path) as src: # Calculate the optimal transform and dimensions for the target CRS transform, width, height = calculate_default_transform( src.crs, target_crs, src.width, src.height, *src.bounds ) # Prepare output profile profile = src.profile.copy() profile.update(crs=target_crs, transform=transform, width=width, height=height) with rasterio.open(output_path, 'w', **profile) as dst: for band_idx in range(1, src.count + 1): reproject( source=rasterio.band(src, band_idx), destination=rasterio.band(dst, band_idx), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=target_crs, resampling=Resampling.bilinear # bilinear for continuous data # Use Resampling.nearest for categorical (land cover class IDs) ) # reproject_raster('utm_image.tif', 'wgs84_image.tif', 'EPSG:4326') ``` ## Band Math — Spectral Indices ```python import rasterio import numpy as np def compute_indices(input_path: str, output_path: str, band_map: dict = None) -> dict: