Cómo poligonalizar ráster a polígonos bien formados

13

Estoy buscando una solución de python de código abierto para convertir el ráster a polígono (sin ArcPy).

Conocía la función GDAL para convertir el ráster 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 bien formados o cualquier objeto, temporalmente en la memoria, no guardado como un archivo. ¿Hay algún paquete o código para manejar este problema?

Si el ráster se ha procesado en una matriz numpy, el enfoque se enumera a continuación.

python-2.7 shapely polygonize Vicky Liau
fuente

Respuestas:

20

Usa rasterio de Sean Gillies. Se puede combinar fácilmente con Fiona (leer y escribir archivos de forma) y bien formado del mismo autor.

En el script rasterio_polygonize.py el comienzo es

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

El resultado es un generador de características GeoJSON.

 geoms = list(results)
 # first feature
 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 bien formadas

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))

Cree geopandas Dataframe y permita funcionalidades fáciles de usar de unión espacial, trazado, guardar como geojson, archivo de forma ESRI, etc.

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)
gene
fuente
Si el ráster se ha procesado como una matriz numpy, ¿hay alguna forma de convertir una matriz numpy como polígonos? ¡Gracias!
Vicky Liau
En teoría, sí
gen
1
La variable de máscara y el parámetro en su ejemplo parecen innecesarios. Sin embargo, recomendaría agregar if value > src.nodataa la lista la comprensión para hacer uso del valor de nodatos de la fuente y descartar cualquier forma que le corresponda. Sin embargo, no estoy seguro de lo que sucedería si no hay un valor de nodata. : o)
bugmenot123
3
mientras tanto cambiaron rasterio.drivers a rasterio.Env y src.affine a src.transform
Leo
3

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() 
#these lines use gdal to import an image. 'myarray' can be any numpy array

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

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

Vicky Liau
fuente
El resultado de rasterio es directamente una matriz numpy, por lo tanto no es necesariomyarray=srcband.ReadAsArray() #or any array
gen
@gene revisé la nota. Esta línea (myarray = srcband.ReadAsArray ()) usa gdal para importar la imagen.
Vicky Liau
importe la imagen como una matriz numpy y rasterio importe directamente la imagen como una matriz numpy
gen
Esto funcionó para mí, aunque tuve que indexar vec porque estaba siendo devuelto como una tupla. shape (vec [0])
usuario2723146