2

考虑一个有四个顶点的不规则多边形(P)。

我想将 P 细分为更小的多边形,如图所示。1.左上角多边形 P,右上角镶嵌两个,左下角镶嵌四个等

在过去,我使用 ArcGIS 9.3 和 Miles Hitchen 编写的 VBA 宏来实现这一点。VBA代码如下:

Option Explicit

Const xSplit As Long = 10
Const ySplit As Long = 5
Dim dCoords(xSplit, ySplit, 1) As Double


Public Sub GridSelectedQuadrilaterals()
Dim pMxDoc As IMxDocument
Dim pInFtrLyr As IFeatureLayer
Dim pFtrSel As IFeatureSelection
Dim pPolygon As IPolygon
Dim pEnumIDs As IEnumIDs
Dim lID As Long

    ' Get the first selected polygon on the first layer
    Set pMxDoc = ThisDocument
    Set pInFtrLyr = pMxDoc.FocusMap.Layer(0)
    Set pFtrSel = pInFtrLyr
    Set pEnumIDs = pFtrSel.SelectionSet.IDs
    pEnumIDs.Reset

    lID = pEnumIDs.Next
    While lID >= 0
        Set pPolygon = pInFtrLyr.FeatureClass.GetFeature(lID).Shape
        GridQuadrilateral pPolygon
        lID = pEnumIDs.Next
    Wend

    MsgBox "Finished"

End Sub


Private Sub GridQuadrilateral(pPolygon As IPolygon)
Dim pSegColl As ISegmentCollection
Dim lIdx As Long
Dim cx(3) As Double, cy(3) As Double
Dim dx(3) As Double, dy(3) As Double
Dim dx2 As Double, dy2 As Double
Dim x1 As Double, x2 As Double, x3 As Double
Dim y1 As Double, y2 As Double, y3 As Double
Dim l As Long, xs As Long, ys As Long

    Set pSegColl = pPolygon

    ' Get the corner coords of the quad
    lIdx = 0
    For l = 0 To 3
        lIdx = GetIndexOfNextCornerSegment(lIdx, pPolygon)
        cx(l) = pSegColl.Segment(lIdx).FromPoint.X
        cy(l) = pSegColl.Segment(lIdx).FromPoint.Y
    Next l

    dx(0) = (cx(1) - cx(0)) / xSplit
    dx(1) = (cx(1) - cx(2)) / ySplit
    dx(2) = (cx(2) - cx(3)) / xSplit
    dx(3) = (cx(0) - cx(3)) / ySplit

    dy(0) = (cy(1) - cy(0)) / xSplit
    dy(1) = (cy(1) - cy(2)) / ySplit
    dy(2) = (cy(2) - cy(3)) / xSplit
    dy(3) = (cy(0) - cy(3)) / ySplit

    For ys = 0 To ySplit
      x1 = cx(3) + dx(3) * ys
      y1 = cy(3) + dy(3) * ys
      x2 = cx(2) + dx(1) * ys
      y2 = cy(2) + dy(1) * ys
      dx2 = (x2 - x1) / xSplit
      dy2 = (y2 - y1) / xSplit

      For xs = 0 To xSplit
        x3 = x1 + dx2 * xs
        y3 = y1 + dy2 * xs
        dCoords(xs, ys, 0) = x3
        dCoords(xs, ys, 1) = y3
      Next xs

    Next ys

    BuildGrid

End Sub


Private Function GetIndexOfNextCornerSegment(lStartIdx As Long, pPolygon As IPolygon) As Long
Dim PI As Double
Dim pSegColl As ISegmentCollection
Dim pLine1 As ILine, pLine2 As ILine
Dim l As Long
Dim lNxtIdx As Long
Dim dAng As Double

    PI = Atn(1) * 4
    Set pSegColl = pPolygon
    For l = 0 To pSegColl.SegmentCount - 2
        lNxtIdx = lStartIdx + l
        If lNxtIdx = pSegColl.SegmentCount Then lNxtIdx = 0
        Set pLine1 = pSegColl.Segment(lNxtIdx)
        lNxtIdx = lNxtIdx + 1
        If lNxtIdx = pSegColl.SegmentCount Then lNxtIdx = 0
        Set pLine2 = pSegColl.Segment(lNxtIdx)
        dAng = Abs(pLine1.Angle - pLine2.Angle) * 180 / PI
        ' *** Ensure we only deal with angles upto 180 ***
        If dAng >= 180 Then dAng = 360 - dAng
        If dAng > 20 Then
            ' The start point of this segment is a corner point
            GetIndexOfNextCornerSegment = lNxtIdx
            Exit Function
        End If
    Next l

GetIndexOfNextCornerSegment = -1
End Function


Private Sub BuildGrid()
' Now create the polygons on 2nd layer
Dim pMxDoc As IMxDocument
Dim pOutFtrLyr As IFeatureLayer
Dim i As Long, j As Long
Dim pFtrCls As IFeatureClass
Dim pFtrCsr As IFeatureCursor
Dim pFtrBfr As IFeatureBuffer
Dim pPtColl As IPointCollection
Dim pPt As IPoint

    Set pMxDoc = ThisDocument
    Set pOutFtrLyr = pMxDoc.FocusMap.Layer(1)
    Set pFtrCls = pOutFtrLyr.FeatureClass
    Set pFtrBfr = pFtrCls.CreateFeatureBuffer
    Set pFtrCsr = pFtrCls.Insert(True)

    For i = 0 To ySplit - 1
        For j = 0 To xSplit - 1

            Set pPtColl = New Polygon
            Set pPt = New Point
            pPt.PutCoords dCoords(j, i, 0), dCoords(j, i, 1)
            pPtColl.AddPoint pPt
            pPt.PutCoords dCoords(j, i + 1, 0), dCoords(j, i + 1, 1)
            pPtColl.AddPoint pPt
            pPt.PutCoords dCoords(j + 1, i + 1, 0), dCoords(j + 1, i + 1, 1)
            pPtColl.AddPoint pPt
            pPt.PutCoords dCoords(j + 1, i, 0), dCoords(j + 1, i, 1)
            pPtColl.AddPoint pPt
            pPt.PutCoords dCoords(j, i, 0), dCoords(j, i, 1)
            pPtColl.AddPoint pPt

            Set pFtrBfr.Shape = pPtColl
            pFtrCsr.InsertFeature pFtrBfr

        Next j
    Next i

    pMxDoc.ActiveView.Refresh

End Sub

VBA 代码顶部的用户输入告诉脚本用户需要对每个轴进行多少细分

Const xSplit As Long = 10
    Const ySplit As Long = 5

谁能帮我用 Python 翻译这段代码?我想使用 Python 和 OGR 来实现同样的目标。

4

1 回答 1

0

编辑(部分 - 已解决):以下 Python 代码解决了部分问题。进行了细分,但仅将其中一个子多边形的结果写入新的形状文件。在 Python 中完成此算法的任何其他帮助将不胜感激。

import os
import ogr
import sys
import numpy as np

#Driver
driver = ogr.GetDriverByName('ESRI Shapefile')

#Working Directory Path
ws = 'C:/data/'

#Make ws the working directory
os.chdir(ws)

#vector path
in_data = 'data.shp'

#user inputs; how many splits per axis
xSplit = 2
ySplit = 2

#This is a three dimensional array (In numpy Z,Y,X values) 
dCoords = np.ones((2,ySplit+1,xSplit+1))

#Input polygons
source_ds = ogr.Open(in_data,0)
source_layer = source_ds.GetLayer(0)
source_features = source_layer.GetFeatureCount()

#Output polygons
target_ds = driver.CreateDataSource('split.shp')

if target_ds is None:
    print 'Creation of output file failed.\n'
    sys.exit(1)

geom_type = ogr.wkbPolygon
target_layer = target_ds.CreateLayer('split',None,geom_type)

if target_layer is None:
    print 'Creating of layer failed.\n'
    sys.exit(1)

field_defn = ogr.FieldDefn('parcel_no',ogr.OFTInteger)
field_defn.SetWidth(30)

if target_layer.CreateField(field_defn) !=0:
    print 'Creating plot_num field failed.\n'
    sys.exit(1)

featureDefn = target_layer.GetLayerDefn()

#Get the polygons from the layer
for i in xrange(0,source_features):
    source_poly = source_layer.GetFeature(i)
    source_plot = source_poly.GetField('id')
    source_geom = source_poly.GetGeometryRef()

    #Get the corner coords of the quad
    cx = []
    cy = []

    pts = source_geom.GetGeometryRef(0)
    points = []

    for p in xrange(pts.GetPointCount()):
        cx.append(pts.GetX(p))
        cy.append(pts.GetY(p))

    dx0 = (cx[1] - cx[0]) / xSplit
    dx1 = (cx[1] - cx[2]) / ySplit
    dx2 = (cx[2] - cx[3]) / xSplit
    dx3 = (cx[0] - cx[3]) / ySplit
    dy0 = (cy[1] - cy[0]) / xSplit
    dy1 = (cy[1] - cy[2]) / ySplit
    dy2 = (cy[2] - cy[3]) / xSplit
    dy3 = (cy[0] - cy[3]) / ySplit

    for ys in xrange(0,ySplit):
        x1 = cx[3] + dx3 * ys
        y1 = cy[3] + dy3 * ys
        x2 = cx[2] + dx1 * ys
        y2 = cy[2] + dy1 * ys
        dx2 = ( x2 - x1 )  / xSplit
        dy2 = ( y2 - y1 )  / xSplit

        for xs in xrange(0,xSplit):
            x3 = x1 + dx2 * xs
            y3 = y1 + dy2 * xs

            #Calculate Grid Coordinates
            dCoords[0][xs][ys] = x3
            dCoords[1][xs][ys] = y3

    #Build Grid
    for i in xrange(0,ySplit - 1):
        for j in xrange(0,xSplit - 1):

            #Create an empty ring geometry
            ring = ogr.Geometry(ogr.wkbLinearRing)

            ring.AddPoint(dCoords[0][j][i],dCoords[1][j][i])
            ring.AddPoint(dCoords[0][j][i+1],dCoords[1][j][i+1])
            ring.AddPoint(dCoords[0][j + 1][i + 1],dCoords[1][j + 1][i + 1])
            ring.AddPoint(dCoords[0][j + 1][i],dCoords[1][j + 1][i])
            ring.AddPoint(dCoords[0][j][i], dCoords[1][j][i])

            target_poly = ogr.Geometry(ogr.wkbPolygon)
            target_poly.AddGeometry(ring)

            target_feature = ogr.Feature(featureDefn)
            target_feature.SetGeometry(target_poly)
            target_feature.SetField('id',source_plot)
            target_layer.CreateFeature(target_feature)
于 2013-05-03T09:05:01.887 回答