考虑一个有四个顶点的不规则多边形(P)。
我想将 P 细分为更小的多边形,如图所示。1.
在过去,我使用 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 来实现同样的目标。