2
Рассмотрим нерегулярный многоугольник (P), который имеет четыре вершины.Вычисление прямоугольных многоугольников в равных частях с использованием Python/OGR
Я хотел бы разделить P на более мелкие многоугольники, как показано на рис. 1. 
В старые времена я использовал ArcGIS 9.3 и макрос 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.