Convert lat/lng from float to double. (#29)

code/degree_coordinate_system.cpp

@@ -33,7 +33,7 @@
 using std::string;
 using std::vector;
 
-DegreeCoordinateSystem::DegreeCoordinateSystem(float minLat, float minLng, float maxLat, float maxLng,
+DegreeCoordinateSystem::DegreeCoordinateSystem(double minLat, double minLng, double maxLat, double maxLng,
                                                int pixelsPerDegreeLat, int pixelsPerDegreeLng) {
   mMinLatitude = minLat;
   mMinLongitude = minLng;
@@ -66,10 +66,10 @@ bool DegreeCoordinateSystem::compatibleWith(const CoordinateSystem &that) const
 
 LatLng DegreeCoordinateSystem::getLatLng(Offsets offsets) const {
   // Positive y is south
-  float latitude = mMaxLatitude -
-    ((float) offsets.y()) / mSamplesPerDegreeLatitude;
-  float longitude = mMinLongitude +
-    ((float) offsets.x()) / mSamplesPerDegreeLongitude;
+  double latitude = mMaxLatitude -
+    ((double) offsets.y()) / mSamplesPerDegreeLatitude;
+  double longitude = mMinLongitude +
+    ((double) offsets.x()) / mSamplesPerDegreeLongitude;
   return LatLng(latitude, longitude);
 }
 
@@ -119,7 +119,7 @@ string DegreeCoordinateSystem::toString() const {
 }
 
 CoordinateSystem *DegreeCoordinateSystem::fromString(const string &str) {
-  float minLat = 0, minLng = 0, maxLat = 0, maxLng = 0;
+  double minLat = 0, minLng = 0, maxLat = 0, maxLng = 0;
   int samplesPerLat = 0, samplesPerLng = 0;
 
   vector<string> elements;
@@ -129,14 +129,14 @@ CoordinateSystem *DegreeCoordinateSystem::fromString(const string &str) {
     return nullptr;
   }
 
-  minLat = stof(elements[1]);
-  minLng = stof(elements[2]);
+  minLat = stod(elements[1]);
+  minLng = stod(elements[2]);
   samplesPerLat = stoi(elements[3]);
   samplesPerLng = stoi(elements[4]);
   // Max lat/lng was added later for non-1x1 tile support
   if (elements.size() >= 7) {
-    maxLat = stof(elements[5]);
-    maxLng = stof(elements[6]);
+    maxLat = stod(elements[5]);
+    maxLng = stod(elements[6]);
   } else {
     // Assume 1x1 tile for old DVT files
     maxLat = minLat + 1;

code/degree_coordinate_system.h

@@ -33,7 +33,7 @@
 class DegreeCoordinateSystem : public CoordinateSystem {
 public:
 
-  DegreeCoordinateSystem(float minLat, float minLng, float maxLat, float maxLng,
+  DegreeCoordinateSystem(double minLat, double minLng, double maxLat, double maxLng,
                          int pixelsPerDegreeLat, int pixelsPerDegreeLng);
   virtual ~DegreeCoordinateSystem();
 
@@ -55,10 +55,10 @@ class DegreeCoordinateSystem : public CoordinateSystem {
   static CoordinateSystem *fromString(const std::string &str);
 
 private:
-  float mMinLatitude;
-  float mMinLongitude;
-  float mMaxLatitude;
-  float mMaxLongitude;
+  double mMinLatitude;
+  double mMinLongitude;
+  double mMaxLatitude;
+  double mMaxLongitude;
   int mSamplesPerDegreeLatitude;
   int mSamplesPerDegreeLongitude;
 };

code/fabdem_loader.cpp

@@ -36,7 +36,7 @@ using std::string;
 
 static const float COPERNICUS_NODATA_ELEVATION = -32767.0f;
 
-Tile *FabdemLoader::loadTile(const std::string &directory, float minLat, float minLng) {
+Tile *FabdemLoader::loadTile(const std::string &directory, double minLat, double minLng) {
   char buf[100];
   snprintf(buf, sizeof(buf), "%c%02d%c%03d_FABDEM_V1-0.flt",
            (minLat >= 0) ? 'N' : 'S',

code/fabdem_loader.h

@@ -35,7 +35,7 @@
 class FabdemLoader : public TileLoader {
 public:
   // minLat and minLng name the SW corner of the tile, in degrees
-  virtual Tile *loadTile(const std::string &directory, float minLat, float minLng);
+  virtual Tile *loadTile(const std::string &directory, double minLat, double minLng);
 
 private:
 };

code/file_format.cpp

@@ -70,22 +70,22 @@ int FileFormat::inMemorySamplesAcross() const {
 }
 
 
-float FileFormat::degreesAcross() const {
+double FileFormat::degreesAcross() const {
   switch (mValue) {
   case Value::NED13:  // fall through
   case Value::NED13_ZIP:
-    return 1.0f;
-  case Value::NED1_ZIP:   return 1.0f;
-  case Value::NED19:      return 0.25f;
-  case Value::HGT:        return 1.0f;
+    return 1.0;
+  case Value::NED1_ZIP:   return 1.0;
+  case Value::NED19:      return 0.25;
+  case Value::HGT:        return 1.0;
   case Value::GLO30:  // Fall through
   case Value::FABDEM:
-    return 1.0f;
-  case Value::LIDAR: return 0.1f;
+    return 1.0;
+  case Value::LIDAR: return 0.1;
   case Value::THREEDEP_1M:
     // This is a misnomer, as these tiles are in UTM coordinates.  The "degrees" across
     // means one x or y unit per tile (where each tile is 10000m in UTM).
-    return 1.0f;
+    return 1.0;
   default:
     LOG(ERROR) << "Couldn't compute degree span of unknown file format";
     exit(1);
@@ -96,7 +96,7 @@ bool FileFormat::isUtm() const {
   return mValue == Value::THREEDEP_1M;
 }
 
-CoordinateSystem *FileFormat::coordinateSystemForOrigin(float lat, float lng, int utmZone) const {
+CoordinateSystem *FileFormat::coordinateSystemForOrigin(double lat, double lng, int utmZone) const {
   switch (mValue) {
   case Value::NED13_ZIP:  // fall through
   case Value::NED13:
@@ -131,7 +131,7 @@ CoordinateSystem *FileFormat::coordinateSystemForOrigin(float lat, float lng, in
                                    static_cast<int>((lat - 1) * 10000),
                                    static_cast<int>((lng + 1) * 10000),
                                    static_cast<int>(lat * 10000),
-                                   1.0f);
+                                   1.0);
 
   default:
     assert(false);

code/file_format.h

@@ -61,14 +61,14 @@ class FileFormat {
   // Return the degrees in lat or lng covered by one tile.
   // Note that this is the logical value (1 degree, 0.25 degree), not necessarily
   // the precise value covered, including border samples.
-  float degreesAcross() const;
+  double degreesAcross() const;
 
   // Does this format use UTM coordinates rather than lat/lng?
   bool isUtm() const;
 
   // Return a new CoordinateSystem describing the section of the Earth that
   // the given tile with the given origin (lower-left corner) covers.
-  CoordinateSystem *coordinateSystemForOrigin(float lat, float lng, int utmZone = 0) const;
+  CoordinateSystem *coordinateSystemForOrigin(double lat, double lng, int utmZone = 0) const;
 
   // Return a FileFormat object for the given human-readable string,
   // or nullptr if none.

code/filter.cpp

@@ -73,8 +73,8 @@ bool Filter::addPolygonsFromKml(const string &filename) {
       for (const string &eachPoint : pointText) {
         split(eachPoint, ',', elements);
         if (elements.size() >= 2) {
-          float lat = stof(elements[1]);  // note order in KML
-          float lng = stof(elements[0]);
+          double lat = stod(elements[1]);  // note order in KML
+          double lng = stod(elements[0]);
           polygon.push_back(LatLng(lat, lng));
         }
       }
@@ -109,8 +109,8 @@ bool Filter::isPointInside(const LatLng &latlng) const {
     }
 
     bool inside = false;
-    float testx = latlng.longitude();
-    float testy = latlng.latitude();
+    auto testx = latlng.longitude();
+    auto testy = latlng.latitude();
 
     // Allow wrapping around antimeridian
     if (testx < mWrapLongitude) {
@@ -134,7 +134,7 @@ bool Filter::isPointInside(const LatLng &latlng) const {
   return false;
 }
 
-void Filter::setWrapLongitude(float wrapLongitude) {
+void Filter::setWrapLongitude(double wrapLongitude) {
   mWrapLongitude = wrapLongitude;
   vector<vector<LatLng>> newPolygons;
   for (auto &polygon : mPolygons) {
@@ -151,7 +151,7 @@ void Filter::setWrapLongitude(float wrapLongitude) {
   mPolygons = newPolygons;
 }
 
-bool Filter::intersects(float minLat, float maxLat, float minLng, float maxLng) const {
+bool Filter::intersects(double minLat, double maxLat, double minLng, double maxLng) const {
   // If any corner of the rectangle is inside, there is an intersection
   LatLng p1(minLat, minLng);
   LatLng p2(minLat, maxLng);
@@ -198,10 +198,10 @@ void Filter::getBounds(LatLng *sw, LatLng *ne) const {
 
   bool anyPoints = false;
 
-  float min_latitude = 999;
-  float max_latitude = -999;
-  float min_longitude = 999;
-  float max_longitude = -999;
+  double min_latitude = 999;
+  double max_latitude = -999;
+  double min_longitude = 999;
+  double max_longitude = -999;
 
   for (auto &polygon : mPolygons) {
     if (polygon.empty()) {
@@ -221,14 +221,14 @@ void Filter::getBounds(LatLng *sw, LatLng *ne) const {
 }
 
 // Returns true if the lines intersect, otherwise false.
-bool Filter::segmentsIntersect(float p0_x, float p0_y, float p1_x, float p1_y,
-                               float p2_x, float p2_y, float p3_x, float p3_y) const {
+bool Filter::segmentsIntersect(double p0_x, double p0_y, double p1_x, double p1_y,
+                               double p2_x, double p2_y, double p3_x, double p3_y) const {
   // See http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
-  float s1_x, s1_y, s2_x, s2_y;
+  double s1_x, s1_y, s2_x, s2_y;
   s1_x = p1_x - p0_x;     s1_y = p1_y - p0_y;
   s2_x = p3_x - p2_x;     s2_y = p3_y - p2_y;
 
-  float s, t;
+  double s, t;
   s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / (-s2_x * s1_y + s1_x * s2_y);
   t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / (-s2_x * s1_y + s1_x * s2_y);
 

code/filter.h

@@ -50,21 +50,21 @@ class Filter {
   // Add 360 to the longitude of (presumably negative) longitudes <
   // wrapLongitude in the polygons.  This allows the polygons to cross
   // the antimeridian in a primitive way.
-  void setWrapLongitude(float wrapLongitude);
+  void setWrapLongitude(double wrapLongitude);
 
   // Determines whether given rectangle intersects any polygon.
-  bool intersects(float minLat, float maxLat, float minLng, float maxLng) const;
+  bool intersects(double minLat, double maxLat, double minLng, double maxLng) const;
 
   // Sets sw and ne to northwest and southeast corners of bounds.
   void getBounds(LatLng *sw, LatLng *ne) const;
 
 private:
   std::vector<std::vector<LatLng>> mPolygons;
-  float mWrapLongitude;
+  double mWrapLongitude;
 
   // Returns true if line segment p0-p1 intersects with segment p2-p3.
-  bool segmentsIntersect(float p0_x, float p0_y, float p1_x, float p1_y,
-                         float p2_x, float p2_y, float p3_x, float p3_y) const;  
+  bool segmentsIntersect(double p0_x, double p0_y, double p1_x, double p1_y,
+                         double p2_x, double p2_y, double p3_x, double p3_y) const;  
 };
 
 #endif  // _FILTER_H_

code/filter_points.cpp

@@ -52,7 +52,7 @@ static void usage() {
 }
 
 int main(int argc, char **argv) {
-  float wrapLongitude = -180;
+  double wrapLongitude = -180;
 
   // Parse options
   START_EASYLOGGINGPP(argc, argv);
@@ -66,7 +66,7 @@ int main(int argc, char **argv) {
   while ((ch = getopt_long(argc, argv, "a:", long_options, nullptr)) != -1) {
     switch (ch) {
     case 'a':
-      wrapLongitude = static_cast<float>(atof(optarg));
+      wrapLongitude = atof(optarg);
       break;
     }
   }
@@ -111,8 +111,8 @@ int main(int argc, char **argv) {
     }
 
     split(line, ',', elements);
-    float lat = stof(elements[0]);
-    float lng = stof(elements[1]);
+    double lat = stod(elements[0]);
+    double lng = stod(elements[1]);
     LatLng point(lat, lng);
 
     if (filter.isPointInside(point)) {

code/flt_loader.cpp

@@ -38,7 +38,7 @@ FltLoader::FltLoader(const FileFormat &format, int utmZone) {
   mUtmZone = utmZone;
 }
 
-Tile *FltLoader::loadTile(const std::string &directory, float minLat, float minLng) {
+Tile *FltLoader::loadTile(const std::string &directory, double minLat, double minLng) {
   switch (mFormat.value()) {
   case FileFormat::Value::NED13_ZIP:  // fall through
   case FileFormat::Value::NED1_ZIP:
@@ -56,7 +56,7 @@ Tile *FltLoader::loadTile(const std::string &directory, float minLat, float minL
   }
 }
 
-Tile *FltLoader::loadFromFltFile(const string &directory, float minLat, float minLng) {
+Tile *FltLoader::loadFromFltFile(const string &directory, double minLat, double minLng) {
   string filename = getFltFilename(minLat, minLng, mFormat);
   if (!directory.empty()) {
     filename = directory + "/" + filename;
@@ -117,7 +117,7 @@ Tile *FltLoader::loadFromFltFile(const string &directory, float minLat, float mi
 }
 
 Tile *FltLoader::loadFromNEDZipFileInternal(const std::string &directory,
-                                            float minLat, float minLng) {
+                                            double minLat, double minLng) {
   // ZIP formats come only in 1x1 degree formats, so OK to cast lat/lng to int.
   // NED uses upper left corner for naming.
   char buf[100];
@@ -161,10 +161,10 @@ Tile *FltLoader::loadFromNEDZipFileInternal(const std::string &directory,
   return tile;  
 }
 
-string FltLoader::getFltFilename(float minLat, float minLng, const FileFormat &format) {
+string FltLoader::getFltFilename(double minLat, double minLng, const FileFormat &format) {
   char buf[100];
   // NED uses upper left corner for naming
-  float upperLat = minLat + format.degreesAcross();
+  auto upperLat = minLat + format.degreesAcross();
   switch (format.value()) {
   case FileFormat::Value::NED13_ZIP:  // fall through
   case FileFormat::Value::NED1_ZIP:
@@ -216,7 +216,7 @@ string FltLoader::getFltFilename(float minLat, float minLng, const FileFormat &f
   return buf;
 }
 
-int FltLoader::fractionalDegree(float degree) const {
-  float excess = fabs(degree - static_cast<int>(degree));
+int FltLoader::fractionalDegree(double degree) const {
+  double excess = fabs(degree - static_cast<int>(degree));
   return static_cast<int>(std::round(100 * excess));
 }

code/flt_loader.h

@@ -40,21 +40,21 @@ class FltLoader : public TileLoader {
   // the latitude will be treated as northing, and the longitude as easting.
   explicit FltLoader(const FileFormat &format, int utmZone);
 
-  virtual Tile *loadTile(const std::string &directory, float minLat, float minLng);
+  virtual Tile *loadTile(const std::string &directory, double minLat, double minLng);
 
 private:
   FileFormat mFormat;
   int mUtmZone;  // For data in UTM coordinates
 
-  Tile *loadFromNEDZipFileInternal(const std::string &directory, float minLat, float minLng);
+  Tile *loadFromNEDZipFileInternal(const std::string &directory, double minLat, double minLng);
 
-  Tile *loadFromFltFile(const std::string &directory, float minLat, float minLng);
+  Tile *loadFromFltFile(const std::string &directory, double minLat, double minLng);
 
   // Return the filename for the .flt file for the given coordinates
-  std::string getFltFilename(float minLat, float minLng, const FileFormat &format);
+  std::string getFltFilename(double minLat, double minLng, const FileFormat &format);
 
   // Return hundredths of a degree from the given value
-  int fractionalDegree(float degree) const;
+  int fractionalDegree(double degree) const;
 };
 
 #endif  // _FLOT_LOADER_H_

code/glo_loader.cpp

@@ -38,7 +38,7 @@ static const float COPERNICUS_NODATA_ELEVATION = -32767.0f;
 // For ALOS 30 tiles impersonating GLO 30 tiles where GLO 30 doesn't have coverage
 static const float ALOSWORLD3D_NODATA_ELEVATION = -9999.0f;
 
-Tile *GloLoader::loadTile(const std::string &directory, float minLat, float minLng) {
+Tile *GloLoader::loadTile(const std::string &directory, double minLat, double minLng) {
   char buf[100];
   snprintf(buf, sizeof(buf), "Copernicus_DSM_COG_10_%c%02d_00_%c%03d_00_DEM.flt",
            (minLat >= 0) ? 'N' : 'S',
@@ -136,7 +136,7 @@ Tile *GloLoader::loadTile(const std::string &directory, float minLat, float minL
   return tile;  
 }
 
-int GloLoader::getWidthForLatitude(float minLat) const {
+int GloLoader::getWidthForLatitude(double minLat) const {
   // See table at
   // https://copernicus-dem-30m.s3.amazonaws.com/readme.html
   if (minLat >= 85 || minLat < -85) {

code/glo_loader.h

@@ -34,11 +34,11 @@
 class GloLoader : public TileLoader {
 public:
   // minLat and minLng name the SW corner of the tile, in degrees
-  virtual Tile *loadTile(const std::string &directory, float minLat, float minLng);
+  virtual Tile *loadTile(const std::string &directory, double minLat, double minLng);
 
 private:
   // Width in samples for a tile at the given latitude
-  int getWidthForLatitude(float minLat) const;
+  int getWidthForLatitude(double minLat) const;
 };
 
 #endif  // _GLO_LOADER_H_