GDAL poligonizuje w pythonie tworząc pusty wielokąt?

12

Mam problem z użyciem funkcji Polygonize w Pythonie. Przykład książki kucharskiej na ten temat można znaleźć tutaj .

Odpowiednia część mojego kodu to:

sourceRaster = gdal.Open('myraster.tif')
band = sourceRaster.GetRasterBand(1)
bandArray = band.ReadAsArray()
outShapefile = "polygonized"
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outShapefile+".shp"):
    driver.DeleteDataSource(outShapefile+".shp")
outDatasource = driver.CreateDataSource(outShapefile+ ".shp")
outLayer = outDatasource.CreateLayer("polygonized", srs=None)
gdal.Polygonize( band, None, outLayer, -1, [], callback=None )
outDatasource.Destroy()
sourceRaster = None

Wiem, że zespół ma odpowiednie informacje, oto fragment bandArray:

array([[ 4.,  4.,  3.,  3.,  3.,  2.,  2.,  2.,  2.,  3.,  3.,  3.,  3.,
         3.,  3.,  3.,  3.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,
         4.,  4.,  4.,  4.],

Kiedy otwieram tabelę atrybutów w QGIS, jest pusta: Screencapture QGIS

Edytować:

Konwersja działa dobrze w QGIS przy użyciu narzędzia Raster -> Konwersja -> Narzędzie poligonizacji

Zrzut ekranu rastra do poligonizacji:

raster do poligonizacji

I zrzut ekranu wynikowej konwersji z narzędzia QGIS:

poligonizowany raster z narzędzia QGIS

Używam dystrybucji Enthought na Windows 7, GDAL wersja 1.10.0-3

Problem polega na tym, że nie mogę poligonizować rastra w pythonie za pomocą GDAL i przykładu książki kucharskiej, mogę poligonizować ten sam raster bez problemu w GUI QGIS

Camdenl
źródło
Jak wygląda twój raster? Czy to naprawdę zawiera wielokąty? Czy to działa, jeśli zamiast tego używasz gdal_polygonize.py?
BradHards,
Edytowano, aby dodać zrzuty ekranu z pracy w QGIS
camdenl
Jaki jest tutaj prawdziwy problem?
Fezter
Dodano konkretny problem
Camdenl
3
Miałem podobny problem (tworzenie pustego pliku kształtu), a utworzenie pola nie pomogło. To, co robiłem źle, to to, że nie zamknąłem pliku kształtu w moim kodzie przed wywołaniem poligonize. Zamykacie to w swoim przykładzie, ja po prostu zamieszczam to dla odniesienia innych.
Stephanie,

Odpowiedzi:

19

Problem polega na tym, że nie tworzyłem pola do przechowywania pasma rastrowego. Po przekopaniu się przez plik gdal_polygonize.py zdałem sobie sprawę, że nie dzieje się to automatycznie podczas wywoływania gdal.Polygonize, który zamiast tego używa funkcji tutaj znalezionej .

Oto dodatkowy krok potrzebny do utworzenia pola i napisania pasma w polu:

newField = ogr.FieldDefn('MYFLD', ogr.OFTInteger)
outLayer.CreateField(newField)

Następnie możemy zapisać pasmo w tym polu o indeksie 0:

gdal.Polygonize(band, None, outLayer, 0, [], callback=None )
Camdenl
źródło
Próbuję również użyć funkcji gdal.Polygonize (), aby uzyskać raster jako wielokąt w pythonie. Ale w linii końcowej pokazuje błąd czasu wykonania !! czemu?
Shiuli Pervin
Działa dobrze z georeferencyjnymi plikami rastrowymi. Rezultatem jest zbyt wiele wielokątów, ale chcę tylko jeden duży wielokąt pokazujący kontur rastra. Czy ktoś ma pojęcie, jak to działa w tym samym czasie?
Shiuli Pervin
Nadal otrzymuję pusty plik shapefile, ale mam wiersze w pliku dbf. Proszę, wyjaśnij mnie!
Satya Chandra
Właśnie miałem ten problem, ale zamiast dodawać fikcyjne pole, możesz wprowadzić indeks -1. Zobacz tutaj, że to pole jest dodawane tylko wtedy, gdy indeks> = 0.
jon_two