Geographic Information Systems Asked on October 3, 2020
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.
The only thing I can think of is to write a bit of code so that:
Am I over-complicating this? Is there an easier way? Or a tool out there that already exists?
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()
Correct answer by James on October 3, 2020
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:
Answered by Llaves on October 3, 2020
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP