import fiona
shp = fiona.open('zonas_farrapos.shp')
res = 72
width = (shp.bounds[2] - shp.bounds[0]) // res
height = (shp.bounds[3] - shp.bounds[1]) // res
out_shape = int(height), int(width)
out_shape
(714, 429)
from rasterio.transform import from_origin
transform = from_origin(shp.bounds[0] - res / 2,
shp.bounds[3] + res / 2, res, res)
transform
Affine(72.0, 0.0, 391451.64563056955, 0.0, -72.0, 6391938.555971239)
shapes = [(geometry['geometry'], k) for k, geometry in shp.items()]
#colors = (255, 0, 50, 0)
#shapes = [(geometry['geometry'], color) for geometry, color in zip(list(shp.values()), colors)]
import numpy as np
import rasterio
from rasterio.features import rasterize
dtype = rasterio.float64
fill = np.nan
raster = rasterize(
shapes,
out_shape=out_shape,
fill=fill,
transform=transform,
dtype=dtype
)
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(9, 9))
ax.imshow(raster)
ax.axis('off');
Salvar o geotiff.
with rasterio.open(
'test.tif', 'w',
crs=shp.crs,
driver='GTiff',
dtype=dtype,
count=3,
width=width,
height=height,
nodata=fill,
transform=transform
) as dst:
dst.write(raster, indexes=1)
Round-trip
with rasterio.open('test.tif') as src:
r, g, b = src.read()
print(src.width, src.height)
print(src.crs)
429 714 CRS({'init': 'epsg:32721'})
total = np.zeros(r.shape)
for band in r, g, b:
total += band
total /= 3
Aparentemente a figura "preta" é pq não tenho um visualizador corrento.
A tif
parece OK mesmo salvando em 3 bandas.
plt.imshow(total)
<matplotlib.image.AxesImage at 0x7f429fc60390>