From d0bafbba230a27510a5f30690dbdd05da6b08582 Mon Sep 17 00:00:00 2001 From: Alvaro Huarte Date: Wed, 23 May 2018 22:49:25 +0200 Subject: [PATCH] Create missing root tiles from scratch --- tools/ctb-tile.cpp | 94 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 86 insertions(+), 8 deletions(-) diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp index 2a9e8cb..1b5f78d 100644 --- a/tools/ctb-tile.cpp +++ b/tools/ctb-tile.cpp @@ -255,9 +255,9 @@ class TerrainBuild : public Command { * to ensure all tiles are iterated over consecutively. It assumes individual * tile iterators point to the same source GDAL dataset. */ +static int globalIteratorIndex = 0; // keep track of where we are globally template int incrementIterator(T &iter, int currentIndex) { - static int globalIteratorIndex = 0; // keep track of where we are globally static mutex mutex; // ensure iterations occur serially between threads lock_guard lock(mutex); @@ -501,6 +501,71 @@ class TerrainMetadata { } }; +/// Create an empty root temporary elevation file (GTiff) +static std::string +createEmptyRootElevationFile(std::string &fileName, const Grid &grid, const TileCoordinate& coord) { + GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + + if (poDriver == NULL) { + throw CTBException("Could not retrieve GTiff GDAL driver"); + } + + // Create the geo transform for this temporary elevation tile. + // We apply an 1-degree negative offset to avoid problems in borders. + CRSBounds tileBounds = grid.tileBounds(coord); + tileBounds.setMinX(tileBounds.getMinX() + 1); + tileBounds.setMinY(tileBounds.getMinY() + 1); + tileBounds.setMaxX(tileBounds.getMaxX() - 1); + tileBounds.setMaxY(tileBounds.getMaxY() - 1); + const i_tile tileSize = grid.tileSize() - 2; + const double resolution = tileBounds.getWidth() / tileSize; + double adfGeoTransform[6] = { tileBounds.getMinX(), resolution, 0, tileBounds.getMaxY(), 0, -resolution }; + + // Create the spatial reference system for the file + OGRSpatialReference oSRS; + if (oSRS.importFromEPSG(4326) != OGRERR_NONE) { + throw CTBException("Could not create EPSG:4326 spatial reference"); + } + char *pszDstWKT = NULL; + if (oSRS.exportToWkt(&pszDstWKT) != OGRERR_NONE) { + CPLFree(pszDstWKT); + throw CTBException("Could not create EPSG:4326 WKT string"); + } + + // Create the GTiff file + fileName += ".tif"; + GDALDataset *poDataset = poDriver->Create(fileName.c_str(), tileSize, tileSize, 1, GDT_Float32, NULL); + + // Set the projection + if (poDataset->SetProjection(pszDstWKT) != CE_None) { + poDataset->Release(); + CPLFree(pszDstWKT); + throw CTBException("Could not set projection on temporary elevation file"); + } + CPLFree(pszDstWKT); + + // Apply the geo transform + if (poDataset->SetGeoTransform(adfGeoTransform) != CE_None) { + poDataset->Release(); + throw CTBException("Could not set projection on temporary elevation file"); + } + + // Finally write the height data + float *rasterHeights = (float *)CPLCalloc(tileSize * tileSize, sizeof(float)); + GDALRasterBand *heightsBand = poDataset->GetRasterBand(1); + if (heightsBand->RasterIO(GF_Write, 0, 0, tileSize, tileSize, + (void *)rasterHeights, tileSize, tileSize, GDT_Float32, + 0, 0) != CE_None) { + CPLFree(rasterHeights); + throw CTBException("Could not write heights on temporary elevation file"); + } + CPLFree(rasterHeights); + + poDataset->FlushCache(); + poDataset->Release(); + return fileName; +} + /// Output GDAL tiles represented by a tiler to a directory static void buildGDAL(GDALSerializer &serializer, const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) { @@ -638,8 +703,8 @@ buildMetadata(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata * * This function is designed to be run in a separate thread. */ static int -runTiler(TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { - GDALDataset *poDataset = (GDALDataset *) GDALOpen(command->getInputFilename(), GA_ReadOnly); +runTiler(const char *inputFilename, TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) { + GDALDataset *poDataset = (GDALDataset *) GDALOpen(inputFilename, GA_ReadOnly); if (poDataset == NULL) { cerr << "Error: could not open GDAL dataset" << endl; return 1; @@ -757,9 +822,9 @@ main(int argc, char *argv[]) { // Instantiate the threads using futures from a packaged_task for (int i = 0; i < threadCount ; ++i) { - packaged_task task(runTiler); // wrap the function - tasks.push_back(task.get_future()); // get a future - thread(move(task), &command, &grid, metadata).detach(); // launch on a thread + packaged_task task(runTiler); // wrap the function + tasks.push_back(task.get_future()); // get a future + thread(move(task), command.getInputFilename(), &command, &grid, metadata).detach(); // launch on a thread } // Synchronise the completion of the threads @@ -785,14 +850,27 @@ main(int argc, char *argv[]) { std::string tileName0 = dirName0 + osDirSep + "0.terrain"; std::string tileName1 = dirName1 + osDirSep + "0.terrain"; + i_zoom missing_zoom = 65535; + ctb::TileCoordinate missingTileCoord(missing_zoom, 0, 0); + std::string missingTileName; + if (fileExists(tileName0) && !fileExists(tileName1)) { VSIMkdir(dirName1.c_str(), 0755); - fileCopy(tileName0, tileName1); + missingTileCoord = ctb::TileCoordinate(0, 1, 0); + missingTileName = tileName1; } else if (!fileExists(tileName0) && fileExists(tileName1)) { VSIMkdir(dirName0.c_str(), 0755); - fileCopy(tileName1, tileName0); + missingTileCoord = ctb::TileCoordinate(0, 0, 0); + missingTileName = tileName0; + } + if (missingTileCoord.zoom != missing_zoom) { + missingTileName = createEmptyRootElevationFile(missingTileName, grid, missingTileCoord); + command.startZoom = command.endZoom = 0; + globalIteratorIndex = 0; // reset global iterator index + runTiler (missingTileName.c_str(), &command, &grid, NULL); + VSIUnlink(missingTileName.c_str()); } // Fix available indexes.