1

In QGIS 3.10.4, I have two shapefiles (polyline and polygon, layer name 'boundaries', 'lines' respectively), look like below:

enter image description here

My goal is to calculate the perpendicular distance from the centroid of each region to each of the lines. I've tried the two methods in this thread: Creating dummy variables between layers based on their geometries in QGIS, which are: using NNJoin, or using query in a Virtual Layer.

However, I met some problems with both of the methods:

  1. NNJoin (I selected the option "Approximate geometries by centroids"): the distances are in decimal degrees because the join layer (polyline layer) uses Geographic CRS. It's hard to interpret the distances and I wasn't able to figure out how to convert Geographic CRS to Projected CRS. Any idea how can I solve this problem?
  2. Then I tried the following query (provided by @Taras) in Add Layer -> Add/Edit Virtual Layer, it successfully generates a Virtual Layer with a "distance" column, with one distance associated with one region:
SELECT "boundaries".*,
        ROUND(ST_Length(ST_ShortestLine("boundaries".geometry, "lines".geometry)),2) AS distance
FROM "boundaries", "lines"
GROUP BY "boundaries".id
ORDER BY MIN(ST_Length(ST_ShortestLine("boundaries".geometry, "lines".geometry)))

So my ultimate question is: how can I calculate the perpendicular distance (in km or mile) from the centroids of each region on the polygon to each line of the polyline?

8
  • 1
    What should be used when it is impossible to draw a perpendicular line between the centroid and the boundaries? (ex: lower left regions with lines pointing to the north-east) – JGH Nov 24 '20 at 20:53
  • Using the SQL, to get distances in meters, use ST_Length(..., true) instead of only ST_Length(...) – JGH Nov 24 '20 at 20:58
  • Ah, you are right @Taras, it should be in the virtual layer, thank you very much. I successfully created a virtual layer with a "distance" column using your query. Just wanted to make sure, can I ask if my understanding is correct? The calculated distance is the distance from the centroid to its nearest line (other lines are not calculated), is this right? Thank you! – gislearner Nov 24 '20 at 22:30
  • Thank you @JGH! I think if a perpendicular line is impossible, I should use the nearest end point of the line. – gislearner Nov 24 '20 at 22:32
  • @Taras I'm so sorry, but can I also ask what's the unit of calculated distances? The numbers are different from those calculated by NNJoin, so I assume not in degrees? But it doesn't look like km/mile either.. – gislearner Nov 24 '20 at 22:50
3

The ST_ShortestLine is what is needed here, as it will be the "perpendicular" distance or to the nearest endpoint.

By default, it considers the polygon edge, not its centroid, so the query needs to compute the centroid first.

If working with simple polygons, you can remove the group by boundaries.id. If you need it, you would also have to group by line.id to keep the distinct polygons and lines and to apply an aggregate function on the distance computation (MIN(round...))

SELECT "boundaries".*,
       "lines".id,
       -- MIN( --if you use group by, you need the min() distance
       ROUND(
         ST_Length(
           ST_ShortestLine(
            ST_Centroid("boundaries".geometry), 
            "lines".geometry), 
         true),
       2)
     -- )  --if you use the group by, you need the min() function
     AS distance
FROM "boundaries", 
     "lines"
--GROUP BY boundaries.id, lines.id

Let's note that the input is in 4326, so in lat-longs. The doc says if you use st_length(geometry) the output is in CRS units (degrees here), but if you use st_length(geometry, true) or st_length(geometry,false), the output will be in meters

3
  • Thank you very much, @JGH! It worked perfectly. I was wondering if I could check the unit with you just to be sure? Thank you! After looking at the numbers, I'm assuming the resulted distance is in meters? The CRS of the layer is EPSG: 4326 - WGS 84. – gislearner Nov 30 '20 at 19:24
  • Yes, in meters. See edited answer – JGH Nov 30 '20 at 19:29
  • Thank you for replying, and for sharing the documentation! It is very helpful. – gislearner Nov 30 '20 at 19:35
1

I may (likely) have misinterpreted your question, reading perpendicular as strictly north-south or east-west facing ...the shortest distance perpendicular to the actual lines is trivially solved with ST_ShortestLine...I leave this here in case anyone ever wants it the other way.


Probably the easiest way to get all possible perpendicular distances (and the corresponding connecting LineStrings) is to work out y = mx+b and x = (y-b)/m of the line equations; this assumes, however, that your LineStrings are two-vertice segments!

The query (to use in a Virtual Layer) is a little convoluted without the option to add custom functions (and it would then be more versatile to implement proper vector dot product intersections instead):

SELECT poly_id, line_id, 'vertical' AS perp,
       ST_Distance(centroid, MakePoint(x, m*x+(ln_y-m*ln_x)), true) AS dist,
       MakeLine(centroid, MakePoint(x, m*x+(ln_y-m*ln_x))) AS geometry
FROM   (
    SELECT lns.<id> AS line_id, ply.<id> AS poly_id,
           (ST_Y(ST_EndPoint(lns.geometry)) - ST_Y(ST_StartPoint(lns.geometry))) / (ST_X(ST_EndPoint(lns.geometry)) - ST_X(ST_StartPoint(lns.geometry))) AS m,
           ST_X(ST_Centroid(ply.geometry) AS x,
           ST_X(ST_StartPoint(lns.geometry)) AS ln_X, ST_Y(ST_StartPoint(lns.geometry)) AS ln_y,
           ST_Centroid(ply.geometry) AS centroid
    FROM   <polys> AS ply
    JOIN   <lines> AS lns
      ON   ST_X(ST_Centroid(ply.geometry)) BETWEEN MIN(ST_X(ST_StartPoint(lns.geometry)), ST_X(ST_EndPoint(lns.geometry))) AND MAX(ST_X(ST_StartPoint(lns.geometry)), ST_X(ST_EndPoint(lns.geometry)))
) q
UNION ALL
SELECT poly_id, line_id, 'horizontal' AS perp,
       ST_Distance(centroid, MakePoint((y-(ln_y-m*ln_x))/m, y), true) AS dist,
       MakeLine(centroid, MakePoint((y-(ln_y-m*ln_x))/m, y)) AS geometry
FROM   (
    SELECT lns.<id> AS line_id, ply.<id> AS poly_id,
           (ST_Y(ST_EndPoint(lns.geometry)) - ST_Y(ST_StartPoint(lns.geometry))) / (ST_X(ST_EndPoint(lns.geometry)) - ST_X(ST_StartPoint(lns.geometry))) AS m,
           ST_X(ST_Centroid(ply.geometry) AS x,
           ST_X(ST_StartPoint(lns.geometry)) AS ln_X, ST_Y(ST_StartPoint(lns.geometry)) AS ln_y,
           ST_Centroid(ply.geometry) AS centroid
    FROM   <polys> AS ply
    JOIN   <lines> AS lns
      ON   ST_Y(ST_Centroid(ply.geometry)) BETWEEN MIN(ST_Y(ST_StartPoint(lns.geometry)), ST_Y(ST_EndPoint(lns.geometry))) AND MAX(ST_Y(ST_StartPoint(lns.geometry)), ST_Y(ST_EndPoint(lns.geometry)))
) q
;

This also assumes your geometries to be referenced in EPSG:4326, and the resulting dist will be in meter.

1
  • Thank you so much, @geozelot! This is really really helpful! Due to complexity, I might be using the distance to the nearest line for now, but I'm sure I'll be using this in the future. – gislearner Nov 30 '20 at 19:26

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.