11

I am quite new to Postgres and PostGIS.

I am trying to do a simple map-matching (or path-inference if you wish) of some GPS points. I am using postgres(postgis).

I have a table in my database containing my GPS points: gpspoints (Lat, Long, time, speed, azimuth, geometry, ...)

I have imported a shapefile of the road network into my database as well: routes (gid, idrte,version, nomrte, norte, clsrte, geom, ...) -- I DO NOT HAVE the azimuth of the links, and I am not sure how can I calculate it using the postgis function ST_azimuth. All I have is a shapefile of the road network containing columns enumerated above.

I want to associate each GPS point to the nearest link (in a buffer of 20 meters around the point) only if the direction of the link agrees with the azimuth of the GPS point (+ or -15 degrees) and retrieve the projected position. Otherwise I want it to search for the NEXT nearest link, within the 20 meter buffer, which has an acceptable azimuth! (Just like the picture!)

I want the new coordinates of the projected GPS points to be added in the table "gpspoints" as "projectedLat" and "projectedLong".

(In the picture below, actual points are demonstrated using a direction while projected points do not have any direction mark) enter image description here

  • It is an interesting problem. I assume that if the next point within the 20m buffer also has an unacceptable azimuth, then you want to continue the search? – John Powell Apr 26 '15 at 8:46
  • Yes John, but just for the road segments which are withing the 20 meter buffer! Actually, I want to identify the road segments which are in the 20 m buffer, then find the closest segment which has an acceptable azimuth! – Nader Apr 27 '15 at 13:48
  • "I DO NOT HAVE the azimuth of the links, and I am not sure how can I calculate it using the postgis function ST_azimuth. All I have is a shapefile of the road network containing columns enumerated above." You can (and should) bring your roads into Postgres/PostGIS with shp2pgsql. Then you can use ST_Azimuth, but bear in mind that it will tell you the azimuth with respect to the direction of digitisation. – alphabetasoup Oct 11 '15 at 3:21
1

I don't have a fully worked out answer, but maybe enough to get you started. These functions might help:

ST_Line_Locate_Point() gives the distance as a portion of the total line length of a point along a line. ST_Line_Interpolate_Point() returns a point geometry for a point at a given distance (again as a portion of the total length) along a line.

What I suggest it to get the point geometry of a location 'new_pts' on your line which are closest to your gpspoints (in a subquery below, formulated as a WITH clause). Then use that point to find the nearest point actually on the line, but at a very slightly smaller distance, so it will be closer to the beginning. Next use this found point, together with the 'new_pt' to get the azimuth between them. This would more or less represent the tangent to the line at that point.

Then you can compare this angle to the azimuth of your gpspoints to determine whether to add the gpspoint or not.

with new_pts AS (
    SELECT ST_ClosestPoint(geom, gpspoints) AS geompt
    FROM routes JOIN gpspoints ON ST_Distance(routes.geom, gpspoints.geom)<20
    )


SELECT ST_Azimuth(new_pts.geompt, 
    ST_Line_Interpolate_Point(routes.geom, 
        ST_Line_Locate_Point(routes.geom, new_pts.geompt)*0.999))

FROM routes JOIN new_pts ON ST_Distance(routes.geom,new_pts.geompt)<20;

Again, not tested, but I hope it will be some help.

0

try something like this:

  1. calculate the orientation of your polylines in degrees, example: e.g. in the ArcGIS Field Calculator: 180+math.atan2(( !Shape.firstpoint.X! - !Shape.lastpoint.X! ),( !Shape.firstpoint.Y! - !Shape.lastpoint.Y! ) ) * (180 / math.pi )

  2. convert your road network to road vertices (points) and import to PostGIS.

  3. apply matching algorithm:

    INSERT into matching(gid, vehicle_id, segment_id, date, time, the_geom)

    SELECT DISTINCT ON(b.gid) b.gid, b.vehicle_id, a.segment_id, b.date, b.time, ST_Closestpoint(ST_Collect(a.the_geom), b.the_geom) as the_geom FROM

    segment_vertices_geom a inner join vehicle_geom b ON ST_Dwithin(ST_Transform(a.the_geom,32632),

    ST_Transform(b.the_geom,32632), 20) WHERE ((CAST(a.azimuth AS float8) - CAST(b.direction AS float8)) < 180 OR (CAST(a.azimuth AS float8) - CAST(b.direction AS float8)) > -180) GROUP BY b.gid, b.vehicle_id, a.segment_id, b.date, b.time, a.the_geom, b.the_geom ORDER BY b.gid, ST_Distance(a.the_geom,b.the_geom);

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.