Fix precision problems with very high-resolution tiles

code/prominence.cpp

@@ -198,9 +198,10 @@ int main(int argc, char **argv) {
   auto threadPool = std::make_unique<ThreadPool>(numThreads);
   int num_tiles_processed = 0;
   vector<std::future<bool>> results;
-  float lat = bounds[0];
+  // Use double precision to avoid accumulating floating-point error during loop
+  double lat = bounds[0];
   while (lat < bounds[1]) {
-    float lng = bounds[2];
+    double lng = bounds[2];
     while (lng < bounds[3]) {
       // Allow specifying longitude ranges that span the antimeridian (lng > 180)
       auto wrappedLng = lng;
@@ -209,17 +210,21 @@ int main(int argc, char **argv) {
       }
 
       std::shared_ptr<CoordinateSystem> coordinateSystem(
-        fileFormat.coordinateSystemForOrigin(lat, wrappedLng, utmZone));
+        fileFormat.coordinateSystemForOrigin(
+          static_cast<float>(lat), static_cast<float>(wrappedLng), utmZone));
 
       // Skip tiles that don't intersect filtering polygon
-      if (!filter.intersects(lat, lat + fileFormat.degreesAcross(),
-                             lng, lng + fileFormat.degreesAcross())) {
+      if (!filter.intersects(static_cast<float>(lat),
+                             static_cast<float>(lat + fileFormat.degreesAcross()),
+                             static_cast<float>(lng),
+                             static_cast<float>(lng + fileFormat.degreesAcross()))) {
         VLOG(3) << "Skipping tile that doesn't intersect polygon " << lat << " " << lng;
       } else {
         // Actually calculate prominence
         ProminenceTask *task = new ProminenceTask(cache.get(), options);
         results.push_back(threadPool->enqueue([=] {
-              return task->run(lat, wrappedLng, *coordinateSystem);
+              return task->run(static_cast<float>(lat), static_cast<float>(wrappedLng),
+                               *coordinateSystem);
             }));
       }
 

code/prominence_task.cpp

@@ -27,6 +27,7 @@
 #include "divide_tree.h"
 #include "island_tree.h"
 #include "tree_builder.h"
+#include "util.h"
 
 #include "easylogging++.h"
 
@@ -55,7 +56,7 @@ bool ProminenceTask::run(float lat, float lng, const CoordinateSystem &coordinat
     VLOG(2) << "Couldn't load tile for " << lat << " " << lng;
     return false;
   }
-
+  
   // Flip tile upside down if we're computing anti-prominence
   if (mOptions.antiprominence) {
     tile->flipElevations();
@@ -115,11 +116,16 @@ bool ProminenceTask::writeStringToOutputFile(const string &filename, const strin
 
 string ProminenceTask::getFilenamePrefix() const {
   char filename[PATH_MAX];
-  int latHundredths = fractionalDegree(mCurrentLatitude);
-  int lngHundredths = fractionalDegree(mCurrentLongitude);
+
+  // Deal with floating point imprecision
+  float lat = adjustCoordinate(mCurrentLatitude);
+  float lng = adjustCoordinate(mCurrentLongitude);
+
+  int latHundredths = fractionalDegree(lat);
+  int lngHundredths = fractionalDegree(lng);
   snprintf(filename, sizeof(filename), "prominence-%02dx%02d-%03dx%02d",
-           static_cast<int>(mCurrentLatitude), latHundredths,
-           static_cast<int>(mCurrentLongitude), lngHundredths);
+           static_cast<int>(lat), latHundredths,
+           static_cast<int>(lng), lngHundredths);
   return mOptions.outputDir + "/" + filename;
 }
 

code/tile_loading_policy.cpp

@@ -103,6 +103,12 @@ Tile *BasicTileLoadingPolicy::loadTile(float minLat, float minLng) const {
 Tile *BasicTileLoadingPolicy::loadInternal(float minLat, float minLng) const {
   TileLoader *loader = nullptr;
 
+  // If the lat/lng can't be represented exactly in floating point,
+  // there's sometimes roundoff error when converting to integers for
+  // tile filenames.  Adding a little slop prevents truncation.  
+  minLat = adjustCoordinate(minLat);
+  minLng = adjustCoordinate(minLng);
+
   switch (mFileFormat.value()) {
   case FileFormat::Value::HGT:
     loader = new HgtLoader();

code/util.cpp

@@ -44,6 +44,13 @@ float metersToFeet(float meters) {
   return meters / 0.3048f;
 }
 
+float adjustCoordinate(float coordinate) {
+  // A tile is not going to be smaller than 0.1 degrees or so, so this
+  // amount should be safe.
+  const float epsilon = 0.001f;
+  return coordinate + ((coordinate >= 0) ? epsilon : -epsilon);
+}
+
 string trim(const string &s) {
   static const std::string whitespace = " \t\f\v\n\r";
   std::size_t start = s.find_first_not_of(whitespace);

code/util.h

@@ -36,6 +36,10 @@ float metersToFeet(float meters);
 // Trim whitespace from the start and end of the string, return a new string
 std::string trim(const std::string &s);
 
+// Adjust the given coordinate by an epsilon value away from 0, so that truncation
+// to int doesn't give incorrect values to due floating-point imprecision.
+float adjustCoordinate(float coordinate);
+
 /*
  * Split the given string by the given delimiter, putting the
  * result in elems and returning it.