[原][osg][gdal]两种方式修改tiff高程

因为对于globalmap不熟悉，不怎么怎么修改高程，好像也没有这功能。

干脆自己手动修改了高程图tiff了

由于自身一直使用osg的 自己使用了osgDB直接读取tiff，修改后保存的。

同事小周一直研究gdal，她使用了gdal库直接改的，事实证明在专业gis处理上还是gdal更合适，现在把两种方式都总结一下：

第一种：通过osgDB修改tiff

#include <osg/Image>
#include <osgViewer/Viewer>
#include <osgEarth/Registry>
#include <osgEarth/ImageToHeightFieldConverter>
#include <osgEarth/ImageUtils>
#include <osgEarth/FileUtils>
#include <osgEarth/MapFrame>
#include <osgDB/FileUtils>
#include <osgDB/FileNameUtils>
#include <osgDB/ReadFile>
#include <osgDB/WriteFile>
#include <math.h>
using namespace osgEarth;


void floating(osg::HeightField* &hf)
{
    float floatingData = 179.0;
    float fBegin = 140.0;
    for (unsigned col = 0; col < hf->getNumColumns(); ++col)
    {
        for (unsigned row = 0; row < hf->getNumRows(); ++row)
        {
            float height = hf->getHeight(col, row);
            if (height < fBegin)
            {
                hf->setHeight(col, row, floatingData);
            }
        }
    }
}


void seaboard(osg::HeightField* &hf)
{

    float fBegin = 0.1;
    //float fEnd = 0.000001;
    float fLowestValue = 1000.0;
    int fWide = 100.0;
    int *fFlag = new int[hf->getNumColumns()*hf->getNumRows()];


    for (unsigned col = 0; col < hf->getNumColumns(); ++col)
    {
        for (unsigned row = 0; row < hf->getNumRows(); ++row)
        {
            fFlag[col*hf->getNumRows() + row] = 0;
            float height = hf->getHeight(col, row);
            if (height < fBegin)
            {
                fFlag[col*hf->getNumRows() + row] = 1;
                hf->setHeight(col, row, - fLowestValue);

                /*简单的计算
                float newh = -1000.0;
                if(height > 0.00001)
                newh = 0.1 - (0.1 - height)/ (0.1-0.00001)*1000.0;
                hf->setHeight(col, row, newh);*/
            }
        }
    }

    for (int i = 0; i < hf->getNumColumns()*hf->getNumRows(); i++)
    {
        if (fFlag[i] == 1)//如果这值在海面以下
        {
            bool isNearSide = false;
            int nowX = i / hf->getNumRows();
            int nowY = i%hf->getNumRows();
            for (int j = 0; j <= fWide; j++)
            {
                //从离此值最近的值开始找附近的岸边，往外延伸
                //向东南西北四个方向找，没层都遍历一圈
                for (int x = 0; x <= j; x++)
                {
                    //如果找到有岸边
                    int fDifValueX = x;
                    int fDifValueY = j - x;
                    int realX = nowX - fDifValueX;
                    if (realX > 0)
                    {
                        int realY = nowY - fDifValueY;
                        if (realY > 0)
                        {
                            if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边
                                isNearSide = true;
                        }
                        realY = nowY + fDifValueY;
                        if (realY < hf->getNumRows())
                        {
                            if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边
                                isNearSide = true;
                        }
                    }

                    realX = nowX + fDifValueX;
                    if (realX < hf->getNumColumns())
                    {
                        int realY = nowY - fDifValueY;
                        if (realY > 0)
                        {
                            if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边
                                isNearSide = true;
                        }
                        realY = nowY + fDifValueY;
                        if (realY < hf->getNumRows())
                        {
                            if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边
                                isNearSide = true;
                        }
                    }
                }

                //查找这个范围内是否有值,如果有值则用此值
                if (isNearSide)
                {
                    float fRealHeight = fBegin - j * fLowestValue / fWide;
                    hf->setHeight((i / hf->getNumRows()), (i % hf->getNumRows()), fRealHeight);
                    break;//退出当前寻找的延伸
                }
            }
        }
    }
    
    delete[]fFlag;
}


int main(int argc, char** argv)
{

    ReadResult result;
    osg::ref_ptr<osgDB::ReaderWriter> reader = osgDB::Registry::instance()->getReaderWriterForExtension("tif");
    std::string name("D:\srtm_19_032.tif");
    osgDB::ReaderWriter::Options* opt= NULL;
    osgDB::ReaderWriter::ReadResult rr = reader->readImage(name, opt);

    if (rr.validImage())
    {
        result = ReadResult(rr.takeImage());
        result.getImage()->setName("srtm_19_03.tif");
    }

    if (result.succeeded())
    {
        result.getObject();
        result.metadata();
        osg::ref_ptr<osg::Image> image = result.getImage();

        osgEarth::ImageToHeightFieldConverter conv;
        osg::HeightField* hf = conv.convert(image.get());

        //seaboard(hf);
        floating(hf);

        osg::Image* newimage = conv.convert(hf);
        std::string nameofnew("D:\srtm_19_03.tif");
        reader->writeImage(*newimage, nameofnew);
        
    }

    return 1;
}

第二种：通过gdal库修改tiff

getImgInfo.h

#include "gdal.h" 
#include "gdal_priv.h" 
#include <string>
using namespace std;

int getImgInfo(string szInFile, GDALDataset **poDataset, int *nbands, double **geoTrans, int *width, int *height, GDALDataType *gdt, const char** projRef, GDALRasterBand *** poBand)
{
    GDALAllRegister();

    GDALDataset *poDatasetTmp = *poDataset;
    poDatasetTmp = (GDALDataset*)GDALOpen(szInFile.c_str(), GA_ReadOnly);

    int widthTmp = *width, heightTmp = *height, nbandsTmp = *nbands;
    widthTmp = poDatasetTmp->GetRasterXSize();
    heightTmp = poDatasetTmp->GetRasterYSize();
    nbandsTmp = poDatasetTmp->GetRasterCount();

    GDALDataType gdtTmp = *gdt;
    gdtTmp = poDatasetTmp->GetRasterBand(1)->GetRasterDataType();  

    double *geoTransTmp = *geoTrans;
    geoTransTmp = new double[6];
    poDatasetTmp->GetGeoTransform(geoTransTmp);//获取地理坐标信息，地理坐标信息是一个含6个double型数据的数组，

    const char* projRefTmp = *projRef;
    projRefTmp = poDatasetTmp->GetProjectionRef();  //获取投影信息

    GDALRasterBand ** poBandTmp = *poBand;
    poBandTmp = new GDALRasterBand *[nbandsTmp];
    if (poBand == NULL)
    {
        cout << "GDALRasterBand ** poBand = new GDALRasterBand *[nBands]; failed!" << endl;
    }
    for (int i = 0; i < nbandsTmp; i++)
    {
        poBandTmp[i] = poDatasetTmp->GetRasterBand(i + 1);
    }

    *poDataset = poDatasetTmp;
    *nbands = nbandsTmp;
    *geoTrans = geoTransTmp;
    *width = widthTmp;
    *height = heightTmp;
    *gdt = gdtTmp;
    *projRef = projRefTmp;
    *poBand = poBandTmp;

    return 0;
}

main.cpp

#include "gdal.h" 
#include "gdal_priv.h" 
#include <iostream>
#include <string>

#include "getImgInfo.h"
using namespace std;

typedef unsigned short UINT;    //(0,65535)


int _tmain(int argc, _TCHAR* argv[])
{
    int ret=0;
    GDALAllRegister();

    string szInFile="F:\IMG_PROCESS\srtm_19_03.tif";

    double *adfGeoTransform;
    GDALDataType gdt;
    GDALDataset *poDataset;
    int Width, Height, nBands;
    GDALRasterBand ** poBand;
    const char* projRef;
    ret=getImgInfo(szInFile,&poDataset,&nBands,&adfGeoTransform,&Width,&Height ,&gdt,&projRef, &poBand);

    UINT **poBandBlock = new UINT*[nBands];
    if (poBandBlock==NULL)
    {
        cout<<"UINT **poBandBlock = new UINT *[outBand]; failed!"<<endl;
        system("pause");
        return -1;
    }
    UINT **poBandBlock_out = new UINT *[nBands];
    if (poBandBlock_out==NULL)
    {
        cout<<"UINT **poBandBlock_out = new byte *[outBand]; failed!"<<endl;
        system("pause");
        return -1;
    }

    int xHei=400;                //一次处理xHei=400行数据  分块处理
    int timePerxH=Height/xHei+1;

    GDALDriver *poDriver=NULL;
    poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

    string output_file="F:\IMG_PROCESS\srtm_19_03_tmp.tif";


    GDALDataset *BiomassDataset;
    BiomassDataset=poDriver->Create(output_file.c_str(),Width,Height,nBands,gdt, NULL); 

    int tmp=0;//防止无符号数越界

    //开始分块处理////
    for (int i=0;i<timePerxH - 1;i++)  //i表示分块处理需处理次数计数
    {
        for (int j=0;j<nBands;j++)   
        {
            poBandBlock[j] = new UINT[Width*xHei];
            if (poBandBlock[j]==NULL)
            {
                cout<<"poBandBlock[j] = new UINT[Width*xHei]; failed!"<<endl;
                system("pause");
                return -1;
            }
            poBandBlock_out[j] = new UINT[Width*xHei];
            if (poBandBlock_out[j]==NULL)
            {
                cout<<"poBandBlock_out[j] = new byte[Width*xHei]; failed!"<<endl;
                system("pause");
                return -1;
            }
        }
        for (int k=0; k<nBands; k++ )
        {
            poBand[k]->RasterIO(GF_Read, 0,i*xHei,Width,xHei, poBandBlock[k], Width,xHei, gdt, 0, 0);
        }

        //poBandBlock像素操作////
        for (int i=0;i<nBands;i++)
        {
            for (int j=0;j<xHei;j++)
            {
                for (int k=0;k<Width;k++)
                {
                    if (poBandBlock[i][j*Width+k]==0)
                    {
                        poBandBlock_out[i][j*Width+k]=0;
                    }
                    else if (poBandBlock[i][j*Width+k]>=1023)
                    {
                        poBandBlock_out[i][j*Width+k]=1023;
                    }
                    else
                    {
                        poBandBlock_out[i][j*Width+k] = poBandBlock[i][j*Width+k];            
                    }
                }
            }
        }
        for(int k = 0; k < nBands; k++)
        {
            BiomassDataset->GetRasterBand (k+1)->RasterIO(GF_Write,0,xHei*i,Width,xHei,poBandBlock_out[k],Width,xHei,gdt,0,0);
        }

        for (int j=0;j<nBands;j++)   
        {
            delete []poBandBlock[j];    
            poBandBlock[j]=NULL;

            delete poBandBlock_out[j];    
            poBandBlock_out[j]=NULL;            
        }
    }

    //////剩余行数处理////
    int foreH=(timePerxH-1)*xHei;
    int leftH=Height-foreH;
    for (int j=0;j<nBands;j++)   
    {
        poBandBlock[j] = new UINT[Width*leftH];
        if (poBandBlock[j]==NULL)
        {
            cout<<"poBandBlock[j] = new UINT[Width*leftH]; failed!"<<endl;
            system("pause");
            return -1;
        }
        poBandBlock_out[j] = new UINT[Width*leftH];
        if (poBandBlock_out[j]==NULL)
        {
            cout<<"poBandBlock_out[j] = new byte[Width*leftH]; failed!"<<endl;
            system("pause");
            return -1;
        }
    }
    for (int k=0; k<nBands; k++ )
    {
        poBand[k]->RasterIO(GF_Read, 0,foreH,Width,leftH, poBandBlock[k], Width,leftH, gdt, 0, 0);
    }
    //poBandBlock像素操作////
    for (int i=0;i<nBands;i++)
    {
        for (int j=0;j<leftH;j++)
        {
            for (int k=0;k<Width;k++)
            {
                if (poBandBlock[i][j*Width+k]==0)
                {
                    poBandBlock_out[i][j*Width+k]=0;
                }
                else if (poBandBlock[i][j*Width+k]>=1023)
                {
                    poBandBlock_out[i][j*Width+k]=1023;
                }
                else
                {
                    poBandBlock_out[i][j*Width+k]=poBandBlock[i][j*Width+k];            
                }
            }
        }
    }
    for(int k = 0; k < nBands; k++)
    {
        BiomassDataset->GetRasterBand(k+1)->RasterIO(GF_Write,0,foreH,Width,leftH,poBandBlock_out[k],Width,leftH,gdt,0,0);
    }
    for (int j=0;j<nBands;j++)   
    {
        delete poBandBlock[j];    
        poBandBlock[j]=NULL;

        delete poBandBlock_out[j];    
        poBandBlock_out[j]=NULL;            
    }
    

    ////////对指定的某个区域数据重新赋值//////////////////////////////////////////////////////////////////
    for (int j = 0; j < nBands; j++)
    {
        poBandBlock[j] = new UINT[1207 * 1215];
        if (poBandBlock[j] == NULL)
        {
            cout << "poBandBlock[j] = new UINT[Width*leftH]; failed!" << endl;
            system("pause");
            return -1;
        }
        poBandBlock_out[j] = new UINT[1207 * 1215];
        if (poBandBlock_out[j] == NULL)
        {
            cout << "poBandBlock_out[j] = new byte[Width*leftH]; failed!" << endl;
            system("pause");
            return -1;
        }
    }
    for (int k = 0; k < nBands; k++)
    {
        poBand[k]->RasterIO(GF_Read, 3594, 2386, 1207, 1215, poBandBlock[k], 1207, 1215, gdt, 0, 0);
    }
    //poBandBlock像素操作//
    for (int i = 0; i < nBands; i++)
    {
        for (int j = 0; j < 1215; j++)
        {
            for (int k = 0; k < 1207; k++)
            {    
                poBandBlock_out[i][j*1207 + k] = 179;
            }
        }
    }
    for (int k = 0; k < nBands; k++)
    {
        BiomassDataset->GetRasterBand(k + 1)->RasterIO(GF_Write, 3594, 2386, 1207, 1215, poBandBlock_out[k], 1207, 1215, gdt, 0, 0);
    }
    for (int j = 0; j < nBands; j++)
    {
        delete poBandBlock[j];
        poBandBlock[j] = NULL;

        delete poBandBlock_out[j];
        poBandBlock_out[j] = NULL;
    }
    //////////////////////////////////////////////////////////////////////////


    delete poBandBlock;
    poBandBlock=NULL;
    delete poBandBlock_out;
    poBandBlock_out=NULL;


    BiomassDataset->SetGeoTransform(adfGeoTransform);
    BiomassDataset->SetProjection(projRef);
    system("pause");

}