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.
