# Convex Quadrilateral Polygons subdivision in equal parts using Python / OGR

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

I would like to subdivide P into smaller polygons as shown here http://i.stack.imgur.com/4SKhB.jpg

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

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
lID = pEnumIDs.Next
Wend

MsgBox "Finished"

End Sub

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)
pPt.PutCoords dCoords(j, i + 1, 0), dCoords(j, i + 1, 1)
pPt.PutCoords dCoords(j + 1, i + 1, 0), dCoords(j + 1, i + 1, 1)
pPt.PutCoords dCoords(j + 1, i, 0), dCoords(j + 1, i, 1)
pPt.PutCoords dCoords(j, i, 0), dCoords(j, i, 1)

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.

-
Just make sure to apply this code to convex quadrilaterals only. (With concave ones it will produce junk.) – whuber Mar 11 '13 at 18:01
Yes, you are absolutely right. Thank you for adding that piece of information. – Dimitris Mar 12 '13 at 9:20
@Dimitris, please add your "edits" as an answer and then mark the answer as "accepted", as opposed to editing the title as solved. Thanks. – RyanDalton Mar 24 '13 at 15:55
Done. Thanks for the advice. – Dimitris Mar 25 '13 at 9:41

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.

``````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 + 1][i + 1],dCoords[1][j + 1][i + 1])