Cortar un polígono en cuadrados más pequeños en Python

Question

Para analizar una imagen GeoTif nos gustaría cortar una gran característica poligonal de un shp en trozos más pequeños (esto se debe a que por alguna razón nuestro análisis no se realiza bien con un polígono de gran extensión).

Lo que queremos obtener exactamente es que para un polígono dado P , ancho_de_ventana y altura_de_la_ventana queremos obtener una lista de características de polígonos donde cada característica es una intersección de un P con un conjunto de ventanas disjuntas y adyacentes que tiene forma de rectángulo con ancho_de_ventana y altura_de_la_ventana .

Esta imagen muestra lo que queremos decir:

Para este polígono nos gustaría obtener una lista de 12 nuevas características del polígono donde el primero es un triángulo, el segundo - un pentágono, el tercero - un cuadrado y así sucesivamente.

Lo importante es que queremos tener esto como Python función. Así que nos gustaría tener la siguiente función:

cut_polygon_into_windows(p, window_height, window_width)

para ser añadido a nuestra API.

Answer 1

Mi solución parte de la creación previa de una cuadrícula (como capa de memoria) que cubre la extensión de su shapefile de entrada: para ello, he adaptado un código tomado de este blog .

Una vez creada la cuadrícula, intersecciono cada característica de la misma con la capa de entrada.

Partiendo de este shapefile de muestra:

y usando este código (lo formateé como una función, aunque no es muy eficiente agrupar todas estas líneas):

from qgis.core import *
from qgis.PyQt.QtCore import QVariant

def cut_polygon_into_windows(p, window_height, window_width):

    crs = p.crs().toWkt()
    extent = p.extent()
    (xmin, xmax, ymin, ymax) = (extent.xMinimum(), extent.xMaximum(), extent.yMinimum(), extent.yMaximum())

    # Create the grid layer
    vector_grid = QgsVectorLayer('Polygon?crs='+ crs, 'vector_grid' , 'memory')
    prov = vector_grid.dataProvider()

    # Create the grid layer
    output = QgsVectorLayer('Polygon?crs='+ crs, 'output' , 'memory')
    outprov = output.dataProvider()

    # Add ids and coordinates fields
    fields = QgsFields()
    fields.append(QgsField('ID', QVariant.Int, '', 10, 0))
    outprov.addAttributes(fields)

    # Generate the features for the vector grid
    id = 0
    y = ymax
    while y >= ymin:
        x = xmin
        while x <= xmax:
            point1 = QgsPoint(x, y)
            point2 = QgsPoint(x + window_width, y)
            point3 = QgsPoint(x + window_width, y - window_height)
            point4 = QgsPoint(x, y - window_height)
            vertices = [point1, point2, point3, point4] # Vertices of the polygon for the current id
            inAttr = [id]
            feat = QgsFeature()
            feat.setGeometry(QgsGeometry().fromPolygon([vertices])) # Set geometry for the current id
            feat.setAttributes(inAttr) # Set attributes for the current id
            prov.addFeatures([feat])
            x = x + window_width
            id += 1
        y = y - window_height

    index = QgsSpatialIndex() # Spatial index
    for ft in vector_grid.getFeatures():
        index.insertFeature(ft)

    for feat in p.getFeatures():
        geom = feat.geometry()
        idsList = index.intersects(geom.boundingBox())
        for gridfeat in vector_grid.getFeatures(QgsFeatureRequest().setFilterFids(idsList)):
            tmp_geom = QgsGeometry(gridfeat.geometry())
            tmp_attrs = gridfeat.attributes()
            if geom.intersects(tmp_geom):
                int = QgsGeometry(geom.intersection(tmp_geom))
                outfeat = QgsFeature()
                outfeat.setGeometry(int)
                outfeat.setAttributes(tmp_attrs)
                outprov.addFeatures([outfeat])

    output.updateFields()

    return output

# Load the layer
p = QgsVectorLayer('C:/path_to_your_file/input.shp', 'display name', 'ogr')

# Set width and height as you want
window_width = 300
window_height = 300

# Run the function
output = cut_polygon_into_windows(p, window_height, window_width)

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayers([output])

Obtengo el resultado deseado: