I have an OpenStreetMapX (mad) of Madrid and around 4 million coordinates that I wish to associate with their corresponding mad.nodes. Unfortunately using OpenStreetMapX.jl’s point_to_nodes out-of-the-box is infeasible (see here for benchmarks Good Julia function / algorithm idea to find the nearest node to a given point location based on some given road network (for example, open street map data)? - #5 by ctrebbau ). So I have to contruct an RTree with mad.nodes and then query that tree. I’ve done so with LibSpatialIndex like so
const LSI = LibSpatialIndex
madNodes = mad.nodes # This is a Dict{Int64, ENU}
rtree = LSI.RTree(3)
for (node, coord) in madNodes
LSI.insert!(rtree, node, [coord.east, coord.north, coord.up], [coord.east, coord.north, coord.up])
end
And then I query the tree with
queryPts = [LSI.knn(rtree, [pt.east, pt.north, pt.up], 1) for pt in enu_pts] # The 1 refers to the number of neighbors I want
where
enu_pts = [ENU(pt, mad.bounds) for pt in lla_pts]
and
lla_pts = [LLA(pt[1], pt[2]) for pt in pts]
while
pts = zip(df.dLat[1:100], df.dlng[1:100])
The problem is that sum(validPtNodes .== queryPts), where validPtNodes = [point_to_nodes(pt, mad) for pt in lla_pts] returns 45 out of the 100 queried. OpenStreetMapX implements euclidean distance and I have not been able to find in the documentation of LibSpatialIndex whether it does so as well, but I think it’s most probably the default(?).
So my question is, am I building the tree right? I’m I right to use point_to_nodes as a sanity check? Do I lose some precision when converting between LLA and ENU? Do you know of some resource that I can profit from?
P.S. here are the 100 pts I’m querying
pts = zip([40.40035, 40.47780807, 40.4377299, 40.4211153, 40.41384654, 40.471597, 40.48230943, 40.4239884, 40.406806, 40.44236981, 40.4262522, 40.4468401, 40.447563, 40.4396692, 40.43474505, 40.42054213, 40.491255, 40.4696799, 40.46801, 40.4296166, 40.42423842, 40.491255, 40.471597, 40.4358823, 40.4354445, 40.4585917, 40.406806, 40.4730513, 40.43313179, 40.4227419, 40.4256262, 40.491255, 40.46801, 40.4806616, 40.406806, 40.467844, 40.4637403, 40.4431511, 40.41807399, 40.4867882, 40.447689, 40.48837, 40.4295441, 40.45766912, 40.46682955, 40.41458, 40.446625, 40.4466221, 40.42212769, 40.422233, 40.432579, 40.438277, 40.4465975, 40.4277378, 40.39943841, 40.420281, 40.491255, 40.4167081, 40.46575663, 40.42791748, 40.4427623, 40.46801, 40.4349325, 40.471754, 40.50185837, 40.4348437, 40.448043, 40.43067, 40.4339336, 40.44618, 40.4236729, 40.4445497, 40.4634108, 40.4517406, 40.4125946, 40.465511, 40.385241, 40.422343, 40.472916, 40.41819319, 40.4235863, 40.375187, 40.446773, 40.46801, 40.44898857, 40.406806, 40.47537293, 40.4063306, 40.500982, 40.4481, 40.4209713, 40.4489254, 40.4556463, 40.44954419, 40.428623, 40.44696788, 40.42867831, 40.406806, 40.45108855, 40.471597], [-3.70306, -3.69849334, -3.6800531, -3.6879227, -3.69493838, -3.682403, -3.67061352, -3.6864495, -3.691947, -3.65546398, -3.7068227, -3.69488566, -3.656106, -3.8073816, -3.66857059, -3.57832935, -3.594576, -3.6698, -3.57002, -3.6876464, -3.68000284, -3.594576, -3.682403, -3.6421008, -3.7334977, -3.6945422, -3.691947, -3.6868434, -3.73357988, -3.6964978, -3.6844591, -3.594576, -3.57002, -3.6862549, -3.691947, -3.626451, -3.6139041, -3.6823055, -3.68979726, -3.6658943, -3.657577, -3.6959, -3.6835559, -3.66049308, -3.65845665, -3.69028, -3.695867, -3.6924595, -3.69565118, -3.6923577, -3.7076952, -3.637089, -3.6975215, -3.7295722, -3.68889302, -3.690747, -3.594576, -3.6910361, -3.61401077, -3.71710682, -3.7871334, -3.57002, -3.6799566, -3.849293, -3.66879389, -3.6968111, -3.624555, -3.6822437, -3.6924939, -3.71253, -3.7072785, -3.6543204, -3.6967773, -3.69269371, -3.6937105, -3.616573, -3.717764, -3.686257, -3.6636556, -3.68981805, -3.6684468, -3.6470258, -3.6936392, -3.57002, -3.6904071, -3.691947, -3.68825234, -3.6900896, -3.867037, -3.78518, -3.7080051, -3.6708406, -3.6800522, -3.69145327, -3.6951358, -3.69304122, -3.69419441, -3.691947, -3.69518925, -3.682403])
And you can download the map of Madrid here: Geofabrik Download Server
And set up like this:
using OpenStreetMapX
using LibSpatialIndex
const LSI = LibSpatialIndex
mad = get_map_data("madrid.osm", use_cache = false, trim_to_connected_graph=true);
Maybe @pszufe can give me further advice?
Please let me know if you need something else to have a minimum viable example and/or further conceptual clarification.
Thank you so much in advance!