I have code that I think you can modify, if all your rectangles are simple like above you need to get rotation, I would do this by average the azimuth of the short sides and long sides then find the nearest azimuth to 0 and then get the rotation. After that the logic below should work with some modification, I use this to find the NEWS sides of land sections.
I don't use QGIS so all these are arcpy methods, I only cursors and geometry so should be able to translate.
This returns azimuth for line segment.
def returnAzimuth(shape):
point1 = shape.firstPoint
point2 = shape.lastPoint
dX = point2.X-point1.X
dY = point2.Y-point1.Y
az = math.atan2(dX,dY)*180/math.pi
if az<0:
az = az+360
return az
return az
This function takes a feature class or shapefile and finds the width angle and length angle of the rectangle.
def getTractSideRatio(landpolyshape):
distances = []
area = landpolyshape.area
if area:
pnts= pnts = landpolyshape.getPart(0)
distance = 0
for x in range(pnts.count-2):
if distance==0: distance = returnInverse(pnts.getObject(x),pnts.getObject(x+1))[1]
az1=returnInverse(pnts.getObject(x),pnts.getObject(x+1))[0]
az2=returnInverse(pnts.getObject(x+1),pnts.getObject(x+2))[0]
length = returnInverse(pnts.getObject(x+1),pnts.getObject(x+2))[1]
topaz = returnInverse(pnts.getObject(x),pnts.getObject(x+1))[0]+2
lowaz = returnInverse(pnts.getObject(x),pnts.getObject(x+1))[0]-2
if az2 <= topaz and az2 >= lowaz:
distance = distance + length
if x==pnts.count-2 : distances.append(distance)
else:
distances.append(distance)
distance = 0
length = max(distances)
width = area/length
w2 = width/2
l2 = length/2
widthangle = 2*int(round(math.degrees(math.atan2(w2, l2))))
lengthangle = 2*int(round(math.degrees(math.atan2(l2, w2))))
return widthangle,lengthangle
This function return inverse from 2 points
def returnInverse(point1,point2):
dX = point2.X-point1.X
dY = point2.Y-point1.Y
dis = sqrt(dX**2+dY**2)
az = math.atan2(dX,dY)*180/math.pi
if az<0:
az = az+360
return az,dis
return az,dis
I have modified this to find the cardinal direction of the points of the polygon from the polygon center, you will need to provide the rotation angle for this to be accurate.
from math import sqrt
def getNEWSByAzimuth(landPoly,rot):
sr = arcpy.SpatialReference(6578)
widthangle,lengthangle = getTractSideRatio(landPoly)
if widthangle and lengthangle:
with arcpy.da.SearchCursor(landPoly.name,["SHAPE@"],spatial_reference=sr) as sc:
for row in sc:
polyShape=row[0]
pnts = polyShape.getPart(0)
for pnt in pnts:
polyCenter = polyShape.centroid
linetopoint = pnt
az=returnInverse(polyCenter,pnt)[0]
az = az + rot
if az>360:
az=az-360
if az < 360 and az >= 360 - lengthangle/2: #this sets the ratio for the rectangle lengthangle/2
side="NORTH"
print(side)
if az <= 0 + lengthangle/2 and az >= 0:
side="NORTH"
print(side)
if az < 270 + widthangle/2 and az >= 270 - widthangle/2:
side="WEST"
print(side)
if az < 180 + lengthangle/2 and az >= 180 - lengthangle/2:
side="SOUTH"
print(side)
if az > 90 - widthangle/2 and az < 90 + widthangle/2:
side="EAST"
print(side)