0
# Raster Processing
1

2
Advanced processing operations including reprojection, masking, merging, and resampling. These functions provide comprehensive raster manipulation capabilities with support for various algorithms and coordinate systems.
3

4
## Capabilities
5

6
### Reprojection and Warping
7

8
Transform raster data between different coordinate reference systems with various resampling algorithms.
9

10
```python { .api }
11
def reproject(source, destination, src_transform=None, src_crs=None, 
12
              src_nodata=None, dst_transform=None, dst_crs=None, 
13
              dst_nodata=None, resampling=Resampling.nearest, num_threads=1, 
14
              **kwargs):
15
    """
16
    Reproject raster data from source to destination CRS.
17
    
18
    Parameters:
19
    - source (numpy.ndarray or DatasetReader): Source raster data
20
    - destination (numpy.ndarray or DatasetWriter): Destination array or dataset
21
    - src_transform (Affine): Source geospatial transformation
22
    - src_crs (CRS): Source coordinate reference system
23
    - src_nodata (number): Source nodata value
24
    - dst_transform (Affine): Destination geospatial transformation  
25
    - dst_crs (CRS): Destination coordinate reference system
26
    - dst_nodata (number): Destination nodata value
27
    - resampling (Resampling): Resampling algorithm
28
    - num_threads (int): Number of processing threads
29
    - **kwargs: Additional GDAL warp options
30
    
31
    Returns:
32
    numpy.ndarray or None: Reprojected data (if destination is array)
33
    """
34

35
def calculate_default_transform(src_crs, dst_crs, width, height, left, bottom, 
36
                               right, top, resolution=None, dst_width=None, 
37
                               dst_height=None):
38
    """
39
    Calculate optimal transform and dimensions for reprojection.
40
    
41
    Parameters:
42
    - src_crs (CRS): Source coordinate reference system
43
    - dst_crs (CRS): Destination coordinate reference system
44
    - width (int): Source width in pixels
45
    - height (int): Source height in pixels
46
    - left, bottom, right, top (float): Source bounds
47
    - resolution (float): Target pixel resolution
48
    - dst_width, dst_height (int): Target dimensions
49
    
50
    Returns:
51
    tuple: (transform, width, height) for destination
52
    """
53

54
def aligned_target(transform, width, height, resolution):
55
    """
56
    Create aligned target parameters for reprojection.
57
    
58
    Parameters:
59
    - transform (Affine): Source transformation
60
    - width (int): Source width
61
    - height (int): Source height
62
    - resolution (float): Target resolution
63
    
64
    Returns:
65
    tuple: (aligned_transform, aligned_width, aligned_height)
66
    """
67
```
68

69
Usage examples:
70

71
```python
72
from rasterio.warp import reproject, calculate_default_transform, Resampling
73
from rasterio.crs import CRS
74
import numpy as np
75

76
# Reproject between datasets
77
with rasterio.open('input_wgs84.tif') as src:
78
    src_data = src.read(1)
79
    
80
    # Calculate target parameters
81
    dst_crs = CRS.from_epsg(3857)  # Web Mercator
82
    transform, width, height = calculate_default_transform(
83
        src.crs, dst_crs, src.width, src.height, *src.bounds
84
    )
85
    
86
    # Create destination array
87
    dst_data = np.empty((height, width), dtype=src_data.dtype)
88
    
89
    # Perform reprojection
90
    reproject(
91
        source=src_data,
92
        destination=dst_data,
93
        src_transform=src.transform,
94
        src_crs=src.crs,
95
        dst_transform=transform,
96
        dst_crs=dst_crs,
97
        resampling=Resampling.bilinear
98
    )
99

100
# Write reprojected data
101
profile = src.profile.copy()
102
profile.update({
103
    'crs': dst_crs,
104
    'transform': transform,
105
    'width': width,
106
    'height': height
107
})
108

109
with rasterio.open('output_mercator.tif', 'w', **profile) as dst:
110
    dst.write(dst_data, 1)
111
```
112

113
### Coordinate Transformations
114

115
Transform coordinates and geometries between different coordinate reference systems.
116

117
```python { .api }
118
def transform_geom(src_crs, dst_crs, geom, antimeridian_cutting=False, 
119
                   antimeridian_offset=10.0, precision=-1):
120
    """
121
    Transform geometry coordinates between CRS.
122
    
123
    Parameters:
124
    - src_crs (CRS): Source coordinate reference system
125
    - dst_crs (CRS): Destination coordinate reference system  
126
    - geom (dict): GeoJSON-like geometry
127
    - antimeridian_cutting (bool): Cut geometries at antimeridian
128
    - antimeridian_offset (float): Offset for antimeridian cutting
129
    - precision (int): Coordinate precision (decimal places)
130
    
131
    Returns:
132
    dict: Transformed geometry
133
    """
134

135
def transform_bounds(src_crs, dst_crs, left, bottom, right, top, densify_pts=21):
136
    """
137
    Transform bounding box coordinates between CRS.
138
    
139
    Parameters:
140
    - src_crs (CRS): Source coordinate reference system
141
    - dst_crs (CRS): Destination coordinate reference system
142
    - left, bottom, right, top (float): Source bounds
143
    - densify_pts (int): Number of points for accurate transformation
144
    
145
    Returns:
146
    tuple: (left, bottom, right, top) in destination CRS  
147
    """
148
```
149

150
Usage examples:
151

152
```python
153
from rasterio.warp import transform_geom, transform_bounds
154
from rasterio.crs import CRS
155

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

160
# GeoJSON polygon
161
polygon = {
162
    'type': 'Polygon',
163
    'coordinates': [[
164
        [-120.0, 35.0], [-119.0, 35.0], 
165
        [-119.0, 36.0], [-120.0, 36.0], 
166
        [-120.0, 35.0]
167
    ]]
168
}
169

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

173
# Transform bounding box
174
src_bounds = (-120.5, 35.2, -119.8, 35.9)
175
dst_bounds = transform_bounds(src_crs, dst_crs, *src_bounds)
176
print(f"Transformed bounds: {dst_bounds}")
177
```
178

179
### Masking Operations
180

181
Apply vector masks to raster data for spatial subsetting and analysis.
182

183
```python { .api }
184
def mask(dataset, shapes, all_touched=False, invert=False, nodata=None, 
185
         filled=True, crop=False, pad=False, pad_width=0, indexes=None):
186
    """
187
    Mask raster using vector geometries.
188
    
189
    Parameters:
190
    - dataset (DatasetReader): Input raster dataset
191
    - shapes (iterable): Geometries for masking (GeoJSON-like)
192
    - all_touched (bool): Include pixels touched by geometries
193
    - invert (bool): Invert mask (exclude geometries instead)
194
    - nodata (number): NoData value for masked areas
195
    - filled (bool): Fill masked areas with nodata
196
    - crop (bool): Crop result to geometry bounds
197
    - pad (bool): Pad cropped result
198
    - pad_width (int): Padding width in pixels
199
    - indexes (sequence): Band indexes to mask
200
    
201
    Returns:
202
    tuple: (masked_array, output_transform)  
203
    """
204

205
def raster_geometry_mask(dataset, shapes, all_touched=False, invert=False, 
206
                        crop=False, pad=False, pad_width=0):
207
    """
208
    Create geometry mask for raster dataset.
209
    
210
    Parameters:
211
    - dataset (DatasetReader): Raster dataset
212
    - shapes (iterable): Masking geometries
213
    - all_touched (bool): Include touched pixels
214
    - invert (bool): Invert mask
215
    - crop (bool): Crop to geometry bounds  
216
    - pad (bool): Pad result
217
    - pad_width (int): Padding width
218
    
219
    Returns:
220
    tuple: (mask_array, output_transform)
221
    """
222

223
def geometry_mask(geometries, out_shape, transform, all_touched=False, 
224
                  invert=False):
225
    """
226
    Create boolean mask from geometries.
227
    
228
    Parameters:
229
    - geometries (iterable): Input geometries
230
    - out_shape (tuple): Output array shape (height, width)
231
    - transform (Affine): Geospatial transformation
232
    - all_touched (bool): Include touched pixels
233
    - invert (bool): Invert mask
234
    
235
    Returns:
236
    numpy.ndarray: Boolean mask array
237
    """
238
```
239

240
Usage examples:
241

242
```python
243
from rasterio.mask import mask, geometry_mask
244
import rasterio
245
import geopandas as gpd
246

247
# Mask raster with shapefile
248
with rasterio.open('landsat_scene.tif') as dataset:
249
    # Load vector mask
250
    shapes_gdf = gpd.read_file('study_area.shp')
251
    shapes = [geom.__geo_interface__ for geom in shapes_gdf.geometry]
252
    
253
    # Apply mask
254
    masked_data, masked_transform = mask(
255
        dataset, shapes, crop=True, filled=True, nodata=-9999
256
    )
257
    
258
    # Save masked result
259
    profile = dataset.profile.copy()
260
    profile.update({
261
        'height': masked_data.shape[1],
262
        'width': masked_data.shape[2], 
263
        'transform': masked_transform,
264
        'nodata': -9999
265
    })
266
    
267
    with rasterio.open('masked_landsat.tif', 'w', **profile) as dst:
268
        dst.write(masked_data)
269

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

274
# Create boolean mask
275
mask_array = geometry_mask(
276
    [polygon], 
277
    out_shape=(1000, 1000),
278
    transform=from_bounds(-121, 34, -118, 37, 1000, 1000),
279
    invert=True  # True inside polygon
280
)
281
```
282

283
### Merging Operations
284

285
Combine multiple raster datasets into a single mosaic with various compositing methods.
286

287
```python { .api }
288
def merge(datasets, bounds=None, res=None, nodata=None, dtype=None, 
289
          precision=7, indexes=None, output_count=None, resampling=Resampling.nearest,
290
          method='first', **kwargs):
291
    """
292
    Merge multiple raster datasets into a mosaic.
293
    
294
    Parameters:
295
    - datasets (sequence): Input datasets (DatasetReader objects)
296
    - bounds (tuple): Output bounds (left, bottom, right, top) 
297
    - res (tuple): Output resolution (x_res, y_res)
298
    - nodata (number): Output nodata value
299
    - dtype (str): Output data type
300
    - precision (int): Coordinate precision
301
    - indexes (sequence): Band indexes to merge
302
    - output_count (int): Number of output bands
303
    - resampling (Resampling): Resampling algorithm
304
    - method (str): Merge method ('first', 'last', 'min', 'max', 'mean')
305
    - **kwargs: Additional options
306
    
307
    Returns:
308
    tuple: (merged_array, output_transform)
309
    """
310
```
311

312
Usage examples:
313

314
```python
315
from rasterio.merge import merge
316
import glob
317

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

322
try:
323
    # Merge with default settings (first method)
324
    mosaic, output_transform = merge(datasets)
325
    
326
    # Create output profile
327
    profile = datasets[0].profile.copy()
328
    profile.update({
329
        'height': mosaic.shape[1],
330
        'width': mosaic.shape[2],
331
        'transform': output_transform
332
    })
333
    
334
    # Save merged result
335
    with rasterio.open('merged_mosaic.tif', 'w', **profile) as dst:
336
        dst.write(mosaic)
337
        
338
    # Merge with specific bounds and resolution
339
    bounds = (-120, 35, -119, 36)  # Custom extent  
340
    res = (0.001, 0.001)          # 0.001 degree pixels
341
    
342
    mosaic_custom, transform_custom = merge(
343
        datasets, 
344
        bounds=bounds,
345
        res=res,
346
        method='mean',  # Average overlapping pixels
347
        resampling=Resampling.bilinear
348
    )
349

350
finally:
351
    # Close all datasets
352
    for dataset in datasets:
353
        dataset.close()
354
```
355

356
### Resampling Algorithms
357

358
Rasterio provides various resampling algorithms for different use cases:
359

360
```python { .api }
361
class Resampling(Enum):
362
    """Resampling algorithms enumeration."""
363
    
364
    nearest = 'nearest'           # Nearest neighbor (fastest, categorical data)
365
    bilinear = 'bilinear'         # Bilinear interpolation (continuous data)
366
    cubic = 'cubic'               # Cubic convolution (smooth continuous data)
367
    cubic_spline = 'cubic_spline' # Cubic spline (very smooth)
368
    lanczos = 'lanczos'           # Lanczos windowed sinc (sharp details)
369
    average = 'average'           # Average of contributing pixels
370
    mode = 'mode'                 # Most common value (categorical data)
371
    gauss = 'gauss'               # Gaussian kernel
372
    max = 'max'                   # Maximum value
373
    min = 'min'                   # Minimum value
374
    med = 'med'                   # Median value
375
    q1 = 'q1'                     # First quartile
376
    q3 = 'q3'                     # Third quartile
377
    sum = 'sum'                   # Sum of values
378
    rms = 'rms'                   # Root mean square
379
```
380

381
Usage guidelines:
382

383
```python
384
from rasterio.enums import Resampling
385

386
# Resampling algorithm selection guide:
387

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

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

397
# Statistical summaries
398
resampling = Resampling.average   # Mean value
399
resampling = Resampling.max       # Maximum value  
400
resampling = Resampling.min       # Minimum value
401
resampling = Resampling.med       # Median value
402

403
# Use in reprojection
404
reproject(
405
    source=src_data,
406
    destination=dst_data,
407
    src_transform=src_transform,
408
    dst_transform=dst_transform,
409
    src_crs=src_crs,
410
    dst_crs=dst_crs,
411
    resampling=Resampling.cubic  # High-quality interpolation
412
)
413
```