Jak poligonizować raster do kształtnych wielokątów

14

Szukam rozwiązania python typu open source do konwersji rastra na wielokąt (bez ArcPy).

Znałem funkcję GDAL do konwersji rastra na wielokąt, a oto instrukcja: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#polygonize-a-raster-band

Niemniej jednak oczekuję, że dane wyjściowe mogą być kształtnymi wielokątami lub dowolnym obiektem, tymczasowo w pamięci, a nie zapisanym jako plik. Czy jest jakiś pakiet lub kod do obsługi tego problemu?

Jeśli raster został przetworzony w tablicy numpy, podejście jest wymienione poniżej.

Vicky Liau
źródło

Odpowiedzi:

21

Użyj rasterio Seana Gilliesa. Można go łatwo łączyć z Fioną (odczytywać i zapisywać pliki kształtów) i kształtować tego samego autora.

W skrypcie rasterio_polygonize.py początek jest

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

Wynikiem jest generator funkcji 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}}

Że możesz przekształcić się w kształtne geometrie

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

Twórz geopandasową ramkę danych i włącz łatwe w obsłudze funkcje łączenia przestrzennego, kreślenia, zapisywania jako geojson, plik kształtu ESRI itp.

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)
gen
źródło
jeśli raster został przetworzony jako tablica numpy, czy istnieje sposób na przekonwertowanie tablicy numpy na wielokąty? Dzięki!
Vicky Liau,
Teoretycznie tak
gen
1
Zmienna maski i parametr w twoim przykładzie wydają się niepotrzebne. if value > src.nodataPoleciłbym jednak dodać do listy zrozumienie, aby skorzystać z wartości źródłowej źródła i odrzucić wszelkie kształty, które jej odpowiadają. Nie jestem jednak pewien, co by się stało, gdyby nie istniały wartości nodata. : o)
bugmenot123
3
w międzyczasie zmienili rasterio.drivers na rasterio.Env i src.affine na src.transform
Leo
4

Oto moja implementacja.

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

Sposób zainstalowania rasterio to „conda install -c https://conda.anaconda.org/ioos rasterio”, jeśli występuje problem z instalacją.

Vicky Liau
źródło
Wynik rasterio jest bezpośrednio tablicą liczb liczbowych, więc nie potrzebujeszmyarray=srcband.ReadAsArray() #or any array
gen
@gene Poprawiłem notatkę. Ta linia (myarray = srcband.ReadAsArray ()) używa gdal do importowania obrazu.
Vicky Liau,
zaimportuj obraz jako tablicę numpy, a rasterio zaimportuj bezpośrednio obraz jako tablicę numpy
gen
To działało dla mnie, chociaż musiałem indeksować vec, ponieważ był zwracany jako krotka. shape (vec [0])
user2723146