Membuat polifon persegi bentuk shapefile dengan Python?

9

Saya memiliki koordinat berikut

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

Saya ingin membuat kotak persegi ukuran 1 m menggunakan 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
sumber
Apakah ini terlampir pada platform GIS tertentu atau merupakan persyaratan untuk melakukan ini dalam python murni tanpa format output yang ditentukan (mis. Shapefile, textfile, dll.)
Terima kasih @Dan, saya ingin tampil dalam python murni dan hasilnya akan dalam format shapefile
Gianni
Tingkat lisensi ArcInfo dari ArcMap memiliki alat Fishnet tetapi Anda belum menunjukkan bagaimana Anda bermaksud membuat shapefile.
Maaf saya tidak menggunakan Perangkat Lunak komersial. Saya lebih suka program dalam bahasa Jawa murni, Python, C ++.
Gianni
1
Tetapi Anda tidak keberatan menggunakan perpustakaan seperti GDAL / OGR ( pypi.python.org/pypi/GDAL ) atau pyshp ( pypi.python.org/pypi/pyshp )?
Snorfalorpagus

Jawaban:

11

Script berikut akan melakukan pekerjaan dengan GDAL dan 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
sumber
6

Skrip Python ini menggunakan pyshp library, seperti yang disarankan oleh 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')

Catatan: kisi persegi ukuran 1 m dengan tingkat yang sama dengan lapisan yang berisi sekitar 1 juta poligon dan karenanya kinerja skrip berkurang secara masuk akal.

Antonio Falciano
sumber
1

Pertanyaan ini dijawab sejak waktu yang lalu, tetapi saya menambahkan solusi lain menggunakan pustaka rupawan dan 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
sumber
-1

Jawaban untuk Membuat Shapnet jala jaring di QGIS? menunjukkan opsi buat kisi di kotak alat pemrosesan QGIS.

pengguna965586
sumber
OP menyatakan bahwa dia lebih suka contoh dalam Python murni daripada dengan perangkat lunak
LaughU
Mengingat jawaban lain mengimpor pustaka, mengimpor modul QGIS adalah cara yang benar untuk menghindari GUI. Seperti dalam jawaban ini gis.stackexchange.com/questions/79916/…
user965586
Jawaban lain memberikan kode jadi jika Anda juga melakukannya maka saya pikir akan lebih baik diterima.
PolyGeo