I wrote a script to compute / interpolate the nodes of a Polyline at a given latitude, input_lat
. It works for my purpose but I am wondering if there is a better, more optimized way of doing this.
The idea is to first check if there is any exact match, i.e. if there are nodes at the exact given latitude so I subtract the input latitude to the list of nodes (the latitudes
list) and if there are any zeros, I can just read the longitude of these.
Most of the time the crossing will be between 2 consecutive nodes, with one being below (or above) and the next one being on the other side of the input_lat
(the node above will retain a positive value after subtracting input_lat
and the node below a negative value) so in my latitudes
list I detect that by checking for a change of sign between two consecutive elements, so I multiply each element with the next: if the result is negative, there is a change of sign.
I do this with np.multiply(latitudes[1:], latitudes[:-1])
but I wonder if there is a better way.
The script is meant to run from within QGIS so all the qgis.*
libraries that are imported and the Qgs
commands are specific to QGIS.
"""Script to compute coordinates of a Polyline at a given Latitude"""
#import QGIS specific libraries
from qgis.core import *
from qgis.gui import *
import qgis.utils
#import numpy
import numpy as np
"""Assuming the following :
- a Line/Polyline Layer is selected"""
#the latitude we want to compute intersections with
input_lat = 28.54456111
ray = [0, input_lat]
longitudes = []
def getIntersectLon(a, b, lat):
"""Compute the Longitude of the intersection between two nodes"""
return a[0] - (a[1] - lat) * (a[0] - b[0]) / (a[1] - b[1])
layer = qgis.utils.iface.activeLayer()
if layer is None:
iface.messageBar().pushMessage("Select a Layer", level=QgsMessageBar.CRITICAL, duration=2)
else:
for feature in layer.getFeatures():
geom = feature.geometry()
if geom.type() == QGis.Line:
# get the nodes coordinates in an array :
# geom.asPolyline() will produce an array of tuples
# that we will convert to an array of arrays so we
# can mutate the values
nodes = np.asarray(geom.asPolyline())
#
# THIS IS THE CORE ROUTINE
#
# get the list of the difference between latitudes and the input
latitudes = (nodes - ray)[:,1]
# if there are zeros, we have points at the exact given latitude
exact = np.where(latitudes == 0)[0]
for e in exact:
longitudes.append(nodes[e][0])
# where line crosses between nodes, there will be a change of sign
# one node being below, and the next above the input
xing = np.multiply(latitudes[1:], latitudes[:-1])
# get the indexes of the sign changes
crossing = np.where(xing < 0)[0]
for c in crossing:
longitudes.append(getIntersectLon(nodes[c], nodes[c+1], input_lat))
#
# THIS IS THE END OF THE CORE ROUTINE
#
# we will now create points the found latitudes (if any)
# that we will load into a new Point layer for display
#
if longitudes:
#create a Point layer to store the intersections found
newlayer = QgsVectorLayer("Point", "Intersections", "memory")
pr = newlayer.dataProvider()
feat = QgsFeature()
for lon in longitudes:
feat.setGeometry(QgsGeometry.fromPoint(QgsPoint(lon, input_lat)))
pr.addFeatures( [ feat ] )
print lon
QgsMapLayerRegistry.instance().addMapLayers([newlayer])
qgis.utils.iface.mapCanvas().refresh()
else:
iface.messageBar().pushMessage("Select a Line layer", level=QgsMessageBar.CRITICAL, duration=2)
qgis.*
libraries are for the script to run within QGIS so the only one that matters isnumpy
really. Not sure I understand what you are saying: isn't that what I am doing already (substracting the latitude from the list of latitudes from the Polyline and then looking for zeros/roots) ? \$\endgroup\$ – YeO Jul 7 '16 at 13:12latitudes = (nodes - ray)[:,1]
why not just usenp.roots(latitudes)
? \$\endgroup\$ – Oscar Smith Jul 7 '16 at 13:16