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 - (a - lat) * (a - b) / (a - b) 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) for e in exact: longitudes.append(nodes[e]) # 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) 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)