0
# Feature Operations
1

2
Conversion between raster and vector data including shape extraction, geometry rasterization, and spatial analysis operations. These functions bridge raster and vector geospatial data formats.
3

4
## Capabilities
5

6
### Shape Extraction
7

8
Extract polygon shapes from raster data, converting raster regions into vector geometries.
9

10
```python { .api }
11
def shapes(image, mask=None, connectivity=4, transform=None):
12
    """
13
    Extract polygon shapes from raster data.
14
    
15
    Parameters:
16
    - image (numpy.ndarray): Input raster array
17
    - mask (numpy.ndarray): Mask array (same shape as image)
18
    - connectivity (int): Pixel connectivity (4 or 8)
19
    - transform (Affine): Geospatial transformation for coordinates
20
    
21
    Yields:
22
    tuple: (geometry, value) pairs where geometry is GeoJSON-like dict
23
    """
24
```
25

26
Usage examples:
27

28
```python
29
from rasterio.features import shapes
30
import rasterio
31
import json
32

33
# Extract shapes from classified raster
34
with rasterio.open('landcover.tif') as dataset:
35
    # Read classification data
36
    landcover = dataset.read(1)
37
    
38
    # Extract all shapes
39
    shape_gen = shapes(
40
        landcover, 
41
        transform=dataset.transform,
42
        connectivity=8  # 8-connected pixels
43
    )
44
    
45
    # Convert to GeoJSON features
46
    features = []
47
    for geometry, value in shape_gen:
48
        feature = {
49
            'type': 'Feature',
50
            'geometry': geometry,
51
            'properties': {'class_id': int(value)}
52
        }
53
        features.append(feature)
54
    
55
    # Save as GeoJSON
56
    geojson = {
57
        'type': 'FeatureCollection',
58
        'features': features
59
    }
60
    
61
    with open('landcover_polygons.geojson', 'w') as f:
62
        json.dump(geojson, f)
63

64
# Extract shapes with custom mask
65
with rasterio.open('elevation.tif') as dataset:
66
    elevation = dataset.read(1)
67
    
68
    # Create mask for areas above threshold
69
    high_elevation_mask = elevation > 2000  # meters
70
    
71
    # Extract high elevation areas as shapes
72
    high_areas = list(shapes(
73
        elevation,
74
        mask=high_elevation_mask,
75
        transform=dataset.transform
76
    ))
77
    
78
    print(f"Found {len(high_areas)} high elevation areas")
79
```
80

81
### Geometry Rasterization
82

83
Convert vector geometries into raster format by "burning" them into pixel arrays.
84

85
```python { .api }
86
def rasterize(shapes, out_shape=None, fill=0, out=None, transform=None, 
87
              all_touched=False, merge_alg=MergeAlg.replace, default_value=1, 
88
              dtype=None):
89
    """
90
    Burn vector geometries into raster array.
91
    
92
    Parameters:
93
    - shapes (iterable): Geometries with optional values [(geom, value), ...]
94
    - out_shape (tuple): Output array shape (height, width)
95
    - fill (number): Fill value for background pixels
96
    - out (numpy.ndarray): Pre-allocated output array
97
    - transform (Affine): Geospatial transformation
98
    - all_touched (bool): Burn all pixels touched by geometry
99
    - merge_alg (MergeAlg): Algorithm for merging values
100
    - default_value (number): Default value for geometries without values
101
    - dtype (numpy.dtype): Output data type
102
    
103
    Returns:
104
    numpy.ndarray: Rasterized array
105
    """
106
```
107

108
Usage examples:
109

110
```python
111
from rasterio.features import rasterize
112
from rasterio.enums import MergeAlg
113
import geopandas as gpd
114
import numpy as np
115

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

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

124
# Prepare geometries with values
125
shapes_with_values = [
126
    (geom, value) for geom, value in 
127
    zip(gdf.geometry, gdf['area_id'])
128
]
129

130
# Rasterize polygons
131
rasterized = rasterize(
132
    shapes_with_values,
133
    out_shape=(height, width),
134
    transform=transform,
135
    fill=0,  # Background value
136
    all_touched=False,  # Only pixels with center inside geometry
137
    dtype='uint16'
138
)
139

140
# Save rasterized result
141
profile = {
142
    'driver': 'GTiff',
143
    'dtype': 'uint16',
144
    'nodata': 0,
145
    'width': width,
146
    'height': height,
147
    'count': 1,
148
    'crs': gdf.crs,
149
    'transform': transform
150
}
151

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

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

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

162
# Add - sum overlapping values
163
raster_add = rasterize(
164
    overlapping_shapes, 
165
    (height, width), 
166
    transform=transform,
167
    merge_alg=MergeAlg.add  
168
)
169
```
170

171
### Geometry Masking
172

173
Create boolean masks from vector geometries for spatial filtering and analysis.
174

175
```python { .api }
176
def geometry_mask(geometries, out_shape, transform, all_touched=False, 
177
                  invert=False):
178
    """
179
    Create boolean mask from geometries.
180
    
181
    Parameters:
182
    - geometries (iterable): Input geometries (GeoJSON-like)
183
    - out_shape (tuple): Output array shape (height, width)
184
    - transform (Affine): Geospatial transformation
185
    - all_touched (bool): Include pixels touched by geometries  
186
    - invert (bool): Invert mask (True inside geometries)
187
    
188
    Returns:
189
    numpy.ndarray: Boolean mask array (True = masked by default)
190
    """
191
```
192

193
Usage examples:
194

195
```python
196
from rasterio.features import geometry_mask
197

198
# Create mask from study area boundary
199
study_area = {
200
    'type': 'Polygon',
201
    'coordinates': [[
202
        [-120.5, 35.2], [-119.8, 35.2],
203
        [-119.8, 35.9], [-120.5, 35.9],
204
        [-120.5, 35.2]
205
    ]]
206
}
207

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

212
# Create mask (True outside polygon, False inside)
213
mask = geometry_mask(
214
    [study_area],
215
    out_shape=(height, width),
216
    transform=transform,
217
    invert=False  # False = not masked (inside geometry)
218
)
219

220
# Create inverted mask (True inside polygon)  
221
inside_mask = geometry_mask(
222
    [study_area],
223
    out_shape=(height, width),
224
    transform=transform,
225
    invert=True  # True = inside geometry
226
)
227

228
# Apply mask to raster data
229
with rasterio.open('satellite_image.tif') as dataset:
230
    # Read data
231
    data = dataset.read(1)
232
    
233
    # Apply study area mask
234
    masked_data = np.where(inside_mask, data, -9999)  # -9999 outside study area
235
    
236
    # Or use numpy masked arrays
237
    import numpy.ma as ma
238
    masked_array = ma.array(data, mask=~inside_mask)  # ~ inverts boolean
239
```
240

241
### Spatial Analysis Utilities
242

243
Additional utilities for spatial analysis and feature processing.
244

245
```python { .api }
246
def sieve(image, size, out=None, mask=None, connectivity=4):
247
    """
248
    Remove small polygons from raster data.
249
    
250
    Parameters:
251
    - image (numpy.ndarray): Input raster array
252
    - size (int): Minimum polygon size (pixels)
253
    - out (numpy.ndarray): Output array
254
    - mask (numpy.ndarray): Mask array
255
    - connectivity (int): Pixel connectivity (4 or 8)
256
    
257
    Returns:
258
    numpy.ndarray: Sieved raster array
259
    """
260

261
def bounds(geometry):
262
    """
263
    Calculate bounding box of geometry.
264
    
265
    Parameters:
266
    - geometry (dict): GeoJSON-like geometry
267
    
268
    Returns:
269
    tuple: (left, bottom, right, top) bounds
270
    """
271
```
272

273
Usage examples:
274

275
```python
276
from rasterio.features import sieve, bounds
277

278
# Remove small polygons from classification
279
with rasterio.open('classification.tif') as dataset:
280
    classified = dataset.read(1)
281
    
282
    # Remove polygons smaller than 100 pixels
283
    cleaned = sieve(
284
        classified, 
285
        size=100,
286
        connectivity=8  # 8-connected neighborhoods
287
    )
288
    
289
    # Save cleaned classification
290
    profile = dataset.profile.copy()
291
    with rasterio.open('cleaned_classification.tif', 'w', **profile) as dst:
292
        dst.write(cleaned, 1)
293

294
# Calculate geometry bounds
295
polygon = {
296
    'type': 'Polygon', 
297
    'coordinates': [[
298
        [-120.5, 35.2], [-119.8, 35.2],
299
        [-119.8, 35.9], [-120.5, 35.9], 
300
        [-120.5, 35.2]
301
    ]]
302
}
303

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

308
### Feature Processing Enumerations
309

310
Enumerations for controlling feature processing algorithms:
311

312
```python { .api }
313
class MergeAlg(Enum):
314
    """Algorithms for merging rasterized values."""
315
    
316
    replace = 'replace'  # Later values overwrite earlier (default)
317
    add = 'add'          # Sum overlapping values
318
```
319

320
Usage in rasterization:
321

322
```python
323
from rasterio.enums import MergeAlg
324

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

328
# Sum weights in overlapping areas
329
density_raster = rasterize(
330
    shapes_with_weights,
331
    out_shape=(height, width),
332
    transform=transform,
333
    merge_alg=MergeAlg.add,  # Sum overlapping values
334
    dtype='float32'
335
)
336

337
# Last feature wins in overlapping areas  
338
priority_raster = rasterize(
339
    shapes_with_weights,
340
    out_shape=(height, width), 
341
    transform=transform,
342
    merge_alg=MergeAlg.replace,  # Default behavior
343
    dtype='uint16'
344
)
345
```
346

347
### Integration with Vector Libraries
348

349
Feature operations work seamlessly with popular Python geospatial libraries:
350

351
```python
352
# Integration with GeoPandas
353
import geopandas as gpd
354

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

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

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

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

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

373
# Convert to GeoJSON for rasterio
374
polygon_geojson = polygon.__geo_interface__
375
point_geojson = point.__geo_interface__
376

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