Get the coordinates of all points having the same given latitude along a line shapefile

Question

I have a shapefile representing a shoreline and I need to extract the coordinates of the points at a certain latitude (or longitude).

Is there a tool or plugin that could do that ?

Edit: to be more precise, I need to compute/interpolate the intersections where the lines actually crosses the polyline between 2 nodes not just locate the nearest node.

Answer 1

You could calculate lat and lon as attributes in you shapefile and then use the "select features using an expression" tool. Important to make sure that the lat and lon are using the same precision and scale as the coordinate that you are searching for.

Amended - included steps for using the "Select features using an expression tool"

You mentioned in your comment that you have extracted nodes and then calculated the geometry. So I am assuming you point attribute table looks like mine below with a "XCOORD" and "YCOORD" column.

You then want to click on the "Select features using an expression tool button" shown below in the red circle.

This opens a dialog where you can create an expression using your attribute columns and values. For example below I am going to query the "XCOORD" attribute for all values that are equal to 146.587500, hit the select button at the bottom right to execute the expression.

The result will be a selection in your attribute table which you can export to a new layer if required. My resultant selection returned two points out of 147 with the same "XCOORD".

With the method make sure that the value you are searching for in the expression tool is an exact match to the value you want to find in your table.

Alternatively if the shapefile is of a large shoreline consider putting it in PostGIS. This could speed things up for you, but probably only worth doing if it something that you'll do more than once. If you have the nodes of the shoreline you could then do a query on the geometry for all values that equal a certain latitude. For example

SELECT * FROM table WHERE ST_Y(the_geom)::decimal(5,3) = -35.883

Answer 2

After much trial and error, I went for a Python script to solve my problem as I want to get the Longitudes of the crossings in between nodes rather than getting the closest node. In retrospect, I should have been more clear about that in my question. Here's my code:

from qgis.core import *
from qgis.gui import *
import qgis.utils
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
            nodes = np.asarray(geom.asPolyline())
            # 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))
            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)