Créer un fichier de formes de polygone de grille carrée avec Python?

9

J'ai les coordonnées suivantes

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

Je souhaite créer une grille carrée de taille 1 m en 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
la source
Est-ce attaché à une plate-forme SIG spécifique ou est-il obligatoire de le faire en python pur sans aucun format de sortie spécifié (par exemple, shapefile, fichier texte, etc.)
Merci @Dan, je souhaite jouer en pur python et la sortie sera au format shapefile
Gianni
Le niveau de licence ArcInfo d'ArcMap comprend l'outil Fishnet, mais vous n'avez pas indiqué comment vous envisagez de créer le fichier de formes.
Désolé, je n'utilise pas de logiciel commercial. Je préfère le programme en langage pur Java, Python, C ++.
Gianni
1
Mais cela ne vous dérange pas d'utiliser une bibliothèque telle que GDAL / OGR ( pypi.python.org/pypi/GDAL ) ou pyshp ( pypi.python.org/pypi/pyshp )?
Snorfalorpagus

Réponses:

11

Le script suivant fera le travail avec GDAL et 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
la source
6

Ce script Python utilise la bibliothèque pyshp , comme suggéré par user16044:

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')

Remarque: une grille carrée de taille 1 m avec une telle étendue équivaut à une couche contenant environ 1 million de polygones et donc les performances du script diminuent sensiblement.

Antonio Falciano
la source
1

Il y a une réponse à cette question depuis longtemps, mais j'ajoute une autre solution en utilisant les bibliothèques galbées et 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
la source
-1

La réponse à la création de Shapefile de grille résille dans QGIS? affiche une option de création de grille dans la boîte à outils de traitement QGIS.

user965586
la source
OP a déclaré qu'il préfère un exemple en Python pur plutôt qu'avec un logiciel
LaughU
Compte tenu des autres réponses, importez des bibliothèques, l'importation de modules QGIS est un moyen valable pour éviter une interface graphique. Comme dans cette réponse gis.stackexchange.com/questions/79916/…
user965586
Les autres réponses fournissent du code, donc si le vôtre le faisait aussi, je pense qu'il serait mieux reçu.
PolyGeo