1

I have a river shapefile, got from r.watershed module, and it looks like this:

enter image description here

I have eight points on the river which I want to analyse, and I only want this river segments:(the second photo below is the river I want)

enter image description here

enter image description here

So I open its attribute table to delete segments I don't want, but I think it's not efficient.

I have tried the intersect tool in Geoprocessing tools, but the results were not what I want.

Is there more efficient way to achieve this goal?

  • 1
    To be clear, what you want is, the river segment next to the point? – huckfinn Jan 29 '14 at 13:03
  • 1
    Have you tried Select By Location to get the segments that intersect the points then saving the selection as a layer or etc. – landocalrissian Jan 29 '14 at 13:09
  • 1
    Check attribute table. "Hydro" tools usually offer some form of identyfication cerain features (ie. rivers chould have catchment IDs). Select all features within certain catchment (same ID). – Tomek Jan 29 '14 at 13:17
  • @huckfinn I've updated my post, thank you for helping me make my post clear. – Heinz Jan 29 '14 at 13:42
  • @ChrisR I've tried Select By Location but got nothing because the points are just beside the line segment and not intersect with it. – Heinz Jan 29 '14 at 13:43
3

Select features that are within a distance of the points rather than intersect (as not all of them will directly intersect). I'm not sure of the scale of your files, so try within 10m at first and if that doesn't work try a larger amount. Then flip the selection and delete the items.

  • thank you, but where is the tool that I can select features within a distance? – Heinz Jan 29 '14 at 13:38
  • 3
    In QGIS you can do this by buffering the points, then selecting features from the river that intersect the buffer – Alan Carr Jan 29 '14 at 14:00
3

If you can use Postgresql/PostGIS I've this code snippet for you. May be this will help you. The strategy is calculate the distance between each river vertex and the station point (adjacence list) and select the two naerest to build the corresponding segment. I think a transcription in Python will be easy.

Preparation of river data set:

    -- create a river table 
    CREATE TABLE river (id SERIAL);

    -- add a PostGIS LINESTRING column
    SELECT addGeometryColumn('public','river','geom',4326,'LINESTRING',2);

    -- bring in a river ...
    INSERT INTO  river(geom) 
    VALUES (ST_GeometryFromText(
    'LINESTRING(
      13.6938847 53.8612068,
      13.6974361 53.8628047,
      13.7013093 53.8634454,
      13.7038579 53.8632003,
      13.7058647 53.8627565,
      13.7085524 53.8616235,
      13.7102725 53.8610502,
      13.7127681 53.8602984,
      13.7166154 53.8596167,
      13.7185271 53.8596682,
      13.7216061 53.8607707)',
      4326));

Preparation of station data set:

    -- create a table for the river_stations
    CREATE TABLE river_station (id SERIAL);

    -- add a PostGIS POINT column
    SELECT addGeometryColumn('public','river_station','geom',4326,'POINT',2);

    -- Set some pointes near the river 
    INSERT INTO river_station(geom) VALUES
    (ST_GeometryFromText('POINT(13.70 53.85)',4326)),
    (ST_GeometryFromText('POINT(13.71 53.86)',4326));

The function that will help to find the segment.

    -- Function to find the neares segment of a STATION type POINT 
    -- to a LINESTRING hold in table river 
    -- as gives you the corresponding WKT LINESTRING
    CREATE OR REPLACE FUNCTION nearest_segment (GEOMETRY)
    RETURNS TEXT AS $$
    DECLARE result   TEXT;
       pos      INTEGER;
       rec      RECORD;
       station  GEOMETRY;
       river_id INTEGER;
       posX     TEXT;
       posY     TEXT;
    BEGIN

     -- Station Geometry parameter 
     station  := $1;

     -- Flag for adding the comma in WKT
     pos := 0;

     -- LINSTRING prefix in result
     result ='LINESTRING(';   

     -- Find the two closest points of the river goemetry to the station
     FOR rec IN SELECT
      (ST_DumpPoints(geom)).geom as point,
       st_distance((ST_DumpPoints(geom)).geom, station) as distance
     FROM river 
     -- WHERE id = $2 if you want to a special select a river object 
     ORDER BY distance ASCI LIMIT 2 

     LOOP 

      -- FLOAT8 to TEXT conversation of the segement part
        posX := st_X(rec.point)::TEXT; 
        posY := st_Y(rec.point)::TEXT;

     -- BUILD the WKT STRING 
        IF pos = 0 THEN
          result := result||posX||' '||posY;
        ELSE 
          result := result||', '||posX||' '||posY;
        END IF;   

       -- set mark to use commas before
         pos := pos + 1;
      END LOOP;

    -- give back the WKT String, as an alternative a geometry
     RETURN result||')';
   END;
   $$  LANGUAGE plpgsql;

Query the segments:

   -- Find the segments for the river 
   SELECT 
     st_astext(geom) as station, 
     nearest_segment(geom) as segment
    FROM
      river_station;

The result:

           station      |                      segment                       
   --------------------+----------------------------------------------------------
     POINT(13.7 53.85)  | LINESTRING(13.6938847 53.8612068, 13.6974361 53.8628047)
     POINT(13.71 53.86) | LINESTRING(13.7102725 53.8610502, 13.7085524 53.8616235)

May be there will be problems, if two segments of different rivers are very close to the station. Then you will need a parameter (a river_id in the station ensemble) to resolve the ambiguous situation.

To add the ID's from the river_stations you can use SQL Statements:

Prepare the table to hold the segments:

     CREATE TABLE river_segments (id SERIAL, ref_station INTEGER);
     SELECT addGeometryColumn('public','river_segments','geom',4326,'LINESTRING',2);

Insert the data sets:

     INSERT INTO river_segments (ref_station, geom)
        SELECT 
          id as ref_station, 
          ST_GeometryFromText(nearest_segment(geom),4326) AS geom 
     FROM river_station;

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.