¿Crear un polígono de forma cuadriculada con Python?

9

Tengo las siguientes coordenadas

minx, maxx, miny ,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073

Deseo crear una cuadrícula cuadrada de tamaño 1 m usando python.

import math


minx,maxx,miny,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073
size = 1

def set_bbox(minx, maxx, miny, maxy, distx, disty):
    nx = int(math.ceil(abs(maxx - minx)/distx))
    ny = int(math.ceil(abs(maxy - miny)/disty))
    new_maxx = minx + (nx*distx)
    new_miny = maxy - (ny*disty)
    return ((minx, new_maxx, new_miny, maxy),ny,nx)

# shift the bottom (right - down)
coord, ny, nx = set_bbox(minx,maxx,miny,maxy,size,size)
# left-up origin
origin = coord[0],coord[3]
# number of tiles
ncell = ny*nx
Gianni
fuente
¿Está conectado a una plataforma SIG específica o es el requisito para hacerlo en Python puro sin ningún formato de salida especificado (por ejemplo, shapefile, textfile, etc.)
Gracias @Dan, deseo actuar en Python puro y la salida estará en formato shapefile
Gianni
El nivel de licencia ArcInfo de ArcMap tiene la herramienta Fishnet, pero no ha indicado cómo piensa crear el archivo de forma.
Lo siento, no uso software comercial. Prefiero el programa en lenguaje puro Java, Python, C ++.
Gianni
1
¿Pero no le importa usar una biblioteca como GDAL / OGR ( pypi.python.org/pypi/GDAL ) o pyshp ( pypi.python.org/pypi/pyshp )?
Snorfalorpagus

Respuestas:

11

El siguiente script hará el trabajo con GDAL y Python:

import os, sys
import ogr
from math import ceil

def main(outputGridfn,xmin,xmax,ymin,ymax,gridHeight,gridWidth):

    # convert sys.argv to float
    xmin = float(xmin)
    xmax = float(xmax)
    ymin = float(ymin)
    ymax = float(ymax)
    gridWidth = float(gridWidth)
    gridHeight = float(gridHeight)

    # get rows
    rows = ceil((ymax-ymin)/gridHeight)
    # get columns
    cols = ceil((xmax-xmin)/gridWidth)

    # start grid cell envelope
    ringXleftOrigin = xmin
    ringXrightOrigin = xmin + gridWidth
    ringYtopOrigin = ymax
    ringYbottomOrigin = ymax-gridHeight

    # create output file
    outDriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputGridfn):
        os.remove(outputGridfn)
    outDataSource = outDriver.CreateDataSource(outputGridfn)
    outLayer = outDataSource.CreateLayer(outputGridfn,geom_type=ogr.wkbPolygon )
    featureDefn = outLayer.GetLayerDefn()

    # create grid cells
    countcols = 0
    while countcols < cols:
        countcols += 1

        # reset envelope for rows
        ringYtop = ringYtopOrigin
        ringYbottom =ringYbottomOrigin
        countrows = 0

        while countrows < rows:
            countrows += 1
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            # add new geom to layer
            outFeature = ogr.Feature(featureDefn)
            outFeature.SetGeometry(poly)
            outLayer.CreateFeature(outFeature)
            outFeature.Destroy

            # new envelope for next poly
            ringYtop = ringYtop - gridHeight
            ringYbottom = ringYbottom - gridHeight

        # new envelope for next poly
        ringXleftOrigin = ringXleftOrigin + gridWidth
        ringXrightOrigin = ringXrightOrigin + gridWidth

    # Close DataSources
    outDataSource.Destroy()


if __name__ == "__main__":

    #
    # example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
    #

    if len( sys.argv ) != 8:
        print "[ ERROR ] you must supply seven arguments: output-shapefile-name.shp xmin xmax ymin ymax gridHeight gridWidth"
        sys.exit( 1 )

    main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7] )
ustroetz
fuente
6

Este script de Python usa la biblioteca pyshp , como lo sugiere el usuario16044:

import shapefile as shp
import math

minx,maxx,miny,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073
dx = 100
dy = 100

nx = int(math.ceil(abs(maxx - minx)/dx))
ny = int(math.ceil(abs(maxy - miny)/dy))

w = shp.Writer(shp.POLYGON)
w.autoBalance = 1
w.field("ID")
id=0

for i in range(ny):
    for j in range(nx):
        id+=1
        vertices = []
        parts = []
        vertices.append([min(minx+dx*j,maxx),max(maxy-dy*i,miny)])
        vertices.append([min(minx+dx*(j+1),maxx),max(maxy-dy*i,miny)])
        vertices.append([min(minx+dx*(j+1),maxx),max(maxy-dy*(i+1),miny)])
        vertices.append([min(minx+dx*j,maxx),max(maxy-dy*(i+1),miny)])
        parts.append(vertices)
        w.poly(parts)
        w.record(id)

w.save('polygon_grid')

Nota: una cuadrícula cuadrada de tamaño 1 m con tal extensión equivale a una capa que contiene aproximadamente 1 millón de polígonos y, por lo tanto, el rendimiento del guión disminuye sensiblemente.

Antonio Falciano
fuente
1

Esta pregunta se respondió desde hace tiempo, pero agrego otra solución usando bibliotecas bien formadas y fiona:

import fiona
from shapely.geometry import mapping, LineString, MultiLineString

file = 'input.shp'
with fiona.open(file, 'r') as ds_in:
    num_tiles = 5
    schema = {
    "geometry": "MultiLineString",
    "properties": {"id": "int"}
     }
minx, miny, maxx, maxy = ds_in.bounds
dx = (maxx - minx) / num_tiles
dy = (maxy - miny) / num_tiles

lines = []
for x in range(num_tiles + 1):
    lines.append(LineString([(minx + x * dx, miny), (minx + x * dx, maxy)]))
for y in range(num_tiles + 1):
    lines.append(LineString([(minx, miny + y * dy), (maxx, miny + y * dy)]))
grid = MultiLineString(lines)
out = 'gridtest.shp'
with fiona.open(out, 'w', driver=ds_in.driver, schema=schema, crs=ds_in.crs) as ds_dst:
    ds_dst.write({'geometry': mapping(grid), "properties": {"id": 0}})
ImanolUr
fuente
-1

¿La respuesta a Crear un Shapefile de cuadrícula de red en QGIS? muestra una opción de crear cuadrícula en la caja de herramientas de procesamiento de QGIS.

usuario965586
fuente
OP declaró que él / ella prefiere un ejemplo en Python puro en lugar de un software
LaughU
Dadas las otras respuestas de importación de bibliotecas, la importación de módulos QGIS es una forma válida de evitar una GUI. Como en esta respuesta gis.stackexchange.com/questions/79916/…
user965586
Las otras respuestas proporcionan código, por lo que si el suyo también lo hizo, creo que sería mejor recibirlo.
PolyGeo