Computing intersections of a Polyline and a line 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 &lt; 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)
In getIntersectLon, I didn't find the docstring very helpful. For y = m*x + b, I think of m as rise over run, but in (a[0] - b[0]) / (a[1] - b[1]) I think I'm seeing run over rise. Ok. So we're multiplying 1/m by some delta_latitude, then we subtract, which seems like it could be related to the docstring. Maybe the docstring should talk about how the line described by two points intersects a line of latitude. Given that (lat,lng) are spherical coordinates, the approximation seems to have accuracy that falls off with increasing distance. Maybe this function should be named getLongitudeAt, with 1st parameter being the line of latitude.

It's not obvious to me that a[1] - b[1] shall be non-zero.

The latitudes == 0 test is maybe more stringent than you would like. Consider using abs(latitudes) &lt; epsilon.

If you're mostly interested in sign changes, you might want to define sgn() and then work with small signum values.

def sgn(n):
    return math.copysign(1, n)


The comments "THIS IS THE {|END OF} THE CORE ROUTINE" suggest that you may want to break out a helper routine.