¿Cómo convertir un raster en polígonos de Shapely?

Question

Estoy buscando una solución de Python de código abierto para convertir raster a polígono (sin ArcPy).

Conocía la función de GDAL para convertir raster a polígono, y aquí está el manual: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#polygonize-a-raster-band

Sin embargo, espero que la salida pueda ser polígonos Shapely u cualquier objeto, temporalmente en la memoria, no guardado como un archivo. ¿Existe algún paquete o código para manejar este problema?

Si el raster ha sido procesado en un array de numpy, el enfoque está listado abajo.

Answer 1

Utiliza rasterio de Sean Gillies. Puede combinarse fácilmente con Fiona (lectura y escritura de shapefiles) y con shapely del mismo autor.

En el script rasterio_polygonize.py el comienzo es

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('a_raster') as src:
        image = src.read(1) # primera banda
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))

El resultado es un generador de entidades GeoJSON

 geoms = list(results)
 # primera entidad
 print geoms[0]
 {'geometry': {'type': 'Polygon', 'coordinates': [[(202086.577, 90534.3504440678), (202086.577, 90498.96207), (202121.96537406777, 90498.96207), (202121.96537406777, 90534.3504440678), (202086.577, 90534.3504440678)]]}, 'properties': {'raster_val': 170.52000427246094}}

Que puedes transformar en geometrías shapely

from shapely.geometry import shape
print shape(geoms[0]['geometry'])
POLYGON ((202086.577 90534.35044406779, 202086.577 90498.96206999999, 202121.9653740678 90498.96206999999, 202121.9653740678 90534.35044406779, 202086.577 90534.35044406779))

Crea un Dataframe geopandas y habilita funcionalidades fáciles de usar como unión espacial, gráficos, guardar como geojson, ESRI shapefile, etc.

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

Answer 2

Aquí está mi implementación.

from osgeo import ogr, gdal, osr
from osgeo.gdalnumeric import *  
from osgeo.gdalconst import * 
import fiona
from shapely.geometry import shape
import rasterio.features

#segimg=glob.glob('Poly.tif')[0]
#src_ds = gdal.Open(segimg, GA_ReadOnly )
#srcband=src_ds.GetRasterBand(1)
#myarray=srcband.ReadAsArray() 
#estas líneas utilizan gdal para importar una imagen. 'myarray' puede ser cualquier matriz numpy

mypoly=[]
for vec in rasterio.features.shapes(myarray):
    mypoly.append(shape(vec))

La forma de instalar rasterio es mediante 'conda install -c https://conda.anaconda.org/ioos rasterio', si hay algún problema de instalación.

Answer 3

# PYTHON 3

import rasterio.features
from shapely.geometry import shape

# Mask es una matriz numpy de máscara binaria cargada en la forma que sea necesaria
maskShape = rasterio.features.shapes(mask.astype('uint8'))
mypoly=[]
for vec in maskShape:
    mypoly.append(vec[0])
print(mypoly)

Esto es lo que funcionó para mí.

Answer 4

En caso de píxeles con alta variabilidad, convierta los arrays al tipo entero. Ayuda a ahorrar tiempo de cálculo.

with rio.Env():
with rio.open(input_file_path) as depth_image:
    depth_array=depth_image.read(1).astype(int)