7

Setup: PostGIS 3.1, QGIS 3.16

I have a bunch of geographical data (geometries). They show more or less the same route, just measured several times. Due to the many measurements and the GPS precision, the measured points differ, of course:

General visualization of the data:

enter image description here

A more detailed view:

enter image description here

Next to the coordinates the data points keep some additional data, for example the speed of the vehicle which delivered the measurements. This is displayed by the color.

I am interested in the average speed at the track. So I did some interpolation (clustering?) and aggregated the points using PostGIS' ST_SnapToGrid() function:

SELECT
    ST_SnapToGrid(geom, 0.001) as c_geom,
    avg(speed) as avg_speed
FROM my_data
GROUP BY 1

This, naturally, yield a grid-like aggregation:

enter image description here

However, I am searching for a function which I can use like the ST_SnapToGrid() function, but which do not simply trunc the digits of the latitude and longitude but does more like an averaging, so the resulting points are more located on the track itself ("snap to track"). Of course, the speed value needs to be averaged as well.

Desired output:

I am searching for a function which averages lat/lon/speed values within a 10 meters radius or something like that:

enter image description here

I am not even sure if this is possible.

4
  • 2
    Snapping those points to railway tracks is far less complex than matching to a road network. Do you happen to have a railway network data set? – geozelot Feb 10 at 21:36
  • 1
    @geozelot The data is railway data, yes. But this is the only data I have. I don't have any reference data. I believe, I can take some OpenStreetMap data if it helps. This could be a backup strategy and I'd appreciate to see how it works. However, I would love to see a way without the reference data as well, since I am not quite sure whether I am allowed to use it or not. – S-Man Feb 11 at 7:12
  • I am currently getting to the point where I will have to sort out the same problem, keen to add bounty on this question, if I find some good solution to this later I will definitely share my findings. Do your point have speed and bearing information? – Miro Feb 14 at 0:03
  • @Miro Yes I have bearing information as well. :) – S-Man Feb 14 at 8:29
4
+50

The solution presented here consists basicaly of the following steps:

  1. Project the points to the line layer
  2. Create a new point layer with points at a regular interval along the line
  3. For each point of step 2, find the closest point from step 1: this one is the "cluster center": the point that remains and that gets the mean value of the other points in a certain distance (those that are deleted afterwards)
  4. For each cluster center (from step 3) create a buffer and get the mean value for a certain attribute (like speed) of all projected points (from step 1) that fall within the buffer.

So this is how you can do it using QGIS:

  1. Project the points to the line, using an expression like the following on the point layer - us Menu Processing / Toolbox / Geometry by expression to create te points as actual geometries - attributes will be kept:
closest_point (
    collect_geometries (
        overlay_nearest ( 
            'strecken_polyline', 
            $geometry
        )
    ), 
    $geometry
)

enter image description here

  1. Menu Processing / Toolbox / Points along geometry: create regular points along the line - select a distance that fits your date and needs. On the next screenshot, the blue dots are the points created here:

enter image description here

  1. For each blue dot, find the closest white dot: this will be used as cluster center (thus the points that will be kept and that get the mean value off all other points). This is the expression to use with Menu Processing / Toolbox / Geometry by expression, where projected_points is the name of the point layer created in step 1:
collect_geometries (
    overlay_nearest( 
        'projected_points', 
        $geometry, 
        limit:=1
    )
)

In the next screenshot, the cluster center is marked by red arows:

enter image description here

  1. Now create a buffer around each cluster center (result from step 3). For the size of the buffer, for demonstration purpose I selected a value of 50 [meters] - use whatever distance fits your data and change the value on line 7 of the following expression. Than, for each point that falls inside this buffer, calculate the mean value for an attribute (here: the field named speed). You can do this by opening the attribute table of the layer with the cluster center (result from step 3) and use the field calculator. Introduce this expression to create a new field with the mean value:
aggregate( 
    'projected_points', 
    'mean', 
    speed, 
    intersects (
        $geometry, 
        buffer (geometry (@parent), 50)
    )
)

See the result on the next screenshot: the original points (red dots) are reduced to the white squares. Compare the values of the labels that represent the speed field (I used random values from 80 to 220, just to see the effect): the small values nearby the red dots is the speed values in the original data, the bold values is the calculated mean for the points that fall inside the blue buffer (visualized here for better understanding, but it is not necessary to actually create these buffers):

enter image description here

Remarks: as you can see, depending on how far away the nearest projected point (step 1) is from the regular points (step 2), the buffers overlap or have gaps in between - so in some cases not all points are "catched" or some points are "catched" twice. I guess in your data, you have much more points and they are closer together, so this problem should not be too big.

However, there is a possibility to make sure that a) all points are taken into consideration and b) every point is taken into consideration just once. If you would like to do it that way, it's a bit more complicated, but not too much - just leave a comment so that I can add this to the solution.


This here is the first variant of the solution - I keep it here as it still might be helpful:

If you have a line-layer (railway, from OpenStreetMap), you can do the following steps. In principle, it works also without an additional line layer, just with the points (simply skip step 1), but than you only have a mean value for each point, based on all point in a certain distance (that you set in the step numbered with 2 below). To create a line from these points would be another question that should be asked separately.

  1. Create a buffer around each point with a distance that fits your data.

  2. Create a new field on the buffer layer that counts the mean of all values of a particular attribute field of those points that fall inside the buffer. Use this expression (adapted to to filed/layernames you use - I used projected_points for the layer from step 1, and value for an attribute in the original points layer with random values - your field is the one containing the speed information):

    aggregate( 
        'projected_points', 
        'mean', 
        value, 
        intersects (
            $geometry, 
            geometry (@parent)
        )
    )

See screenshot, where the expression is used for labeling the buffer: enter image description here

  1. If you want, using a table join, you can join the value calculated in the last step to your original point layer.
13
  • Great, I am giving it a try. Must fetch some appropriate railway data first :) – S-Man Feb 15 at 16:48
  • Hi, thanks for help so far. However, unfortunately I am struggling at the very first step. I have two layers: My points and the shapefile from data.deutschebahn.com/dataset/geo-strecke Now I calculated a new layer using Geometry Expression (instead of 'line' I used the name 'strecken_polyline'). After some calculating, a new layer is created, but it is empty. What I am doing wrong? Could you help please? :) – S-Man Feb 17 at 16:09
  • I updated the first step of the answer to get you an expression that should also work with the polyline you use. The expression I used was just for demonstration purpose and needs to be adapted to the situation you have. As you use QGIS 3.16 (as stated in your question), you will be able to use the function overlay_nearest that was introduced in this version. – Babel Feb 17 at 21:06
  • Thanks. I am using Toolbox -> Vector geometry -> Buffer to create the buffers. As distance, only geodetic degrees is available, which creates no circles. Is there a way to define "10 meters" instead? Furthermore, could you please describe how to "create a new field on the buffer layer". This is not clear to me, since I am completely new to QGIS... Thanks in advance! :) – S-Man Feb 18 at 11:23
  • Sorry for answering... I cannot see the point, were the clustering happens. Currently I have (for example) 1000 points around the line. These were snapped to 1000 points ON the line. Now I created 1000 buffers (incl. some AVG values). Now I am missing the most important step for me: How could I reduce the 1000 values to, maybe 100. My original problem was: Creating a point which contains the average of all surroundings to remove the surroundings... – S-Man Feb 18 at 12:10
5
+100

A solution using PostGIS alone.
The Steps are as follows:

  1. Generate clusters from the points
  2. Get the centre of each point cluster
  3. Generate a line that connects each cluster centre point (will also approximate the centreline of the route)
  4. Create sampling points at 20m intervals along the line (allowing for 10m radius buffer from each sampling point)
  5. Get an average value for 10m buffer around each sampling point.
--1. Create the point clusters
DROP TABLE IF EXISTS points_clustered;
CREATE TABLE points_clustered AS
SELECT geom, st_clusterkmeans(geom, 10) over () AS cluster
FROM points;

--2. Get the centre of each point cluster
DROP TABLE IF EXISTS centers;
CREATE TABLE centers AS
SELECT cluster, ST_Centroid(ST_collect(geom)) AS geom
FROM points_clustered 
GROUP BY cluster;

--3. Get centreline of connecting points
DROP TABLE IF EXISTS clustersline;
CREATE TABLE clustersline AS
SELECT st_makeline(geom) geom FROM 
    (SELECT geom AS geom FROM centers GROUP BY geom ORDER BY st_closestpoint(geom, geom)) AS f;

--4. Create a point on the line every 20 metres **(CHANGE THE EPSG)**
DROP TABLE IF EXISTS clustersline_20m;
CREATE TABLE clustersline_20m AS
WITH line AS 
    (SELECT
        (ST_Dump(geom)).geom AS geom
    FROM clustersline),
linemeasure AS
    (SELECT
        ST_AddMeasure(line.geom, 0, ST_Length(line.geom)) AS linem,
        generate_series(0, ST_Length(line.geom)::int, 20) AS i
    FROM line),
geometries AS (
    SELECT
        i,
        0.0 avg,
        (ST_Dump(ST_GeometryN(ST_LocateAlong(linem, i), 1))).geom AS geom 
    FROM linemeasure)
SELECT
    i,
    avg,
    ST_SetSRID(ST_MakePoint(ST_X(geom), ST_Y(geom)),28355) AS geom
FROM geometries;

--5. Get an average value for 10m buffer of each sampling point.
UPDATE clustersline_20m 
SET avg = subq.avg
FROM 
(
SELECT id, AVG(p.value) avg
FROM points p, clustersline_20m a
WHERE ST_INTERSECTS(ST_BUFFER(a.geom,10), p.geom)
GROUP BY id
) AS subq
WHERE id = subq.id;

Process Steps Image

Credit to these questions which helped build the answer:
http://blog.cleverelephant.ca/2018/06/polygon-splitting.html
PostGIS or QGIS: Convert unsorted points to a single line by connecting each 2 closest points
How can I transform polylines into points every n metres in PostGIS?

11
  • That looks really nice, great work so far. However at the first step: You give a fixed input for the kMeans cluster of 10. I am not really into this topic, but AFAIK this defines the number of clusters, correct? Well, the problem is that I don't know how many clusters I will need. I just know, that they should be 10m away from each other... Is there a more generic way for this parameter? – S-Man Feb 15 at 16:16
  • 1
    Yes the fixed value of 10 for kMeans is to define the number of clusters. In this case the clusters are being used to define the 'nodes' of a line from which we will create sampling points, not to define the sampling points themselves. The sampling points are created in step 4 from the line. If the kmeans value is increased, the more edges the line will have. Saying that, an alternate method would be to use something like ST_clusterdbscan where the inputs are the desired distance and density: SELECT geom, ST_clusterdbscan(geom, 20, 1) OVER() AS cluster – Cushen Feb 15 at 22:24
  • I've update the image to show the sampling points more clearly (in blue steps 4 and 5). – Cushen Feb 16 at 8:08
  • 1
    @bugmenot123: I guess the normal clusterdbscan would not work because all of my points in the data set are "reachable" from each other according to the definition of DBScan (en.wikipedia.org/wiki/DBSCAN). Indeed, when I run the function, it returns exact one cluster. This is because I have data for every single meter of the railway track. So, I need to limit the cluster size or something. The example given in the answer would also lead to exactly one cluster because all points are reachable from each other. DBSCAN would not deliver several clusters. Am I wrong? – S-Man Feb 16 at 13:15
  • 1
    Because both solutions are very interesting and great responses, I'd like to bounty both of it. So far I am not decided, which one should be accepted because both have advantages and disadvantages. However, thanks for your great input so far to both of you! – S-Man Feb 20 at 8:16
0

You actually seems to try to do something called map-matching. I don't know at wich scale do you want to use your algorithm, but if you just group GPS track by cluster or even snap it on the closest road you will face a lot of problems with real life data. GPS data can easily be noisy, especially in urban environnement, where there is a lot of chances to snap on the wrong road, because the network is more dense.

If you goal is to group data which have passed on the same road, either to count, compare vehicles with each other or gather info on the roads, I really encourage you to use map-matching first.

There is a lot of different ways to do this, for exemple:

  • Using a paying service: Mapbox, TomTom, Here...
    • Sometimes free for demo
    • Can be costly, and dependant of their map
  • Using a free software to install yourself and run with OSM data: GraphHopper, OSRM, ...
    • Free, sometimes there is a demo server too, work with the data you give it (OSM format generally), can run the whole world
    • Needs to be installed and managed, need a big server if you want to run the world efficiently
  • Using small scale map-matching tools

Once your data is properly snapped to your network, it's simpler to cut it like you want, either by using the geometry of your network or by using the snapped points.

Also, this way of trying to cut roads into regular part is also a complicated problem (intersections, map changing over time, ...) and is usually linked to the idea of trying to find concentration of something along the road (event, speed, ...). Your can find more info on this subject for exemple by looking here (where they introduce the notion of lixel - linear pixel): https://www.sciencedirect.com/science/article/pii/S0198971508000318

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.