1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
|
#!/usr/bin/env python
#
# sieve: demonstrate sieving and polygonizing of raster features.
import subprocess
import numpy as np
import rasterio
from rasterio.features import sieve, shapes
# Register GDAL and OGR drivers.
with rasterio.Env():
# Read a raster to be sieved.
with rasterio.open('tests/data/shade.tif') as src:
shade = src.read(1)
# Print the number of shapes in the source raster.
print("Slope shapes: %d" % len(list(shapes(shade))))
# Sieve out features 13 pixels or smaller.
sieved = sieve(shade, 13, out=np.zeros(src.shape, src.dtypes[0]))
# Print the number of shapes in the sieved raster.
print("Sieved (13) shapes: %d" % len(list(shapes(sieved))))
# Write out the sieved raster.
kwargs = src.meta
kwargs['transform'] = rasterio.transform.guard_transform(kwargs['transform'])
with rasterio.open('example-sieved.tif', 'w', **kwargs) as dst:
dst.write(sieved, indexes=1)
# Dump out gdalinfo's report card and open (or "eog") the TIFF.
print(subprocess.check_output(
['gdalinfo', '-stats', 'example-sieved.tif']))
subprocess.call(['open', 'example-sieved.tif'])
|
