5

I have a point layer (xyz) and want to find all points that have elevation difference of 200m or more. I would like to then draw a line between each pair with the largest elevation difference.

The goal is to identify areas that have large elevation difference and (eventually) nearest distance between them in QGIS.

8

The easiest way is to import your points into a format that can be queried with SQL, like PostGIS, SQLite or Shapefile (using OGR).

Then you can query:

SELECT * 
FROM [table] a, [table] b 
WHERE a.[featureid] <> b.[featureid] 
    AND ST_Z(a.[geometry]) - ST_Z(b.[geometry]) >= 200;

Or you can query and make lines in one step:

SELECT ST_MakeLine(a.[geometry], b.[geometry]) AS geom 
FROM [table] a, [table] b 
WHERE a.[featureid] <> b.[featureid] 
    AND ST_Z(a.[geometry]) - ST_Z(b.[geometry]) >= 200;

It might help to filter out any points that cannot be >= 200m from another point first, because this is an O(n^2) query -- if you have a lot of points, it'll take a long time. You can also restrict the matching to within a certain distance to employ the index and speed it up.

Note: replace the [names] here with your own table/column names.

Edit:

I couldn't find where to execute an SQL query in QGIS, but here's how to do it in ogr2ogr on the command line:

ogr2ogr \
    -f "ESRI Shapefile" \
    -dialect "SQLITE" \
    -sql "SELECT * FROM input a, input b WHERE ST_Z(a.GEOMETRY)-ST_Z(b.GEOMETRY) < 200 AND a.FID<>b.FID" \
    output.shp input.shp

The -f switch tells the program what kind of file to produce. The -dialect switch says to use SQLite's SQL engine to perform the query. The -sql argument is the query. GEOMETRY and FID are built-in field names that you can use in your query.

Now, this query is going to be extremely slow because there's no index. It's going to try to match every point to every other point. Using the documentation here I tried to get the query to use a spatial index, but it doesn't seem to work and I don't have time to mess with it at the moment.

Loading the same point set into SpatiaLite goes something like this:

spatialite> .loadshp input input UTF-8 2956
spatialite> select createspatialindex('input', 'geometry');
spatialite> create table output as select * from input a, input b where st_z(a.geometry)-st_z(b.geometry) < 200 and a.rowid<>b.rowid and a.rowid in (select rowid from spatialindex f_table_name='stations' and search_frame=b.geometry);
spatialite> .dumpshp output 'geometry' output UTF-8

The first line loads the shp file, the second creates a spatial index, the third creates a new table from the query and the fourth dumps the new table to a new file. (Docs here.)

This is much faster than the ogr2ogr solution (for me).

  • Thank you for your help. Please bare with me and my beginner questions: – PeterB Aug 4 '15 at 22:14
  • I converted the point layer to a shapefile, are you suggesting to query it in QGIS' "select by expression" function or query it in some other SQL capable program? Or is there a SQL plugin I need for QGIS? Right now, all elevation points are in one table - are you suggesting to duplicate them into table a and table b? Also, for points, what [geometry] should I use? – PeterB Aug 4 '15 at 22:24
  • 1
    If your data were in postgis or spatialite, the right place to query it would be in databases > db manager. – Alexandre Neto Aug 5 '15 at 7:55
  • 1
    Actually, you don't duplicate anything. That just tells the SQL interpeter that you're doing a cross join, so you're comparing input to itself. (This syntax is actually functionally equivalent to input a CROSS JOIN input b.) The a and b are just aliases so you can refer to the sides of the join unambiguously. – Rob Skelly Aug 7 '15 at 16:12
  • 1
    To clarify further, that syntax input a is a bit confusing. It's not "here's an input called 'a'", it's "here's an input called 'input' with the alias 'a'". You could also write input AS a, which would be clearer. – Rob Skelly Aug 7 '15 at 16:19

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.