So polygonisieren Sie Raster in formschöne Polygone

13

Ich suche eine Open-Source-Python-Lösung zum Konvertieren von Raster in Polygon (kein ArcPy).

Ich kannte die GDAL-Funktion zum Konvertieren von Raster in Polygon, und hier ist das Handbuch: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#polygonize-a-raster-band

Trotzdem erwarte ich, dass es sich bei der Ausgabe um formschöne Polygone oder ein beliebiges Objekt handeln kann, das sich vorübergehend im Speicher befindet und nicht als Datei gespeichert wird. Gibt es ein Paket oder einen Code, um dieses Problem zu beheben?

Wenn das Raster in einem Numpy-Array verarbeitet wurde, ist der Ansatz unten aufgeführt.

Vicky Liau
quelle

Antworten:

18

Verwenden Sie das Raster von Sean Gillies. Es kann leicht mit Fiona (Lesen und Schreiben von Shapefiles) und Shapely desselben Autors kombiniert werden .

Im Skript rasterio_polygonize.py ist der Anfang

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

Das Ergebnis ist ein Generator für GeoJSON-Funktionen

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

Dass Sie sich in formschöne Geometrien verwandeln können

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

Erstellen Sie Geopandas Dataframe und aktivieren Sie benutzerfreundliche Funktionen zum räumlichen Verbinden, Plotten, Speichern als Geojson, ESRI-Shapefile usw.

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)
Gen
quelle
Gibt es eine Möglichkeit, ein Numpy-Array als Polygone zu konvertieren, wenn das Raster als Numpy-Array verarbeitet wurde? Vielen Dank!
Vicky Liau
Theoretisch ja-
Gen
1
Die Maskenvariable und der Parameter in Ihrem Beispiel scheinen nicht erforderlich zu sein. Ich würde jedoch empfehlen, if value > src.nodatadas Listenverständnis zu erweitern, um den Nodata-Wert der Quelle zu nutzen und alle Formen zu verwerfen, die diesem Wert entsprechen. Ich bin mir nicht sicher, was passieren würde, wenn es keine nodata vaue gibt. : o)
bugmenot123
3
In der Zwischenzeit haben sie rasterio.drivers zu rasterio.Env und src.affine zu src.transform geändert
Leo
3

Hier ist meine Implementierung.

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

Die Installation von rasterio erfolgt über 'conda install -c https://conda.anaconda.org/ioos rasterio', wenn ein Installationsproblem vorliegt.

Vicky Liau
quelle
Das Ergebnis von rasterio ist direkt ein Numpy-Array, das Sie also nicht brauchenmyarray=srcband.ReadAsArray() #or any array
gene
@gene Ich habe den Hinweis überarbeitet. Diese Zeile (myarray = srcband.ReadAsArray ()) verwendet gdal, um das Bild zu importieren.
Vicky Liau
Importieren Sie das Bild als Numpy-Array und rasterio importieren Sie das Bild direkt als Numpy-Array
Gen
Dies funktionierte für mich, obwohl ich vec indizieren musste, weil es als Tupel zurückgegeben wurde. Form (vec [0])
Benutzer2723146