Определение перекрытия шейп-файла и растра в Python с использованием OGR / GDAL? [закрыто]

9

Я строю сценарий на Python с использованием OGR / GDAL.

У меня есть набор шейп-файлов и набор растровых файлов GeoTiff.

Я хотел бы, чтобы мой скрипт игнорировал шейп-файлы, если они не пересекаются с растровой областью.

Шейп-файл не является прямоугольником, поэтому я не могу просто сравнить значения xmin / xmax, ymin / ymax, возвращаемые layer.GetExtent (). Мне нужен фактический многоугольник, представляющий его общую форму, а затем какой-то способ определить, пересекается ли этот многоугольник с растровым квадратом.

Я думал, что смогу каким-то образом объединить все полигоны в шейп-файле в один объект, а затем прочитать геометрию этого объекта и затем сравнить эту информацию с растровым экстентом. Тем не менее, я не уверен, как именно это выполнить.

  1. Как извлечь информацию о полигоне границы из шейп-файла?
  2. Как определить, пересекает ли этот многоугольник данный квадрат?
python raster shapefile gdal ogr JFerg
источник
Я не знаком с osgeo, но дугообразный эквивалент будет (может) включать: считывание растрового экстента, создание полигона, охватывающего экстент в памяти, циклический переход по шейп-файлам, обрезание каждого до экстента прямоугольника, тестирование, если что-нибудь получится.
флоэма

Ответы:

17

Следующий скрипт определяет ограничивающую рамку растра и создает на основе ограничивающей рамки геометрию. 

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)

Далее определяется геометрия векторного многоугольника. Это отвечает на ваш первый вопрос.

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

Наконец, геометрия вектора и растра проверяются на пересечение (возврат Trueили False). Это отвечает на ваш второй вопрос. 

print rasterGeometry.Intersect(vectorGeometry)
ustroetz
источник
2
Спасибо, это было именно то, что я искал. Этот ответ ясно показывает, как создавать, извлекать и запускать функции между геометрическими объектами, и это именно то, что я искал.
JFerg
Это решение проверяет, пересекается ли полигон FID = 0 с растром. Как получить геометрию итогового шейп-файла, когда ни один полигон не представляет это?
JFerg
1
РЕДАКТИРОВАТЬ: Увеличение времени вычислений несущественно, поэтому я проверяю, пересекается ли каждый многоугольник в шейп-файле.
JFerg
4

Я нахожу решение @ustroetz очень полезным, но его нужно исправить в двух местах. Во-первых, pixelHeight = transform [5] уже является отрицательным значением, поэтому уравнение должно быть

yBottom = yTop+rows*pixelHeight

Во-вторых, порядок точек в кольце должен быть против часовой стрелки. У меня были проблемы с этим. Правильный порядок:

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)
Пандза
источник