Computing intersections of a Polyline and a line

Question

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)

Answer 1

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) < 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.

Answer 2

            # 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])

If the polyline is made up of great circle arcs (i.e. "straight" lines) then this isn't true (unless input_lat == 0, because the equator is the only line of constant latitude which is a great circle). An extreme example to make the point: consider a great circle arc from 10N 90W to 10N 90E. Both endpoints are near the equator, but the line between them passes through the North Pole.

I'm surprised that QGis doesn't have libraries which take care of this for you, but if you really have to do it yourself then it's easier to work with the embedding into 3D Euclidean geometry. To first order we model the Earth as a sphere, so a great circle is the intersection of a sphere with a plane through the centre, and a line of constant latitude is the intersection of a sphere with a normal through a pole. The intersection of two planes is a line (although note the special case that the planes don't intersect), and the intersection of a line with a sphere is 0, 1, or 2 points, (although obvious they might not be in the right range for the great circle arc).

Ignoring issues of numerical analysis, if the endpoints are \$p_1 = (\theta_1, \phi_1) = (x_1, y_1, z_1)\$ (respectively in spherical and Euclidean coordinates), and \$p_2\$ similarly, the plane of their great circle is the locus \$q \cdot (p_1 \times p_2) = 0\$. The plane through latitude \$\lambda\$ is \$q \cdot (0, 0, 1) = \sin \lambda \$. If we take the route of saying that \$q = (\lambda, \mu) = (\cos \mu \cos \lambda, \sin \mu \cos\lambda, \sin\lambda)\$ and let \$p_1 \times p_2 = (x_N, y_N, z_N)\$ then we have \$x_N \cos \mu + y_N \sin \mu = - z_N \tan\lambda\$, which can be solved for \$\mu\$. Then to check whether the resulting values (plural!) are in the great circle arc we calculate the midpoint \$m = \frac{p_1 + p_2}{2}\$ (we don't even need to normalise) and verify that \$q \cdot m \ge p_1 \cdot m\$.