Double precision for lon lat (valhalla#2693)

* Introduce Point2d and Vector2d

* Use double precision for GeoPoint

* Read ll from json as double

* Fix unit tests

* Fix double precision for height action

* Continue fix tests

* hopefully reduce pointll hash collisions, fix long standing but benign bugs in decode7/encode7

* allow differnt precision on encode/decode7

* use 7 digits of precision for finding other places precision is lost

* make 7digit precision possible with a preprocessor defintion

* store 7bits of precision in OSMNode daat tile building time

* Fix hash collisions for PointXY

* Fix tiles unit test

* Don't loose precision while setting lon,lat to OSMNode

* Fix test_trimming

* Fix thor_worker test

* Fix matrix test

* 7 digits of precision for nodeinfo ll

* refactor TileHierarchy while tracking down a typo :(

* update isochrone test

* update summary shape

* fix map match test

* fix astar test

* lint

* update summary

* use double precision for path edges in loki search

* use our own hash_combine since older boost doesnt ahve it at same place

* more doubles in place of floats

* more doubles to shut up clang

* working on bss tests

* fixed one test in the suite

* slight differences to bss testing

* roll our own hash_combine for since boost relocates it

* forgot to remove header

* explicit instantiations

* tidy

* smarter isochrone test

* cleanup comments etc

* changelog derp

Co-authored-by: Kevin Kreiser <[email protected]>

Author: genadz

CHANGELOG.md

@@ -199,6 +199,7 @@
    * ADDED: Added functions for predicted speeds encoding-decoding [#2674](https://github.com/valhalla/valhalla/pull/2674)
    * ADDED: Time invariant routing via the bidirectional algorithm. This has the effect that when time dependent routes (arrive_by and depart_at) fall back to bidirectional due to length restrictions they will actually use the correct time of day for one of the search directions [#2660](https://github.com/valhalla/valhalla/pull/2660)
    * ADDED: If the length of the edge is greater than kMaxEdgeLength, then consider this a catastrophic error if the should_error bool is true in the set_length function. [2678](https://github.com/valhalla/valhalla/pull/2678)
+   * ADDED: Moved lat,lon coordinates structures from single to double precision. Improves geometry accuracy noticibly at zooms above 17 as well as coordinate snapping and any other geometric operations. Addes about a 2% performance pentalty for standard routes. Graph nodes now have 7 digits of precision.  [#2693](https://github.com/valhalla/valhalla/pull/2693)
 
 ## Release Date: 2019-11-21 Valhalla 3.0.9
 * **Bug Fix**

proto/tripcommon.proto

@@ -3,8 +3,8 @@ option optimize_for = LITE_RUNTIME;
 package valhalla;
 
 message LatLng {
-  optional float lat = 1;
-  optional float lng = 2;
+  optional double lat = 1;
+  optional double lng = 2;
 }
 
 message BoundingBox {
@@ -34,10 +34,10 @@ message Location {
 
   message PathEdge {
     optional uint64 graph_id = 1;
-    optional float percent_along = 2;
+    optional double percent_along = 2;
     optional LatLng ll = 3;
     optional SideOfStreet side_of_street = 4;
-    optional float distance = 5;
+    optional double distance = 5;
     optional int32 minimum_reachability = 6; //deprecated
     optional bool begin_node = 7;
     optional bool end_node = 8;

src/baldr/connectivity_map.cc

@@ -116,7 +116,7 @@ connectivity_map_t::connectivity_map_t(const boost::property_tree::ptree& pt) {
   // See what kind of tiles we are dealing with here by getting a graphreader
   GraphReader reader(pt);
   auto tiles = reader.GetTileSet();
-  transit_level = TileHierarchy::levels().rbegin()->second.level + 1;
+  transit_level = TileHierarchy::GetTransitLevel().level;
 
   // Quick hack to remove connectivity between known unconnected regions
   // The only land connection from north to south america is through
@@ -142,12 +142,12 @@ connectivity_map_t::connectivity_map_t(const boost::property_tree::ptree& pt) {
   // (build the ColorMap). Transit level uses local hierarchy tiles
   for (auto& color : colors) {
     if (color.first == transit_level) {
-      TileHierarchy::levels().rbegin()->second.tiles.ColorMap(color.second, not_neighbors);
+      TileHierarchy::GetTransitLevel().tiles.ColorMap(color.second, not_neighbors);
     } else {
-      TileHierarchy::levels()
-          .find(color.first)
-          ->second.tiles.ColorMap(color.second,
-                                  color.first == 2 ? not_neighbors : decltype(not_neighbors){});
+      TileHierarchy::levels()[color.first].tiles.ColorMap(color.second,
+                                                          color.first == 2
+                                                              ? not_neighbors
+                                                              : decltype(not_neighbors){});
     }
   }
 }
@@ -173,7 +173,7 @@ std::unordered_set<size_t> connectivity_map_t::get_colors(uint32_t hierarchy_lev
   if (level == colors.cend()) {
     return result;
   }
-  const auto& tiles = TileHierarchy::levels().find(hierarchy_level)->second.tiles;
+  const auto& tiles = TileHierarchy::levels()[hierarchy_level].tiles;
   std::vector<const decltype(location.edges)*> edge_sets{&location.edges, &location.filtered_edges};
   for (const auto* edges : edge_sets) {
     for (const auto& edge : *edges) {
@@ -199,10 +199,10 @@ std::unordered_set<size_t> connectivity_map_t::get_colors(uint32_t hierarchy_lev
 std::string connectivity_map_t::to_geojson(const uint32_t hierarchy_level) const {
   // bail if we dont have the level
   uint32_t tile_level = (hierarchy_level == transit_level) ? transit_level - 1 : hierarchy_level;
-  auto bbox = TileHierarchy::levels().find(tile_level);
-  if (bbox == TileHierarchy::levels().cend()) {
+  if (tile_level >= TileHierarchy::levels().size()) {
     throw std::runtime_error("hierarchy level not found");
   }
+  const auto& tiles = TileHierarchy::levels()[tile_level].tiles;
 
   // make a region map (inverse mapping of color to lists of tiles)
   // could cache this but shouldnt need to call it much
@@ -230,7 +230,7 @@ std::string connectivity_map_t::to_geojson(const uint32_t hierarchy_level) const
   std::unordered_map<size_t, polygon_t> boundaries;
   for (const auto& arity : arities) {
     auto& region = *regions.find(arity.second);
-    boundaries.emplace(arity.second, to_boundary(region.second, bbox->second.tiles));
+    boundaries.emplace(arity.second, to_boundary(region.second, tiles));
   }
 
   // turn it into geojson
@@ -239,12 +239,12 @@ std::string connectivity_map_t::to_geojson(const uint32_t hierarchy_level) const
 
 std::vector<size_t> connectivity_map_t::to_image(const uint32_t hierarchy_level) const {
   uint32_t tile_level = (hierarchy_level == transit_level) ? transit_level - 1 : hierarchy_level;
-  auto bbox = TileHierarchy::levels().find(tile_level);
-  if (bbox == TileHierarchy::levels().cend()) {
+  if (tile_level >= TileHierarchy::levels().size()) {
     throw std::runtime_error("hierarchy level not found");
   }
+  const auto& level_tiles = TileHierarchy::levels()[tile_level];
 
-  std::vector<size_t> tiles(bbox->second.tiles.nrows() * bbox->second.tiles.ncolumns(),
+  std::vector<size_t> tiles(level_tiles.tiles.nrows() * level_tiles.tiles.ncolumns(),
                             static_cast<uint32_t>(0));
   auto level = colors.find(hierarchy_level);
   if (level != colors.cend()) {

src/baldr/graphreader.cc

@@ -606,7 +606,7 @@ GraphId GraphReader::GetShortcut(const GraphId& id) {
   };
 
   // No shortcuts on the local level or transit level.
-  if (id.level() >= TileHierarchy::levels().rbegin()->second.level) {
+  if (id.level() >= TileHierarchy::levels().back().level) {
     return {};
   }
 
@@ -799,7 +799,7 @@ std::unordered_set<GraphId> GraphReader::GetTileSet() const {
   } // or individually on disk
   else if (!tile_dir_.empty()) {
     // for each level
-    for (uint8_t level = 0; level <= TileHierarchy::levels().rbegin()->first + 1; ++level) {
+    for (uint8_t level = 0; level <= TileHierarchy::GetTransitLevel().level; ++level) {
       // crack open this level of tiles directory
       filesystem::path root_dir(tile_dir_ + filesystem::path::preferred_separator +
                                 std::to_string(level) + filesystem::path::preferred_separator);

src/baldr/graphtile.cc

@@ -410,8 +410,7 @@ GraphTile::FileSuffix(const GraphId& graphid, const std::string& fname_suffix, b
   */
 
   // figure the largest id for this level
-  auto found = TileHierarchy::levels().find(graphid.level());
-  if (found == TileHierarchy::levels().cend() &&
+  if (graphid.level() >= TileHierarchy::levels().size() &&
       graphid.level() != TileHierarchy::GetTransitLevel().level) {
     throw std::runtime_error("Could not compute FileSuffix for non-existent level: " +
                              std::to_string(graphid.level()));
@@ -420,7 +419,7 @@ GraphTile::FileSuffix(const GraphId& graphid, const std::string& fname_suffix, b
   // get the level info
   const auto& level = graphid.level() == TileHierarchy::GetTransitLevel().level
                           ? TileHierarchy::GetTransitLevel()
-                          : found->second;
+                          : TileHierarchy::levels()[graphid.level()];
 
   // figure out how many digits
   auto max_id = level.tiles.ncolumns() * level.tiles.nrows() - 1;
@@ -513,8 +512,7 @@ GraphId GraphTile::GetTileId(const std::string& fname) {
   }
 
   // if the first thing isnt a valid level bail
-  auto found = TileHierarchy::levels().find(digits.back());
-  if (found == TileHierarchy::levels().cend() &&
+  if (digits.back() >= TileHierarchy::levels().size() &&
       digits.back() != TileHierarchy::GetTransitLevel().level) {
     throw std::runtime_error("Invalid tile path: " + fname);
   }
@@ -524,7 +522,7 @@ GraphId GraphTile::GetTileId(const std::string& fname) {
   digits.pop_back();
   const auto& tile_level = level == TileHierarchy::GetTransitLevel().level
                                ? TileHierarchy::GetTransitLevel()
-                               : found->second;
+                               : TileHierarchy::levels()[level];
 
   // get the number of sub directories that we should have
   auto max_id = tile_level.tiles.ncolumns() * tile_level.tiles.nrows() - 1;
@@ -558,15 +556,9 @@ GraphId GraphTile::GetTileId(const std::string& fname) {
 
 // Get the bounding box of this graph tile.
 AABB2<PointLL> GraphTile::BoundingBox() const {
-
-  // figure the largest id for this level
-  auto level = TileHierarchy::levels().find(header_->graphid().level());
-  if (level == TileHierarchy::levels().end() &&
-      header_->graphid().level() == ((TileHierarchy::levels().rbegin())->second.level + 1)) {
-    level = TileHierarchy::levels().begin();
-  }
-
-  auto tiles = level->second.tiles;
+  const auto& tiles = header_->graphid().level() == TileHierarchy::GetTransitLevel().level
+                          ? TileHierarchy::GetTransitLevel().tiles
+                          : TileHierarchy::levels()[header_->graphid().level()].tiles;
   return tiles.TileBounds(header_->graphid().tileid());
 }
 

src/baldr/merge.cc

@@ -12,15 +12,6 @@ namespace baldr {
 namespace merge {
 
 namespace {
-
-uint64_t count_tiles_in_levels(GraphReader& reader) {
-  uint64_t tile_count = 0;
-  for (const auto& level : TileHierarchy::levels() | bra::map_values) {
-    tile_count += level.tiles.ncolumns() * level.tiles.nrows();
-  }
-  return tile_count;
-}
-
 namespace iter {
 
 // the "edges" struct provides a container wrapper for the edges leaving a node

src/baldr/nodeinfo.cc

@@ -72,7 +72,7 @@ NodeInfo::NodeInfo() {
 }
 
 NodeInfo::NodeInfo(const PointLL& tile_corner,
-                   const std::pair<float, float>& ll,
+                   const PointLL& ll,
                    const uint32_t access,
                    const NodeType type,
                    const bool traffic_signal) {
@@ -84,14 +84,21 @@ NodeInfo::NodeInfo(const PointLL& tile_corner,
 }
 
 // Sets the latitude and longitude.
-void NodeInfo::set_latlng(const PointLL& tile_corner, const std::pair<float, float>& ll) {
+void NodeInfo::set_latlng(const PointLL& tile_corner, const PointLL& ll) {
   // Protect against a node being slightly outside the tile (due to float roundoff)
-  lat_offset_ = ll.second < tile_corner.lat()
-                    ? 0
-                    : static_cast<uint32_t>((ll.second - tile_corner.lat()) / kDegreesPrecision);
-  lon_offset_ = ll.first < tile_corner.lng()
-                    ? 0
-                    : static_cast<uint32_t>((ll.first - tile_corner.lng()) / kDegreesPrecision);
+  lat_offset_ = 0;
+  if (ll.lat() > tile_corner.lat()) {
+    auto lat = std::round((ll.lat() - tile_corner.lat()) / 1e-7);
+    lat_offset_ = lat / 10;
+    lat_offset7_ = lat - lat_offset_ * 10;
+  }
+
+  lon_offset_ = 0;
+  if (ll.lng() > tile_corner.lng()) {
+    auto lon = std::round((ll.lng() - tile_corner.lng()) / 1e-7);
+    lon_offset_ = lon / 10;
+    lon_offset7_ = lon - lon_offset_ * 10;
+  }
 }
 
 // Set the index in the node's tile of its first outbound edge.

src/baldr/pathlocation.cc

@@ -5,9 +5,9 @@ namespace valhalla {
 namespace baldr {
 
 PathLocation::PathEdge::PathEdge(const GraphId& id,
-                                 const float percent_along,
+                                 const double percent_along,
                                  const midgard::PointLL& projected,
-                                 const float score,
+                                 const double score,
                                  const SideOfStreet sos,
                                  const unsigned int outbound_reach,
                                  const unsigned int inbound_reach)
@@ -47,8 +47,8 @@ bool PathLocation::shares_edges(const PathLocation& other) const {
     bool found = false;
     for (const auto& other_edge : other.edges) {
       if (edge.id == other_edge.id && edge.sos == other_edge.sos &&
-          midgard::equal<float>(edge.percent_along, other_edge.percent_along) &&
-          midgard::equal<float>(edge.distance, other_edge.distance, .1f) &&
+          midgard::equal<double>(edge.percent_along, other_edge.percent_along) &&
+          midgard::equal<double>(edge.distance, other_edge.distance, .1) &&
           edge.projected.ApproximatelyEqual(other_edge.projected)) {
         found = true;
         break;

src/baldr/tilehierarchy.cc

@@ -6,21 +6,25 @@ using namespace valhalla::midgard;
 namespace valhalla {
 namespace baldr {
 
-const std::map<uint8_t, TileLevel>& TileHierarchy::levels() {
+const std::vector<TileLevel>& TileHierarchy::levels() {
   // Static tile levels
-  static const std::map<uint8_t, TileLevel> levels_ =
-      {{2, TileLevel{2, stringToRoadClass("ServiceOther"), "local",
-                     midgard::Tiles<midgard::PointLL>{{{-180, -90}, {180, 90}},
-                                                      .25,
-                                                      static_cast<unsigned short>(kBinsDim)}}},
-       {1, TileLevel{1, stringToRoadClass("Tertiary"), "arterial",
-                     midgard::Tiles<midgard::PointLL>{{{-180, -90}, {180, 90}},
-                                                      1,
-                                                      static_cast<unsigned short>(kBinsDim)}}},
-       {0, TileLevel{0, stringToRoadClass("Primary"), "highway",
-                     midgard::Tiles<midgard::PointLL>{{{-180, -90}, {180, 90}},
-                                                      4,
-                                                      static_cast<unsigned short>(kBinsDim)}}}};
+  static const std::vector<TileLevel> levels_ = {
+
+      TileLevel{0, stringToRoadClass("Primary"), "highway",
+                midgard::Tiles<midgard::PointLL>{{{-180, -90}, {180, 90}},
+                                                 4,
+                                                 static_cast<unsigned short>(kBinsDim)}},
+
+      TileLevel{1, stringToRoadClass("Tertiary"), "arterial",
+                midgard::Tiles<midgard::PointLL>{{{-180, -90}, {180, 90}},
+                                                 1,
+                                                 static_cast<unsigned short>(kBinsDim)}},
+
+      TileLevel{2, stringToRoadClass("ServiceOther"), "local",
+                midgard::Tiles<midgard::PointLL>{{{-180, -90}, {180, 90}},
+                                                 .25,
+                                                 static_cast<unsigned short>(kBinsDim)}},
+  };
 
   return levels_;
 }
@@ -46,9 +50,8 @@ midgard::AABB2<midgard::PointLL> TileHierarchy::GetGraphIdBoundingBox(const Grap
 // If the level is not supported an invalid id will be returned.
 GraphId TileHierarchy::GetGraphId(const midgard::PointLL& pointll, const uint8_t level) {
   GraphId id;
-  const auto& tl = levels().find(level);
-  if (tl != levels().end()) {
-    auto tile_id = tl->second.tiles.TileId(pointll);
+  if (level < levels().size()) {
+    auto tile_id = levels()[level].tiles.TileId(pointll);
     if (tile_id >= 0) {
       id = {static_cast<uint32_t>(tile_id), level, 0};
     }
@@ -58,9 +61,9 @@ GraphId TileHierarchy::GetGraphId(const midgard::PointLL& pointll, const uint8_t
 
 // Gets the hierarchy level given the road class.
 uint8_t TileHierarchy::get_level(const RoadClass roadclass) {
-  if (roadclass <= levels().find(0)->second.importance) {
+  if (roadclass <= levels()[0].importance) {
     return 0;
-  } else if (roadclass <= levels().find(1)->second.importance) {
+  } else if (roadclass <= levels()[1].importance) {
     return 1;
   } else {
     return 2;
@@ -77,9 +80,8 @@ uint8_t TileHierarchy::get_max_level() {
 std::vector<GraphId> TileHierarchy::GetGraphIds(const midgard::AABB2<midgard::PointLL>& bbox,
                                                 const uint8_t level) {
   std::vector<GraphId> ids;
-  auto itr = levels().find(level);
-  if (itr != levels().end()) {
-    auto tile_ids = itr->second.tiles.TileList(bbox);
+  if (level < levels().size()) {
+    auto tile_ids = levels()[level].tiles.TileList(bbox);
     ids.reserve(tile_ids.size());
 
     for (auto tile_id : tile_ids) {
@@ -94,7 +96,7 @@ std::vector<GraphId> TileHierarchy::GetGraphIds(const midgard::AABB2<midgard::Po
 std::vector<GraphId> TileHierarchy::GetGraphIds(const midgard::AABB2<midgard::PointLL>& bbox) {
   std::vector<GraphId> ids;
   for (const auto& entry : levels()) {
-    auto level_ids = GetGraphIds(bbox, entry.first);
+    auto level_ids = GetGraphIds(bbox, entry.level);
     ids.reserve(ids.size() + level_ids.size());
     ids.insert(ids.end(), level_ids.begin(), level_ids.end());
   }
@@ -107,9 +109,8 @@ std::vector<GraphId> TileHierarchy::GetGraphIds(const midgard::AABB2<midgard::Po
  * @return Returns a const reference to the tiling system for this level.
  */
 const midgard::Tiles<midgard::PointLL>& TileHierarchy::get_tiling(const uint8_t level) {
-  const auto& tl = levels().find(level);
-  if (tl != levels().end()) {
-    return tl->second.tiles;
+  if (level < levels().size()) {
+    return levels()[level].tiles;
   }
   throw std::runtime_error("Invalid level Id for TileHierarchy::get_tiling");
 }

src/loki/height_action.cc

@@ -87,7 +87,7 @@ std::string loki_worker_t::height(Api& request) {
   }
 
   // get the distances between the postings if desired
-  std::vector<float> ranges;
+  std::vector<double> ranges;
   if (request.options().range()) {
     ranges.reserve(shape.size());
     ranges.emplace_back(0);

src/loki/matrix_action.cc

@@ -129,8 +129,7 @@ void loki_worker_t::matrix(Api& request) {
       if (!connectivity_map) {
         continue;
       }
-      auto colors =
-          connectivity_map->get_colors(TileHierarchy::levels().rbegin()->first, projection, 0);
+      auto colors = connectivity_map->get_colors(TileHierarchy::levels().back().level, projection, 0);
       for (auto& color : colors) {
         auto itr = color_counts.find(color);
         if (itr == color_counts.cend()) {

src/loki/node_search.cc

@@ -368,8 +368,8 @@ std::vector<baldr::GraphId> nodes_in_bbox(const vm::AABB2<vm::PointLL>& bbox,
                                           baldr::GraphReader& reader) {
   std::vector<vb::GraphId> nodes;
 
-  auto tiles = vb::TileHierarchy::levels().rbegin()->second.tiles;
-  const uint8_t bin_level = vb::TileHierarchy::levels().rbegin()->second.level;
+  const auto& tiles = vb::TileHierarchy::levels().back().tiles;
+  const uint8_t bin_level = vb::TileHierarchy::levels().back().level;
 
   // if the bbox only touches the edge of the tile or bin, then we need to
   // include neighbouring bins as well, in case both the edge and its opposite
@@ -425,8 +425,8 @@ std::vector<baldr::GraphId> nodes_in_bbox(const vm::AABB2<vm::PointLL>& bbox,
 std::vector<baldr::GraphId> edges_in_bbox(const vm::AABB2<vm::PointLL>& bbox,
                                           baldr::GraphReader& reader) {
 
-  auto tiles = vb::TileHierarchy::levels().rbegin()->second.tiles;
-  const uint8_t bin_level = vb::TileHierarchy::levels().rbegin()->second.level;
+  auto tiles = vb::TileHierarchy::levels().back().tiles;
+  const uint8_t bin_level = vb::TileHierarchy::levels().back().level;
 
   // if the bbox only touches the edge of the tile or bin, then we need to
   // include neighbouring bins as well, in case both the edge and its opposite