Cecisum terrain builder的错误修正

从.terrain->.tif

以及从.tif->.terrain

 

，但是.tif格式是特定的，是浮点型，所以对读取不同类型的.tif就要先判断.tif格式，再读取。

.terrain是16位整型，这个是不变的。

 

具体如下：

 

1,读.tif时，在terrainTiler.cpp中改写以下函数，这样可以读取任意格式的TIF

TerrainTile *
ctb::TerrainTiler::createTile(const TileCoordinate &coord) const {
    // Get a terrain tile represented by the tile coordinate
    TerrainTile *terrainTile = new TerrainTile(coord);
    GDALTile *rasterTile = createRasterTile(coord); // the raster associated with this tile coordinate
    GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1);
    GDALDataType theType = heightsBand->GetRasterDataType();
    //这里将栅格按照浮点型拷贝，栅格是否是浮点型？有问题
    // Copy the raster data into an array
    float rasterHeights[TerrainTile::TILE_CELL_SIZE];
    if (heightsBand->RasterIO(GF_Read, 0, 0, TILE_SIZE, TILE_SIZE,
        (void *)rasterHeights, TILE_SIZE, TILE_SIZE, theType,
        0, 0) != CE_None)
        /*
        if (heightsBand->RasterIO(GF_Read, 0, 0, TILE_SIZE, TILE_SIZE,
            (void *)rasterHeights, TILE_SIZE, TILE_SIZE, GDT_Float32,
            0, 0) != CE_None)
            */
    {
        throw CTBException("Could not read heights from raster");
    }

    delete rasterTile;

    // Convert the raster data into the terrain tile heights.  This assumes the
    // input raster data represents meters above sea level. Each terrain height
    // value is the number of 1/5 meter units above -1000 meters.
    // TODO: try doing this using a VRT derived band:
    // (http://www.gdal.org/gdal_vrttut.html)
    for (unsigned short int i = 0; i < TerrainTile::TILE_CELL_SIZE; i++)
    {
        terrainTile->mHeights[i] = (i_terrain_height)((rasterHeights[i] + 1000) * 5); （高程值公式)
     
    }

    // If we are not at the maximum zoom level we need to set child flags on the
    // tile where child tiles overlap the dataset bounds.
    if (coord.zoom != maxZoomLevel()) {
        CRSBounds tileBounds = mGrid.tileBounds(coord);

        if (!(bounds().overlaps(tileBounds))) {
            terrainTile->setAllChildren(false);
        }
        else {
            if (bounds().overlaps(tileBounds.getSW())) {
                terrainTile->setChildSW();
            }
            if (bounds().overlaps(tileBounds.getNW())) {
                terrainTile->setChildNW();
            }
            if (bounds().overlaps(tileBounds.getNE())) {
                terrainTile->setChildNE();
            }
            if (bounds().overlaps(tileBounds.getSE())) {
                terrainTile->setChildSE();
            }
        }
    }

    return terrainTile;
}

2,写.terrain时，不能用gzwrite，那样写不出正确的.terrain，而是用fwrite方式.修正如下：在terraintile.cpp中改写

void Terrain::writeFile(const char *fileName) const
{
    FILE* terrainFile = fopen(fileName, "wb");

    if (terrainFile == NULL) {
        throw CTBException("Failed to open file");
    }

    // Write the height data,float应该是4 * TILE_CELL_SIZE
    //if (gzwrite(terrainFile, mHeights.data(), TILE_CELL_SIZE * 2) == 0)
    if (fwrite( mHeights.data(),  sizeof(i_terrain_height),TILE_CELL_SIZE, terrainFile) == 0)
    {
        fclose(terrainFile);
        throw CTBException("Failed to write height data");
    }

    // Write the child flags
    if (fputc(mChildren,terrainFile ) == -1) {
        fclose(terrainFile);
        throw CTBException("Failed to write child flags");
    }

    // Write the water mask
    if (fwrite(mMask,sizeof(char), mMaskLength,terrainFile ) == 0) {
        fclose(terrainFile);
        throw CTBException("Failed to write water mask");
    }

    fclose(terrainFile);
}
/*
void Terrain::writeFile(const char *fileName) const
{
    gzFile terrainFile = gzopen(fileName, "wb");

    if (terrainFile == NULL) {
        throw CTBException("Failed to open file");
    }

    // Write the height data,float应该是4 * TILE_CELL_SIZE
    //if (gzwrite(terrainFile, mHeights.data(), TILE_CELL_SIZE * 2) == 0)
    if (gzwrite(terrainFile, mHeights.data(), TILE_CELL_SIZE * sizeof(i_terrain_height)) == 0)
    {
        gzclose(terrainFile);
        throw CTBException("Failed to write height data");
    }

    // Write the child flags
    if (gzputc(terrainFile, mChildren) == -1) {
        gzclose(terrainFile);
        throw CTBException("Failed to write child flags");
    }

    // Write the water mask
    if (gzwrite(terrainFile, mMask, mMaskLength) == 0) {
        gzclose(terrainFile);
        throw CTBException("Failed to write water mask");
    }

    // Try and close the file
    switch (gzclose(terrainFile)) {
    case Z_OK:
        break;
    case Z_STREAM_ERROR:
    case Z_ERRNO:
    case Z_MEM_ERROR:
    case Z_BUF_ERROR:
    default:
        throw CTBException("Failed to close file");
    }
}

3,在环境变量中，如果设置GDAL_DATA，则没问题，如果在程序中按照CPLSetConfigOption("GDAL_DATA", strPath.c_str());的方式添加，则会出问题，原因出在

    OGRCoordinateTransformation *transformer = OGRCreateCoordinateTransformation(&srcSRS, &gridSRS);
      if (transformer == NULL) {
        throw CTBException("The source dataset to tile grid coordinate transformation could not be created");
      } else if (transformer->Transform(4, x, y) != true) {
        delete transformer;
        throw CTBException("Could not transform dataset bounds to tile spatial reference system");
      }

调试发现transformer = NULL;

因为gridSRS为空，OGRSpatialReference gridSRS = mGrid.getSRS();

mGrid构造函数中没有设置空间参考，，只要在构造函数中加上WGS84即可。

代码如下：

 Grid(i_tile tileSize,
       const CRSBounds extent,
       const OGRSpatialReference srs,
       unsigned short int rootTiles = 1,
       float zoomFactor = 2):
    mTileSize(tileSize),
    mExtent(extent),
    //mSRS(srs),
    mInitialResolution((extent.getWidth() / rootTiles) / tileSize ),
    mXOriginShift(extent.getWidth() / 2),
    mYOriginShift(extent.getHeight() / 2),
    mZoomFactor(zoomFactor)
  {
      mSRS.SetWellKnownGeogCS("WGS84");
  }