1

https://onedrive.live.com/embed?cid=AEC7953479B58DC0&resid=AEC7953479B58DC0%218329&authkey=AKHlm_SovPNQOYw

I have a series of points (ObservationPoints.shp) along a line (Centreline.shp). I am looking to find the longest line from these ObservationPoints to the Centreline without crossing the boundary of a polygon (PolygonBoundary.shp).

The image below shows 3 examples, the first two are correct examples (cyan - polygon, magenta - points and dark blue - centreline) where the longest lines (LongestLines.shp) do not cross the polygon boundary. The third example is incorrect as the longest line crosses the polygon boundary.

Examples

The only thing I can think of is to write a bit of code so that:

  • for every point, draw radial lines at a given bearing spacing and length then,
  • clip those lines to the polygon boundary then,
  • remove all but the longest of those radial lines then,
  • clip this longest line to the centreline.

Am I over-complicating this? Is there an easier way? Or a tool out there that already exists?

  • Could you please clarify your question? Your title says point within a polygon, but your three examples all involve points on the boundary. Are asking how to find the two most distant points on the boundary for which the segment connecting them lies within the polygon? Also, what's the difference between your second and third examples? Both cross the boundary, unless you are treating the line thickness as having meaning. – Llaves Jul 23 at 3:37
  • Thanks Llaves, I have edited/updated the post with link to .shp files- hopefully it is clearer. – James Jul 23 at 10:18
  • @BERA - all of the points (the point are 50m spacing along the line), so I can get multiple 'longest lines'. It would be like a line of sight assessment. – James Jul 23 at 12:16
  • Your OneDrive link yields "won't load now". Maybe you haven't set sharing permissions or something? (The "won't load" refers to browser - I can't even download to look at them) – Llaves Jul 23 at 16:51
  • 1
    I'm fascinated. What is the background for this analysis? – Simbamangu Jul 25 at 5:02
1

Your proposed procedure would work if you are willing to accept an approximation, since there is no guarantee that your one of your radial bearing lines is the true longest line. It's also horribly inefficient.

A slightly better procedure would be:

  1. For each observation point Pi, draw the line to Pi + 1, Pi+2 until the line crosses the boundary. Let's say the last point where the line remains inside the polygon is P i + n
  2. The longest line will intersect the centerline between points i+n and i+n+1. Now use binary search to create new points along the centerline, drawing a line from Pi until you reach the precision you desire.
| improve this answer | |
  • Thanks Llave, this seems a LOT more efficient than the method I proposed. I am curious - could I start off with a binary search (draw line from Pi to last point, 1/2 way, 1/4 way until it doesn't crosses the boundary? Would this be more efficient? or I suppose it is based on the geometry. Thanks again, much appreciated. – James Jul 27 at 8:11
  • If the curve is "reasonably smooth", then for point P_i+1 I'd start the search at the last point P_i could "see". For polygons where the visibility is long and the curvature is small, this could save a lot of work – Llaves Aug 3 at 18:12
0

See my finalised code.

BinarySearch - Visibility

name1='Interpolated points' #QGIS Layer - the closer the spacing the greater the accuracy
name2='VergeExtents' #QGIS Layer - Extents to which the visibility should not cross

layer1 = QgsProject.instance().mapLayersByName(name1)[0]
layer2= QgsProject.instance().mapLayersByName(name2)[0]

#initialise variables
pi=0 #first initial observation point
pn=((layer1.featureCount())) #end initial observation point
 
for feature in layer2.getFeatures():
    geom3=feature.geometry()

#create the longest line vector layer
v_layer = QgsVectorLayer('LineString?crs=epsg:27700', 'HMVLines', 'memory')
pr1 = v_layer.dataProvider()

#For loop finding the longest line for every points at a given interval
for x in range(1,pn,50):
    pi=x 
    l=pi+1 #low point
    h=pn # high point


#While loop to binary search the longest line from a given start point from the for loop
while l < h: 
    #initialise start of line
    for feature in layer1.getFeatures([pi]):
        geom1=feature.geometry()
    
    #iteration on the midpoint of the point set
    mid = l + (h - l) // 2; 
    
    #initialise end of line
    for feature in layer1.getFeatures([mid]):
        geom2=feature.geometry()
   
    #create line between pi and mid
    line=QgsFeature()
    line.setGeometry(QgsGeometry.fromPolylineXY([geom1.asPoint(), geom2.asPoint()]))
              
    #1ST CHECK - if difference between h and l is less 1 then stop while loop(i.e. minimum spacing between points)
    if (h-l)<=1: 
        break
    
    #2ND CHECK - if the visibility line and the boundary intersect, set high point to middle
    elif geom3.intersects(line.geometry())==1:
        h=mid
       
    #3RD CHECK - if the visibility line and the boundary do not intersect, set low point to middle
    elif geom3.intersects(line.geometry())==0: 
        l=mid
    
#Once binary search has complete, l will be the longest line that does not intersects the boundary - get this geometry.
for feature in layer1.getFeatures([l]):
    geom4=(feature.geometry())
        
#create the longest line from pi to l
visibility = QgsFeature()
visibility.setGeometry(QgsGeometry.fromPolylineXY([geom1.asPoint(), geom4.asPoint()]))
pr1.addFeatures([visibility])

QgsProject.instance().addMapLayers([v_layer])
iface.zoomToActiveLayer()
| improve this answer | |

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.