4

I am using the latest QGIS 3.4.4

I have a layer 1 which has a bunch of lines in all directions. I have another layer 2 intersecting those lines.

I want to delete all the lines in layer 1 that do not cross perpendicularly to layer 2.

Data: You can download the two line layers here

This seems simple but I cant seem to get it right.

Here is an example of the data I am working with

Improve this question

1 Answer 1

Reset to default
7

The Expression Way

I cannot write spatial SQL queries nor Python code yet. But there are a way to select perpendicular lines by expression. Then, deleting selected features is a trivial task.

There are, at least, two ways to determine if two 2D lines are perpendicular. One is comparing their azimuthes, another is comparig their slopes (in terms of growth ratio between y and x coordinates). Slopes have a problem that we should handle: a vertical line have infinite slope. Azimuths also have a problem: angles multiples of 2*k*PI() are the same azimuth, but it is easier to handle. Therefore, better if we use azimuths.

The azimuth entered in the expression is calculated as accurately as possible, so the lines must be exactly perpendicular. When testing this solution I found that the advanced digitization tool of QGIS draws lines perpendicular to one another that the expression does not recognize as such. When they are digitized at fixed angles it does, however. In case the solution does not meet your needs, calculate the azimuths in a field of the attribute tables, with the precision you need for the calculation, and refer to those fields in the expression (expressions will result even simpler, but less precise).

  • Two lines are perpendicular when their azimuthes differ by 90 degrees.

But the lines can be digitalized in two directions (forward or reversed), therefore:

  • Two lines are perpendicular when their azimuthes differ by -90 degrees, too.

Also, if adding an azimuth of a line 90 degrees, you get an angle greater than 360 degrees, you have to subtract 360. The same thing, if subtracting 90 you get a negative angle, you must add 360. Therefore:

  • Two lines are perpendicular when their azimuthes differ by +270 and by -270 degrees, too.

Some other preliminary considerations:

  • layer1 and layer2 have the same CRS.

  • In layer2 there is a field named id, and there is only one line with "id" = 1

Now, we can select the lines of layer1 that does not satisfy any of the conditions of perpendicularity against that line, with the following expression: 

( azimuth(  start_point( $geometry),  end_point( $geometry)) +  radians( 90) <> azimuth ( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))))
AND
( azimuth(  start_point( $geometry),  end_point( $geometry)) -  radians( 90) <> azimuth ( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))))
AND
( azimuth(  start_point( $geometry),  end_point( $geometry)) +  radians( 270) <> azimuth ( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))))
AND
( azimuth(  start_point( $geometry),  end_point( $geometry)) -  radians( 270) <> azimuth ( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))))

UPDATE

QGIS creating algorithms appear to be using 5 decimals of precision. I don't know if it depends on the reference system of the layer units. Test it please. Try with the following expression: 

( round ( azimuth( start_point( $geometry), end_point( $geometry)), 5) +  round( radians(90), 5) <> round ( azimuth( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))), 5))
AND
( round ( azimuth( start_point( $geometry), end_point( $geometry)), 5) -  round( radians(90), 5) <> round ( azimuth( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))), 5))
AND
( round ( azimuth( start_point( $geometry), end_point( $geometry)), 5) +  round( radians(270), 5) <> round( azimuth( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))), 5))
AND
( round ( azimuth( start_point( $geometry), end_point( $geometry)), 5) -  round( radians(270), 5) <> round ( azimuth( start_point( geometry(  get_feature( 'layer2', 'id', '1'))), end_point( geometry(  get_feature( 'layer2', 'id', '1')))), 5))
Improve this answer
6
  • 3
    Excellent! And it may become a little bit shorter if you use modulo, like (azimuth_layer1 - azimuth_Layer2) % radians(90) <> 0
    – Kazuhito
    Commented Jan 24, 2019 at 4:14
  • 1
    I am not sure what I am doing wrong, it seems to select all the lines in layer 1. The lines in layer 1 were created using the transect algorithm in qgis. I hope thats not the issue. Additionally is there anywhere I can maybe read up on modulo as I am new to this.
    – MapInfoNewbie
    Commented Jan 25, 2019 at 2:47
  • The modulo operator returns the remainder of a division. It is useful to determine if two numbers are multiples of each other (the remainder of their division is zero). In this case, we would have to handle the issue that the difference between the azimuth of two parallel lines is also a multiple of 90 degrees.
    – Gabriel De Luca
    Commented Jan 25, 2019 at 4:01
  • I updated my answer based on the tests on my layers (those I drew with the advanced digitization tool of QGIS), finding that the expression can exclude those lines if it is used by rounding all the values to 5 decimals. Test it, please.
    – Gabriel De Luca
    Commented Jan 25, 2019 at 4:15
  • Hi, I tried our solution but it did not work for me. It appears on the right track but maybe I cant get it to work. For reference, I have uploaded the two line layers and you can check yourself. Thanks for all your efforts.
    – MapInfoNewbie
    Commented Jan 29, 2019 at 3:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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