Add "CUSTOM" input format for abitrary resolution.

This format allows the intermediate tile size and resolution
to be specified, so arbitrary data can be reprojected efficiently.
The default is 0.1 degrees per side with 10000 samples per side,
which corresponds roughly to 1m data.

This replaces the "LIDAR" format.

code/file_format.cpp

@@ -25,11 +25,13 @@
 #include "file_format.h"
 #include "degree_coordinate_system.h"
 #include "utm_coordinate_system.h"
+#include "util.h"
 #include "easylogging++.h"
 
 #include <cassert>
 #include <cmath>
 #include <map>
+#include <vector>
 
 using std::string;
 
@@ -43,7 +45,6 @@ int FileFormat::rawSamplesAcross() const {
   case Value::HGT:        return 1201;
   case Value::HGT30:      return 3601;
   case Value::THREEDEP_1M:  return 10012;
-  case Value::LIDAR: return 10000;
   default:
     // In particular, fail on GLO, because this number is variable with latitude.
     LOG(ERROR) << "Couldn't compute tile size of unknown file format";
@@ -65,7 +66,6 @@ int FileFormat::inMemorySamplesAcross() const {
   case Value::GLO30: // Fall through
   case Value::FABDEM:
     return 3601;
-  case Value::LIDAR: return 10000;
   default:
     LOG(ERROR) << "Couldn't compute tile size of unknown file format";
     exit(1);
@@ -86,7 +86,6 @@ double FileFormat::degreesAcross() const {
   case Value::GLO30:  // Fall through
   case Value::FABDEM:
     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).
@@ -120,7 +119,7 @@ CoordinateSystem *FileFormat::coordinateSystemForOrigin(double lat, double lng,
                                       samplesPerDegreeLat, samplesPerDegreeLng);
   }
 
-  case Value::LIDAR: {
+  case Value::CUSTOM: {
     int samplesPerDegreeLat = static_cast<int>(std::round(inMemorySamplesAcross() / degreesAcross()));
     int samplesPerDegreeLng = static_cast<int>(std::round(inMemorySamplesAcross() / degreesAcross()));
     return new DegreeCoordinateSystem(lat, lng,
@@ -129,7 +128,6 @@ CoordinateSystem *FileFormat::coordinateSystemForOrigin(double lat, double lng,
                                       samplesPerDegreeLat, samplesPerDegreeLng);
   }
 
-
   case Value::THREEDEP_1M:
     // 10km x 10km tiles, NW corner
     return new UtmCoordinateSystem(utmZone,
@@ -157,9 +155,40 @@ FileFormat *FileFormat::fromName(const string &name) {
     { "GLO30",     Value::GLO30, },
     { "FABDEM",    Value::FABDEM, },
     { "3DEP-1M",   Value::THREEDEP_1M, },
-    { "LIDAR",     Value::LIDAR, },
   };
 
+  // Handle CUSTOM-<degrees per tile>-<samples across>
+  if (0 == name.compare(0, 6, "CUSTOM")) {
+    std::vector<string> fields;
+    split(name, '-', fields);
+    if (fields.size() != 3) {
+      LOG(ERROR) << "Custom format must have 3 components";
+      return nullptr;
+    }
+
+    double degreesAcross = std::stod(fields[1]);
+    int samplesAcross = std::stoi(fields[2]);
+
+    if (degreesAcross <= 0) {
+      LOG(ERROR) << "Illegal value for degrees per tile: " << degreesAcross;
+      return nullptr;
+    }
+
+    // degreesAcross must divide 1
+    double tilesPerDegree = 1.0 / degreesAcross;
+    if (tilesPerDegree - static_cast<int>(tilesPerDegree) > 0.001) {
+      LOG(ERROR) << "Value for degrees per tile must divide 1 evenly: " << degreesAcross;
+      return nullptr;
+    }
+
+    if (samplesAcross <= 0) {
+      LOG(ERROR) << "Illegal value for samples across tile: " << samplesAcross;
+      return nullptr;
+    }
+
+    return new CustomFileFormat(degreesAcross, samplesAcross);
+  }
+
   auto it = fileFormatNames.find(name);
   if (it == fileFormatNames.end()) {
     return nullptr;

code/file_format.h

@@ -43,7 +43,7 @@ class FileFormat {
     THREEDEP_1M,  // FLT file containing one-meter LIDAR from 3D Elevation Program (3DEP)
     GLO30,  // Copernicus GLO-30 30m data
     FABDEM, // Tree-free Copernicus GLO-30 30m data
-    LIDAR,  // FLT file converted from LIDAR, 0.1 arcsecond
+    CUSTOM, // FLT file with customizable resolution
   };
 
   FileFormat() = default;
@@ -53,16 +53,16 @@ class FileFormat {
 
   // Return the number of samples in one row or column of the file format,
   // including any border samples.
-  int rawSamplesAcross() const;
+  virtual int rawSamplesAcross() const;
 
   // Return the number of samples in a tile after its loaded and the borders
   // have been modified.
-  int inMemorySamplesAcross() const;
+  virtual int inMemorySamplesAcross() const;
 
   // 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.
-  double degreesAcross() const;
+  virtual double degreesAcross() const;
 
   // Does this format use UTM coordinates rather than lat/lng?
   bool isUtm() const;
@@ -79,4 +79,32 @@ class FileFormat {
   Value mValue;
 };
 
+// A file format where the size of the tile in degrees and samples
+// can be specified.  There must be an integer number of tiles per
+// degree.  The data is assumed to be in FLT format.
+class CustomFileFormat : public FileFormat {
+public:
+  CustomFileFormat(double degreesAcross, int samplesAcross) :
+      FileFormat(Value::CUSTOM),
+      mDegreesAcross(degreesAcross),
+      mSamplesAcross(samplesAcross) {
+  }
+
+  virtual int rawSamplesAcross() const {
+    return mSamplesAcross;
+  }
+
+  virtual int inMemorySamplesAcross() const {
+    return mSamplesAcross;
+  }
+
+  virtual double degreesAcross() const {
+    return mDegreesAcross;
+  }
+
+private:
+  double mDegreesAcross;
+  int mSamplesAcross;
+};
+
 #endif  // _FILE_FORMAT_H_

code/flt_loader.cpp

@@ -33,9 +33,9 @@ using std::string;
 // while NED19 uses some very negative value (< -1e38).
 static const float NED_NODATA_MIN_ELEVATION = -9998;
 
-FltLoader::FltLoader(const FileFormat &format, int utmZone) {
-  mFormat = format;
-  mUtmZone = utmZone;
+FltLoader::FltLoader(const FileFormat &format, int utmZone) :
+    mFormat(format),
+    mUtmZone(utmZone) {
 }
 
 Tile *FltLoader::loadTile(const std::string &directory, double minLat, double minLng) {
@@ -47,7 +47,7 @@ Tile *FltLoader::loadTile(const std::string &directory, double minLat, double mi
   case FileFormat::Value::NED13:
   case FileFormat::Value::NED19:
   case FileFormat::Value::THREEDEP_1M:
-  case FileFormat::Value::LIDAR:
+  case FileFormat::Value::CUSTOM:
     return loadFromFltFile(directory, minLat, minLng);
 
   default:
@@ -194,7 +194,7 @@ string FltLoader::getFltFilename(double minLat, double minLng, const FileFormat
              fractionalDegree(minLng));
     break;
 
-  case FileFormat::Value::LIDAR:
+  case FileFormat::Value::CUSTOM:
     snprintf(buf, sizeof(buf), "tile_%02dx%02d_%03dx%02d.flt",
              static_cast<int>(upperLat),
              fractionalDegree(upperLat),

code/flt_loader.h

@@ -43,7 +43,7 @@ class FltLoader : public TileLoader {
   virtual Tile *loadTile(const std::string &directory, double minLat, double minLng);
 
 private:
-  FileFormat mFormat;
+  const FileFormat &mFormat;
   int mUtmZone;  // For data in UTM coordinates
 
   Tile *loadFromNEDZipFileInternal(const std::string &directory, double minLat, double minLng);

code/hgt_loader.cpp

@@ -38,7 +38,7 @@ static uint16 swapByteOrder16(uint16 us) {
   return (us >> 8) | (us << 8);
 }
 
-HgtLoader::HgtLoader(FileFormat fileFormat) {
+HgtLoader::HgtLoader(const FileFormat &fileFormat) {
   switch (fileFormat.value()) {
   case FileFormat::Value::HGT:
     mTileSize = 1201;

code/hgt_loader.h

@@ -32,7 +32,7 @@
 
 class HgtLoader : public TileLoader {
 public:
-  explicit HgtLoader(FileFormat fileFormat);
+  explicit HgtLoader(const FileFormat &fileFormat);
 
   // minLat and minLng name the SW corner of the tile, in degrees
   virtual Tile *loadTile(const std::string &directory, double minLat, double minLng);

code/isolation_finder.cpp

@@ -36,11 +36,11 @@
 using std::vector;
 
 IsolationFinder::IsolationFinder(TileCache *cache, const Tile *tile,
-                                 const CoordinateSystem &coordinateSystem, FileFormat format) {
+                                 const CoordinateSystem &coordinateSystem, const FileFormat &format) :
+    mFormat(format) {
   mTile = tile;
   mCache = cache;
   mCoordinateSystem = std::unique_ptr<CoordinateSystem>(coordinateSystem.clone());
-  mFormat = format;
 }
 
 IsolationRecord IsolationFinder::findIsolation(Offsets peak) const {

code/isolation_finder.h

@@ -60,7 +60,7 @@ struct IsolationRecord {
 class IsolationFinder {
 public:
   IsolationFinder(TileCache *cache, const Tile *tile,
-                  const CoordinateSystem &coordinateSystem, FileFormat format);
+                  const CoordinateSystem &coordinateSystem, const FileFormat &format);
 
   IsolationRecord findIsolation(Offsets peak) const;
 
@@ -69,7 +69,7 @@ class IsolationFinder {
   const Tile *mTile;
   TileCache *mCache;
   std::unique_ptr<CoordinateSystem> mCoordinateSystem;
-  FileFormat mFormat;
+  const FileFormat &mFormat;
 
   // Search tile for a point higher than seedElevation.
   //

code/prominence.cpp

@@ -73,7 +73,7 @@ int main(int argc, char **argv) {
 
   Elevation minProminence = 100;
   int numThreads = 1;
-  FileFormat fileFormat = FileFormat(FileFormat::Value::HGT);
+  std::unique_ptr<FileFormat> fileFormat(new FileFormat(FileFormat::Value::HGT));
 
   // Parse options
   START_EASYLOGGINGPP(argc, argv);
@@ -104,7 +104,7 @@ int main(int argc, char **argv) {
         usage();
       }
 
-      fileFormat = *format;
+      fileFormat.reset(format);
       break;
     }
 
@@ -142,7 +142,7 @@ int main(int argc, char **argv) {
     usage();
   }
 
-  if (fileFormat.isUtm() && utmZone == NO_UTM_ZONE) {
+  if (fileFormat->isUtm() && utmZone == NO_UTM_ZONE) {
     LOG(ERROR) << "You must specify a UTM zone with this format";
     exit(1);
   }
@@ -158,7 +158,7 @@ int main(int argc, char **argv) {
   }
 
   // Validate that bounds are on tile boundaries
-  double degreesAcross = fileFormat.degreesAcross();
+  double degreesAcross = fileFormat->degreesAcross();
   for (auto bound : bounds) {
     if (fabs(bound / degreesAcross - static_cast<int>(std::round(bound / degreesAcross))) > 0.001) {
       LOG(ERROR) << "Coordinates must be multiples of " << degreesAcross;
@@ -170,7 +170,7 @@ int main(int argc, char **argv) {
   // Load filtering polygon
   Filter filter;
   if (!polygonFilename.empty()) {
-    if (fileFormat.isUtm()) {
+    if (fileFormat->isUtm()) {
       LOG(ERROR) << "Can't specify a filter polygon with UTM data";
       exit(1);
     }
@@ -186,7 +186,7 @@ int main(int argc, char **argv) {
                                static_cast<bool>(writeKml)};
 
   // Caching doesn't do anything for our calculation and the tiles are huge
-  BasicTileLoadingPolicy policy(terrain_directory, fileFormat);
+  BasicTileLoadingPolicy policy(terrain_directory, *fileFormat.get());
   policy.enableNeighborEdgeLoading(true);
   if (utmZone != NO_UTM_ZONE) {
     policy.setUtmZone(utmZone);
@@ -210,16 +210,16 @@ int main(int argc, char **argv) {
     while (lng >= bounds[2]) {
       // Allow specifying longitude ranges that span the antimeridian (lng > 180)
       auto wrappedLng = lng;
-      if (!fileFormat.isUtm() && wrappedLng >= 180) {
+      if (!fileFormat->isUtm() && wrappedLng >= 180) {
         wrappedLng -= 360;
       }
 
       std::shared_ptr<CoordinateSystem> coordinateSystem(
-        fileFormat.coordinateSystemForOrigin(lat, wrappedLng, utmZone));
+        fileFormat->coordinateSystemForOrigin(lat, wrappedLng, utmZone));
 
       // Skip tiles that don't intersect filtering polygon
-      if (!filter.intersects(lat, lat + fileFormat.degreesAcross(),
-                             lng, lng + fileFormat.degreesAcross())) {
+      if (!filter.intersects(lat, lat + fileFormat->degreesAcross(),
+                             lng, lng + fileFormat->degreesAcross())) {
         VLOG(3) << "Skipping tile that doesn't intersect polygon " << lat << " " << lng;
       } else {
         // Actually calculate prominence
@@ -229,9 +229,9 @@ int main(int argc, char **argv) {
             }));
       }
 
-      lng -= fileFormat.degreesAcross();
+      lng -= fileFormat->degreesAcross();
     }
-    lat += fileFormat.degreesAcross();
+    lat += fileFormat->degreesAcross();
   }
 
   for (auto && result : results) {

code/tile_loading_policy.cpp

@@ -87,7 +87,7 @@ Tile *BasicTileLoadingPolicy::loadTile(double minLat, double minLng) const {
       break;
 
     case FileFormat::Value::GLO30:  // fall through
-    case FileFormat::Value::LIDAR:
+    case FileFormat::Value::CUSTOM:
     case FileFormat::Value::FABDEM: {
       // GLO30 "helpfully" removes the last row and column from each tile,
       // so we need to stick them back on.
@@ -126,7 +126,7 @@ Tile *BasicTileLoadingPolicy::loadInternal(double minLat, double minLng) const {
   case FileFormat::Value::NED19:
   case FileFormat::Value::NED1_ZIP:
   case FileFormat::Value::THREEDEP_1M:
-  case FileFormat::Value::LIDAR:
+  case FileFormat::Value::CUSTOM:
     loader = new FltLoader(mFileFormat, mUtmZone);
     break;
 

code/tile_loading_policy.h

@@ -69,7 +69,7 @@ class BasicTileLoadingPolicy : public TileLoadingPolicy {
 
 private:
   std::string mDirectory;  // Directory for loading tiles
-  FileFormat mFileFormat;  
+  const FileFormat &mFileFormat;  
   bool mNeighborEdgeLoadingEnabled;
   int mUtmZone;
   TileCache *mTileCache;