Find the NW vertex of a rotated polygon

I am using QGIS 2.18 and I am wondering how to find the NW corner of a polygon and create an attribute of that point number. See the attached image. For example, please node 8 is the NW corner for Grid #1, 13 is the NW corner for Grid 5 and so on. any help with this would be wonderful. The two files I am working with is

Grid = 'grid' and

Vertex = 'nodes' (Bold)

• Where should I see nodes 8 and 13? – user30184 Sep 24 '18 at 20:59
• to the left. but the name of the node is not all that important – lowsparked Sep 24 '18 at 21:06
• Not so much but I would like to be sure where it is. – user30184 Sep 24 '18 at 21:14
• nodes are labelled east to west, Thaks for asking for clarification, every little bit helps is asking a question. – lowsparked Sep 24 '18 at 21:18

How do you define NW corner? Does the image below have NW corner at all? • 2 would be the NW corner, draw a line down the middle of the polygon east to west and what is the most NW corner. for this one 2 and 3 are the two most north points and 2 is the further west most point – lowsparked Sep 24 '18 at 21:48
• Call your points differently. Hard to imagine point 2 being north west. +1 for answer. – FelixIP Sep 25 '18 at 3:20
• not sure what you mean by calling your points differently? – lowsparked Sep 25 '18 at 5:06
• I agree with @FelixIP. NW is very well known compass direction en.wikipedia.org/wiki/File:Compass_Rose_English_North.svg and it points to 315° N. What you are after is something else. – user30184 Sep 25 '18 at 6:14
• @lowsparked according to what you proposed here, you've got to implement a rule, as results will most likely be inconsistend for different border cases; get the two northern-most points and the western-most of those could be one, another could be find the NW point of the bbox and get the nearest vertex to that point...which would yield different, but consistent and directionally weighted results. – ThingumaBob Sep 25 '18 at 8:41

Having seen ongoing discussions about how to define "NW", I am not sure if this is valid answer. The idea here is to calculate the azimuth between the nodes and the polygon centroid, then to find the most closest to NW (315 degrees) direction.

QGIS 2.18.20 (probably 2.18.21, too)

SELECT grid.*,
nodes.Node_Id,
abs(degrees(st_azimuth(st_centroid(grid.geometry), nodes.geometry))- 315) AS D
FROM grid, nodes
WHERE st_touches(grid.geometry, nodes.geometry) = 1
GROUP BY grid.id
ORDER BY abs(degrees(st_azimuth(st_centroid(grid.geometry), nodes.geometry))- 315) ASC
LIMIT 1

If you want to find NE corner, change 315 to 45.

NB. This does not give solution when the rectangle lies lateral to NE line (45-deg) and corner nodes are located at equi-distance.

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))
az1=returnInverse(pnts.getObject(x),pnts.getObject(x+1))
az2=returnInverse(pnts.getObject(x+1),pnts.getObject(x+2))
length = returnInverse(pnts.getObject(x+1),pnts.getObject(x+2))
topaz = returnInverse(pnts.getObject(x),pnts.getObject(x+1))+2
lowaz = returnInverse(pnts.getObject(x),pnts.getObject(x+1))-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
pnts = polyShape.getPart(0)
for pnt in pnts:
polyCenter = polyShape.centroid
linetopoint = pnt
az=returnInverse(polyCenter,pnt)
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)