1

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.enter image description here 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
2

How do you define NW corner? Does the image below have NW corner at all?

enter image description here

  • 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
1

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.

0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.