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

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

Answer 1

[Chad Greene]

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.

[pietro]

That is a clever solution. The GPS data I have are expressed in geodetic coordinates. How could I convert them to xyz? I have tried the deg2utm function (from filexchange), but roadx and roady data are far from the converted GPS data. Probably the ref point of the convertion carried out by 'usegeocoords' is from that of deg2utm function.

[Chad Greene]

Use projfwd.

[pietro]

but how can I get the scalar map projection structure from my shapefile?

[pietro]

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?