Convex Quadrilateral Polygons subdivision in equal parts using Python / OGR

Question

Consider an irregular polygon (P) which has four vertices.

I would like to subdivide P into smaller polygons as shown in fig. 1.

In the old days I was using ArcGIS 9.3 and a VBA macro written by Miles Hitchen to achieve this. The VBA code is following:

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

The user input at the top of the VBA code tells the script how many tessellations the user requires for each axis

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

Could anyone help me translate this code in Python? I want to achieve the same thing using Python and OGR.

Answer 1

EDIT (PARTLY - SOLVED): Python code that is following solves part of the problem. The subdivision is made but the result is written for only one of the sub polygons to the new shape file. Any additional help to finish this algorithm in Python would be greatly appreciated.

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)