How to check a GPS coordinate is located on a road?

13 views (last 30 days)
Hi all,
I have GPS data collected on a vehicle and I need to check the portion of time the vehicle was used on road. To this goal, I have downloaded the shape file of the road network of my area. My algorithm is quite simple: for each GPS coordinate, I compute a circle of 15m and I compute the intersection between the lines in the shape file and the circle. No intersection means the vehicle was not running on a road. Here you can see my code:
function OnRoad=CheckOnRoad(Lat,Lon,ShapeLat,ShapeLon)
rad = 15/1e3; % [m]edit
OnRoad=false(size(Lat));
for i=1:1:length(Lat)
[lonc,latc] = scircle1(Lon(i),Lat(i),km2deg(rad),[],[],[],40);
[lati, loni] = polyxpoly(latc,lonc,ShapeLat,ShapeLon);
OnRoad(i)=not(isempty(loni));
end
However, this code is pretty slow. How could I speed up the code?
Thanks. Best regards,
Pietro

Accepted Answer

Chad Greene
Chad Greene on 17 Oct 2018
Here's a quick way to do it. Below I'm using standard sample data that comes with the Mapping Toolbox, so you should be able to follow along.
The first thing I'd do is convert all your lat,lon data to some projected x,y data. You can either do that by specifying 'usegeocoords',false when you call shaperead, or you can use projfwd. Then you'll be working in meters, and that should be much faster than calling scircle1 every time you go through the loop.
Run this:
% Load some road data:
S = shaperead('concord_roads.shp','usegeocoords',false);
roadx = [S(:).X];
roady = [S(:).Y];
% Generate 1000 random points in the area:
gpsx = mean(roadx,'omitnan') + 1000*randn(1000,1);
gpsy = mean(roady,'omitnan') + 500*randn(1000,1);
plot(roadx,roady,'k')
hold on
h = plot(gpsx,gpsy,'r.');
axis image
Those are a bunch of roads and a thousand random points in the domain equivalent to your gps points. Now go through every GPS point and calculate the distance from that point to all the road coordinates:
% Preallocate:
d = NaN(size(gpsx));
% Loop through each GPS point:
for k = 1:length(gpsx)
% Calculate the distance from this gps point to all road
% points and log just the minimum distance:
d(k) = min(hypot(roadx-gpsx(k),roady-gpsy(k)));
end
Now we have the distance from each GPS point to the closest road. We can plot those distances as a scatter plot like this:
delete(h) % removes previous red dots
scatter(gpsx,gpsy,10,d,'filled')
cb = colorbar;
ylabel(cb,'distance to nearest road (m)')
caxis([0 200])
Now you can find all the gps points within a given distance from the nearest road by one simple logical expression. Let's mark all the GPS points within 15 m of a road using red x markers:
OnRoad = d<15;
plot(gpsx(OnRoad),gpsy(OnRoad),'rx')
This approach should be faster and more straightforward than the scircle1 approach. If speed is still an issue, it can be sped up further by removing the hypot call in the loop. That's because hypot takes a square root, and taking square roots can be slow. Here's a way you could just take the square root once instead of using hypot every time through the loop:
% Preallocate:
d_sq = NaN(size(gpsx));
% Loop through each GPS point:
for k = 1:length(gpsx)
% Calculate the distance from this gps point to all road
% points and log just the minimum distance:
d_sq(k) = min( (roadx-gpsx(k)).^2 + (roady-gpsy(k)).^2 );
end
% Just take the square root once:
d2 = sqrt(d_sq);
Hope that helps.
  4 Comments
pietro
pietro on 20 Oct 2018
but how can I get the scalar map projection structure from my shapefile?
pietro
pietro on 10 Feb 2019
Edited: pietro on 11 Feb 2019
Hi,
I have solved the problem, I have converted the geodetic coordinateds to cartersian using the geodetic2ecef function (as (*) . However, using the scircle1 method (where the circle radius was only 15m ), I should set a threshold distance of more than 250m. That is a bit confusing. In your opinion, why this large difference in the two methods?
Thanks
(*)
wgs84 = wgs84Ellipsoid('meters');
[gpsx,gpsy,gpsz] = geodetic2ecef(wgs84,LatN,LonN,AltN);
edit, I have found the problem. The road shape file has only the point of road change (e.g. a curve, a cross and so on). So for a long straigth road, there are only the two road ends and when the vehicle is located between the two road ends, the distance between the road and vehicle results very large. Instead, when a vehicle is leaving the road, it is located close to a road point so the distance results very small.
So, point-by-point is every efficient but it cannot work, at least fit the road data or calculate the distance between a point (location of the vehicle) and the line (the road stretch). But, how can I implement it?

Sign in to comment.

More Answers (0)

Categories

Find more on Vehicle Scenarios in Help Center and File Exchange

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!