Déterminer si le fichier de formes et le raster se chevauchent en Python en utilisant OGR / GDAL? [fermé]

9

Je construis un script en python en utilisant OGR / GDAL.

J'ai un ensemble de fichiers de formes et un ensemble de fichiers raster GeoTiff.

J'aimerais que mon script ignore les fichiers de formes s'ils ne se croisent pas avec la zone raster.

Le fichier de formes n'est pas un rectangle, donc je ne peux pas simplement comparer les valeurs xmin / xmax, ymin / ymax retournées par layer.GetExtent (). J'ai besoin du polygone réel représentant sa forme globale, puis d'un moyen de déterminer si ce polygone croise le carré raster.

Je pensais que je pouvais en quelque sorte fusionner tous les polygones du fichier de formes en une seule entité, puis lire la géométrie de cette entité, puis comparer ces informations à l'étendue raster. Cependant, je ne sais pas exactement comment exécuter cela.

  1. Comment extraire des informations de polygone de bordure à partir d'un fichier de formes?
  2. Comment déterminer si ce polygone coupe une zone carrée donnée?
python raster shapefile gdal ogr JFerg
la source
Je ne connais pas osgeo, mais l'équivalent arcpy impliquerait (pourrait) impliquer: lire l'étendue raster, créer une étendue de polygone couvrant en mémoire, parcourir les fichiers de formes, découper chacun pour étendre le rectangle, tester si quelque chose en résulte.
phloème

Réponses:

17

Le script suivant détermine le cadre de sélection d'un raster et crée en fonction du cadre de sélection une géométrie. 

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)

Ensuite, la géométrie du polygone vectoriel est déterminée. Cela répond à votre première question.

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

Enfin, la géométrie du vecteur et du raster est testée pour l'intersection (renvoie Trueou False). Cela répond à votre deuxième question. 

print rasterGeometry.Intersect(vectorGeometry)
ustroetz
la source
2
Merci, c'était exactement ce que je cherchais. Cette réponse montre clairement comment créer, extraire et exécuter des fonctions entre des objets géométriques, ce qui est exactement ce que je cherchais.
JFerg
Cette solution teste si le polygone FID = 0 recoupe le raster. Comment obtenir la géométrie du total du fichier de formes quand aucun polygone ne représente cela?
JFerg
1
EDIT: L'augmentation du temps de calcul est sans conséquence, donc je vérifie si chaque polygone du fichier de formes se coupe maintenant.
JFerg
4

Je trouve la solution @ustroetz très utile mais elle devait être corrigée à deux endroits. Premièrement, pixelHeight = transform [5] est déjà une valeur négative, donc l'équation doit être

yBottom = yTop+rows*pixelHeight

Deuxièmement, l'ordre des points dans l'anneau doit être dans le sens inverse des aiguilles d'une montre. J'avais des problèmes avec ça. L'ordre correct est:

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
la source