1

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.

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.

enter image description here

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

enter image description here

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.

enter image description here

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

enter image description here

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

  • Hi, I was able to get the shapefile nodes coordinates (Vector\Geometry Tools\Extract Nodes... and Vector\Geometry Tools\Export/Add Geometry Columns...) but I am stuck there. I have tried Select Feature using an expression and the Query Builder in the Properties but to no avail... could you expand on how to use PostGIS ? – YeO Jul 6 '16 at 12:49
  • If the data is in the correct CRS for example you are looking for lat/lon and your data is in lat/lon, then what you have done so far is 95% the way there. See the workflow in my amended answer you sound very close. – Andrew Jeffrey Jul 6 '16 at 21:36
  • Nice amendment, thanks. I was there actually but most of the time, there will be no exact match (crossing in between two nodes) so I went for a python script to compute these intersections – YeO Jul 6 '16 at 22:07
1

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.