[原][osg][gdal]两种方式修改tiff高程-白红宇

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

阅读量：6297 次

发布时间：2019-06-22

本文共 8757 字，大约阅读时间需要 29 分钟。

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

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

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

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

第一种：通过osgDB修改tiff

1 #include 
     
        2 #include 
      
         3 #include 
       
          4 #include 
        
           5 #include 
         
            6 #include 
          
            7 #include 
           
             8 #include 
            
              9 #include 
             
               10 #include 
              
                11 #include 
               
                 12 #include 
                
                  13 using namespace osgEarth; 14 15 16 void floating(osg::HeightField* &hf) 17 { 18 float floatingData = 179.0; 19 float fBegin = 140.0; 20 for (unsigned col = 0; col < hf->getNumColumns(); ++col) 21 { 22 for (unsigned row = 0; row < hf->getNumRows(); ++row) 23 { 24 float height = hf->getHeight(col, row); 25 if (height < fBegin) 26 { 27 hf->setHeight(col, row, floatingData); 28 } 29 } 30 } 31 } 32 33 34 void seaboard(osg::HeightField* &hf) 35 { 36 37 float fBegin = 0.1; 38 //float fEnd = 0.000001; 39 float fLowestValue = 1000.0; 40 int fWide = 100.0; 41 int *fFlag = new int[hf->getNumColumns()*hf->getNumRows()]; 42 43 44 for (unsigned col = 0; col < hf->getNumColumns(); ++col) 45 { 46 for (unsigned row = 0; row < hf->getNumRows(); ++row) 47 { 48 fFlag[col*hf->getNumRows() + row] = 0; 49 float height = hf->getHeight(col, row); 50 if (height < fBegin) 51 { 52 fFlag[col*hf->getNumRows() + row] = 1; 53 hf->setHeight(col, row, - fLowestValue); 54 55 /*简单的计算 56 float newh = -1000.0; 57 if(height > 0.00001) 58 newh = 0.1 - (0.1 - height)/ (0.1-0.00001)*1000.0; 59 hf->setHeight(col, row, newh);*/ 60 } 61 } 62 } 63 64 for (int i = 0; i < hf->getNumColumns()*hf->getNumRows(); i++) 65 { 66 if (fFlag[i] == 1)//如果这值在海面以下 67 { 68 bool isNearSide = false; 69 int nowX = i / hf->getNumRows(); 70 int nowY = i%hf->getNumRows(); 71 for (int j = 0; j <= fWide; j++) 72 { 73 //从离此值最近的值开始找附近的岸边，往外延伸 74 //向东南西北四个方向找，没层都遍历一圈 75 for (int x = 0; x <= j; x++) 76 { 77 //如果找到有岸边 78 int fDifValueX = x; 79 int fDifValueY = j - x; 80 int realX = nowX - fDifValueX; 81 if (realX > 0) 82 { 83 int realY = nowY - fDifValueY; 84 if (realY > 0) 85 { 86 if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边 87 isNearSide = true; 88 } 89 realY = nowY + fDifValueY; 90 if (realY < hf->getNumRows()) 91 { 92 if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边 93 isNearSide = true; 94 } 95 } 96 97 realX = nowX + fDifValueX; 98 if (realX < hf->getNumColumns()) 99 {100 int realY = nowY - fDifValueY;101 if (realY > 0)102 {103 if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边104 isNearSide = true;105 }106 realY = nowY + fDifValueY;107 if (realY < hf->getNumRows())108 {109 if (fFlag[realX*hf->getNumRows() + realY] == 0)//如果是岸边110 isNearSide = true;111 }112 }113 }114 115 //查找这个范围内是否有值,如果有值则用此值116 if (isNearSide)117 {118 float fRealHeight = fBegin - j * fLowestValue / fWide;119 hf->setHeight((i / hf->getNumRows()), (i % hf->getNumRows()), fRealHeight);120 break;//退出当前寻找的延伸121 }122 }123 }124 }125 126 delete[]fFlag;127 }128 129 130 int main(int argc, char** argv)131 {132 133 ReadResult result;134 osg::ref_ptr
                 
                   reader = osgDB::Registry::instance()->getReaderWriterForExtension("tif");135 std::string name("D:\\srtm_19_032.tif");136 osgDB::ReaderWriter::Options* opt= NULL;137 osgDB::ReaderWriter::ReadResult rr = reader->readImage(name, opt);138 139 if (rr.validImage())140 {141 result = ReadResult(rr.takeImage());142 result.getImage()->setName("srtm_19_03.tif");143 }144 145 if (result.succeeded())146 {147 result.getObject();148 result.metadata();149 osg::ref_ptr
                  
                    image = result.getImage();150 151 osgEarth::ImageToHeightFieldConverter conv;152 osg::HeightField* hf = conv.convert(image.get());153 154 //seaboard(hf);155 floating(hf);156 157 osg::Image* newimage = conv.convert(hf);158 std::string nameofnew("D:\\srtm_19_03.tif");159 reader->writeImage(*newimage, nameofnew);160 161 }162 163 return 1;164 }

第二种：通过gdal库修改tiff

getImgInfo.h

1 #include "gdal.h"  2 #include "gdal_priv.h"  3 #include 
     
       4 using namespace std; 5  6 int getImgInfo(string szInFile, GDALDataset **poDataset, int *nbands, double **geoTrans, int *width, int *height, GDALDataType *gdt, const char** projRef, GDALRasterBand *** poBand) 7 { 8     GDALAllRegister(); 9 10     GDALDataset *poDatasetTmp = *poDataset;11     poDatasetTmp = (GDALDataset*)GDALOpen(szInFile.c_str(), GA_ReadOnly);12 13     int widthTmp = *width, heightTmp = *height, nbandsTmp = *nbands;14     widthTmp = poDatasetTmp->GetRasterXSize();15     heightTmp = poDatasetTmp->GetRasterYSize();16     nbandsTmp = poDatasetTmp->GetRasterCount();17 18     GDALDataType gdtTmp = *gdt;19     gdtTmp = poDatasetTmp->GetRasterBand(1)->GetRasterDataType();  20 21     double *geoTransTmp = *geoTrans;22     geoTransTmp = new double[6];23     poDatasetTmp->GetGeoTransform(geoTransTmp);//获取地理坐标信息，地理坐标信息是一个含6个double型数据的数组，24 25     const char* projRefTmp = *projRef;26     projRefTmp = poDatasetTmp->GetProjectionRef();  //获取投影信息27 28     GDALRasterBand ** poBandTmp = *poBand;29     poBandTmp = new GDALRasterBand *[nbandsTmp];30     if (poBand == NULL)31     {32         cout << "GDALRasterBand ** poBand = new GDALRasterBand *[nBands]; failed!" << endl;33     }34     for (int i = 0; i < nbandsTmp; i++)35     {36         poBandTmp[i] = poDatasetTmp->GetRasterBand(i + 1);37     }38 39     *poDataset = poDatasetTmp;40     *nbands = nbandsTmp;41     *geoTrans = geoTransTmp;42     *width = widthTmp;43     *height = heightTmp;44     *gdt = gdtTmp;45     *projRef = projRefTmp;46     *poBand = poBandTmp;47 48     return 0;49 }

main.cpp

1 #include "gdal.h"   2 #include "gdal_priv.h"   3 #include 
     
        4 #include 
      
         5   6 #include "getImgInfo.h"  7 using namespace std;  8   9 typedef unsigned short UINT;    //(0,65535) 10  11  12 int _tmain(int argc, _TCHAR* argv[]) 13 { 14     int ret=0; 15     GDALAllRegister(); 16  17     string szInFile="F:\\IMG_PROCESS\\srtm_19_03.tif"; 18  19     double *adfGeoTransform; 20     GDALDataType gdt; 21     GDALDataset *poDataset; 22     int Width, Height, nBands; 23     GDALRasterBand ** poBand; 24     const char* projRef; 25     ret=getImgInfo(szInFile,&poDataset,&nBands,&adfGeoTransform,&Width,&Height ,&gdt,&projRef, &poBand); 26  27     UINT **poBandBlock = new UINT*[nBands]; 28     if (poBandBlock==NULL) 29     { 30         cout<<"UINT **poBandBlock = new UINT *[outBand]; failed!"<
       
        GetDriverByName("GTiff"); 47  48     string output_file="F:\\IMG_PROCESS\\srtm_19_03_tmp.tif"; 49  50  51     GDALDataset *BiomassDataset; 52     BiomassDataset=poDriver->Create(output_file.c_str(),Width,Height,nBands,gdt, NULL);  53  54     int tmp=0;//防止无符号数越界 55  56     //开始分块处理 57     for (int i=0;i
        
         RasterIO(GF_Read, 0,i*xHei,Width,xHei, poBandBlock[k], Width,xHei, gdt, 0, 0); 79         } 80  81         //poBandBlock像素操作 82         for (int i=0;i
         
          =1023) 93                     { 94                         poBandBlock_out[i][j*Width+k]=1023; 95                     } 96                     else 97                     { 98                         poBandBlock_out[i][j*Width+k] = poBandBlock[i][j*Width+k];             99                     }100                 }101             }102         }103         for(int k = 0; k < nBands; k++)104         {105             BiomassDataset->GetRasterBand (k+1)->RasterIO(GF_Write,0,xHei*i,Width,xHei,poBandBlock_out[k],Width,xHei,gdt,0,0);106         }107 108         for (int j=0;j
          
           RasterIO(GF_Read, 0,foreH,Width,leftH, poBandBlock[k], Width,leftH, gdt, 0, 0);141 }142 //poBandBlock像素操作143 for (int i=0;i
           
            =1023)154 {155 poBandBlock_out[i][j*Width+k]=1023;156 }157 else158 {159 poBandBlock_out[i][j*Width+k]=poBandBlock[i][j*Width+k]; 160 }161 }162 }163 }164 for(int k = 0; k < nBands; k++)165 {166 BiomassDataset->GetRasterBand(k+1)->RasterIO(GF_Write,0,foreH,Width,leftH,poBandBlock_out[k],Width,leftH,gdt,0,0);167 }168 for (int j=0;j
            
             RasterIO(GF_Read, 3594, 2386, 1207, 1215, poBandBlock[k], 1207, 1215, gdt, 0, 0);199 }200 //poBandBlock像素操作//201 for (int i = 0; i < nBands; i++)202 {203 for (int j = 0; j < 1215; j++)204 {205 for (int k = 0; k < 1207; k++)206 { 207 poBandBlock_out[i][j*1207 + k] = 179;208 }209 }210 }211 for (int k = 0; k < nBands; k++)212 {213 BiomassDataset->GetRasterBand(k + 1)->RasterIO(GF_Write, 3594, 2386, 1207, 1215, poBandBlock_out[k], 1207, 1215, gdt, 0, 0);214 }215 for (int j = 0; j < nBands; j++)216 {217 delete poBandBlock[j];218 poBandBlock[j] = NULL;219 220 delete poBandBlock_out[j];221 poBandBlock_out[j] = NULL;222 }223 //224 225 226 delete poBandBlock;227 poBandBlock=NULL;228 delete poBandBlock_out;229 poBandBlock_out=NULL;230 231 232 BiomassDataset->SetGeoTransform(adfGeoTransform);233 BiomassDataset->SetProjection(projRef);234 system("pause");235 236 }