Ustalanie, czy plik kształtu i raster pokrywają się w Pythonie przy użyciu OGR / GDAL? [Zamknięte]

9

Buduję skrypt w Pythonie za pomocą OGR / GDAL.

Mam zestaw plików kształtów i zestaw plików rastrowych GeoTiff.

Chciałbym, aby mój skrypt ignorował pliki kształtów, jeśli nie przecinają się one z obszarem rastrowym.

Plik shapefile nie jest prostokątem, więc nie mogę po prostu porównać wartości xmin / xmax, ymin / ymax zwróconych przez layer.GetExtent (). Potrzebuję rzeczywistego wielokąta reprezentującego jego ogólny kształt, a następnie jakiś sposób ustalenia, czy wielokąt przecina się z kwadratem rastrowym.

Myślałem, że mogę jakoś połączyć wszystkie wielokąty w pliku kształtu w jeden element, a następnie odczytać geometrię tego elementu, a następnie porównać te informacje w zakresie rastrowym. Nie jestem jednak pewien, jak to zrobić.

  1. Jak wyodrębnić informacje o wielokącie granicznym z pliku kształtu?
  2. Jak ustalić, czy ten wielokąt przecina dany obszar kwadratowy?
python raster shapefile gdal ogr JFerg
źródło
Nie znam osgeo, ale ekwiwalent arkadowy wymagałby (może): odczytać zasięg rastra, stworzyć zasięg pamięci obejmujący wielobok, cyklicznie zmieniać pliki kształtów, przycinać każdy do prostokąta, sprawdzać, czy coś się pojawi.
łyko

Odpowiedzi:

17

Poniższy skrypt określa obwiednię rastra i tworzy na podstawie obwiedni geometrię. 

import ogr, gdal

raster = gdal.Open('sample.tif')
vector = ogr.Open('sample.shp')

# Get raster geometry
transform = raster.GetGeoTransform()
pixelWidth = transform[1]
pixelHeight = transform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize

xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop-rows*pixelHeight

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)

Następnie określa się geometrię wielokąta wektorowego. To odpowiada na twoje pierwsze pytanie.

# Get vector geometry
layer = vector.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()

Na koniec geometria wektora i rastra jest testowana pod kątem przecięcia (powrotu Truelub False). To odpowiada na twoje drugie pytanie. 

print rasterGeometry.Intersect(vectorGeometry)
ustroetz
źródło
2
Dzięki, właśnie tego szukałem. Ta odpowiedź wyraźnie pokazuje, jak tworzyć, wyodrębniać i uruchamiać funkcje między obiektami geometrycznymi, czego dokładnie szukałem.
JFerg,
To rozwiązanie sprawdza, czy wielokąt FID = 0 przecina się z rastrem. Jak uzyskać geometrię sumy pliku kształtu, jeśli nie reprezentuje tego wielokąt?
JFerg
1
EDYCJA: Wydłużenie czasu obliczeń jest nieistotne, więc sprawdzam, czy każdy wielokąt w pliku kształtu przecina się teraz.
JFerg
4

Uważam, że rozwiązanie @ustroetz jest bardzo pomocne, ale musiało zostać poprawione w dwóch miejscach. Po pierwsze, pixelHeight = transform [5] jest już wartością ujemną, więc równanie powinno być

yBottom = yTop+rows*pixelHeight

Po drugie, kolejność punktów na ringu musi być przeciwna do ruchu wskazówek zegara. Miałem z tym problemy. Prawidłowa kolejność to:

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xLeft, yTop)
Pandza
źródło