Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
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)
How do you define NW corner? Does the image below have NW corner at all?
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))[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)
