¿Determinar si el archivo de forma y el ráster se superponen en Python usando OGR / GDAL? [cerrado]

9

Estoy construyendo un script en python usando OGR / GDAL.

Tengo un conjunto de archivos de forma y un conjunto de archivos ráster GeoTiff.

Me gustaría que mi script ignore los archivos de forma si no se cruzan con el área de la trama.

El shapefile no es un rectángulo, por lo que no puedo simplemente comparar los valores xmin / xmax, ymin / ymax devueltos por layer.GetExtent (). Necesito que el polígono real represente su forma general, y luego alguna forma de determinar si ese polígono se cruza con el cuadrado ráster.

Estaba pensando que de alguna manera podría fusionar todos los polígonos en el archivo de forma en una entidad, y luego leer la geometría de esa entidad, y luego comparar esa información con la extensión ráster. Sin embargo, no estoy seguro de cómo ejecutar esto específicamente.

  1. ¿Cómo extraer información de polígono de borde del archivo shape?
  2. ¿Cómo determinar si ese polígono se cruza con un área cuadrada dada?
python raster shapefile gdal ogr JFerg
fuente
No estoy familiarizado con osgeo, pero el equivalente de arcpy implicaría (podría) involucrar: leer la extensión de la trama, crear la extensión de cobertura del polígono en la memoria, recorrer los archivos de forma, recortar cada rectángulo de extensión, probar si algo resulta.
floema

Respuestas:

17

El siguiente script determina el cuadro delimitador de un ráster y crea una geometría basada en el cuadro delimitador. 

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)

A continuación, se determina la geometría del polígono vectorial. Esto responde a tu primera pregunta.

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

Por último, se prueba la intersección de la geometría del vector y el ráster (retornos Trueo False). Esto responde a tu segunda pregunta. 

print rasterGeometry.Intersect(vectorGeometry)
ustroetz
fuente
2
Gracias, esto era exactamente lo que estaba buscando. Esta respuesta muestra claramente cómo crear, extraer y ejecutar funciones entre objetos de geometría, que es exactamente lo que estaba buscando.
JFerg
Esta solución prueba si el polígono FID = 0 se cruza con el ráster. ¿Cómo se obtiene la geometría del total del archivo de forma cuando ningún polígono representa esto?
JFerg
1
EDITAR: El aumento en el tiempo de cálculo es intrascendente, así que verifico si cada polígono en el archivo de forma se cruza ahora.
JFerg
4

Encuentro la solución @ustroetz muy útil, pero necesitaba corregirse en dos lugares. En primer lugar, pixelHeight = transform [5] ya es un valor negativo, por lo que la ecuación debería ser

yBottom = yTop+rows*pixelHeight

En segundo lugar, el orden de los puntos en el anillo debe ser en sentido antihorario. Estaba teniendo problemas con eso. El orden correcto es:

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
fuente