Menentukan apakah shapefile dan raster tumpang tindih dalam Python menggunakan OGR / GDAL? [Tutup]

9

Saya sedang membangun skrip dengan python menggunakan OGR / GDAL.

Saya memiliki satu set shapefile dan satu set file raster GeoTiff.

Saya ingin skrip abaikan shapefile saya jika tidak bersinggungan dengan area raster.

Shapefile bukan persegi panjang, jadi saya tidak bisa hanya membandingkan nilai xmin / xmax, ymin / ymax yang dikembalikan oleh layer.GetExtent (). Saya membutuhkan poligon aktual yang merepresentasikan bentuk keseluruhannya, dan kemudian beberapa cara untuk menentukan apakah poligon itu bersilangan dengan raster square.

Saya berpikir saya entah bagaimana bisa menggabungkan semua poligon di shapefile menjadi satu fitur, dan kemudian membaca geometri pada fitur itu, dan kemudian membandingkan informasi itu dengan tingkat raster. Namun, saya tidak yakin secara khusus bagaimana melakukan ini.

  1. Bagaimana cara mengekstrak informasi poligon perbatasan dari shapefile?
  2. Bagaimana cara menentukan apakah poligon itu memotong area persegi yang diberikan?
python raster shapefile gdal ogr JFerg
sumber
Saya tidak terbiasa dengan osgeo, tetapi setara arcpy akan (mungkin) melibatkan: membaca tingkat raster, membuat poligon yang mencakup tingkat dalam memori, siklus melalui shapefile, memotong masing-masing ke tingkat persegi panjang, menguji apakah ada hasil.
Floem

Jawaban:

17

Script berikut menentukan kotak pembatas dari raster dan membuat berdasarkan kotak pembatas sebuah 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)

Selanjutnya, geometri poligon vektor ditentukan. Ini menjawab pertanyaan pertama Anda.

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

Terakhir, geometri vektor dan raster diuji untuk persimpangan (pengembalian Trueatau False). Ini menjawab pertanyaan kedua Anda. 

print rasterGeometry.Intersect(vectorGeometry)
ustroetz
sumber
2
Terima kasih, ini persis apa yang saya cari. Jawaban ini jelas menunjukkan cara membuat, mengekstrak, dan menjalankan fungsi antara objek geometri, yang persis apa yang saya cari.
JFerg
Solusi ini menguji apakah FID = 0 poligon bersilangan dengan raster. Bagaimana Anda mendapatkan geometri total shapefile ketika tidak ada poligon yang mewakili ini?
JFerg
1
EDIT: Peningkatan waktu komputasi tidak penting, jadi saya memeriksa apakah setiap poligon di shapefile berpotongan sekarang.
JFerg
4

Saya menemukan solusi @ustroetz sangat membantu tetapi perlu diperbaiki di dua tempat. Pertama, pixelHeight = transform [5] sudah bernilai negatif, jadi persamaannya seharusnya

yBottom = yTop+rows*pixelHeight

Kedua, urutan titik-titik dalam cincin harus berlawanan arah jarum jam. Saya mengalami masalah dengan itu. Pesanan yang benar adalah:

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
sumber