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

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

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

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

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

第一种：通过osgDB修改tiff

  1 #include <osg/Image>
  2 #include <osgViewer/Viewer>
  3 #include <osgEarth/Registry>
  4 #include <osgEarth/ImageToHeightFieldConverter>
  5 #include <osgEarth/ImageUtils>
  6 #include <osgEarth/FileUtils>
  7 #include <osgEarth/MapFrame>
  8 #include <osgDB/FileUtils>
  9 #include <osgDB/FileNameUtils>
 10 #include <osgDB/ReadFile>
 11 #include <osgDB/WriteFile>
 12 #include <math.h>
 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<osgDB::ReaderWriter> 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<osg::Image> 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 <string>
 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 <iostream>
  4 #include <string>
  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!"<<endl;
 31         system("pause");
 32         return -1;
 33     }
 34     UINT **poBandBlock_out = new UINT *[nBands];
 35     if (poBandBlock_out==NULL)
 36     {
 37         cout<<"UINT **poBandBlock_out = new byte *[outBand]; failed!"<<endl;
 38         system("pause");
 39         return -1;
 40     }
 41 
 42     int xHei=400;                //一次处理xHei=400行数据  分块处理
 43     int timePerxH=Height/xHei+1;
 44 
 45     GDALDriver *poDriver=NULL;
 46     poDriver = GetGDALDriverManager()->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<timePerxH - 1;i++)  //i表示分块处理需处理次数计数
 58     {
 59         for (int j=0;j<nBands;j++)   
 60         {
 61             poBandBlock[j] = new UINT[Width*xHei];
 62             if (poBandBlock[j]==NULL)
 63             {
 64                 cout<<"poBandBlock[j] = new UINT[Width*xHei]; failed!"<<endl;
 65                 system("pause");
 66                 return -1;
 67             }
 68             poBandBlock_out[j] = new UINT[Width*xHei];
 69             if (poBandBlock_out[j]==NULL)
 70             {
 71                 cout<<"poBandBlock_out[j] = new byte[Width*xHei]; failed!"<<endl;
 72                 system("pause");
 73                 return -1;
 74             }
 75         }
 76         for (int k=0; k<nBands; k++ )
 77         {
 78             poBand[k]->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<nBands;i++)
 83         {
 84             for (int j=0;j<xHei;j++)
 85             {
 86                 for (int k=0;k<Width;k++)
 87                 {
 88                     if (poBandBlock[i][j*Width+k]==0)
 89                     {
 90                         poBandBlock_out[i][j*Width+k]=0;
 91                     }
 92                     else if (poBandBlock[i][j*Width+k]>=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<nBands;j++)   
109         {
110             delete []poBandBlock[j];    
111             poBandBlock[j]=NULL;
112 
113             delete poBandBlock_out[j];    
114             poBandBlock_out[j]=NULL;            
115         }
116     }
117 
118     //////剩余行数处理////
119     int foreH=(timePerxH-1)*xHei;
120     int leftH=Height-foreH;
121     for (int j=0;j<nBands;j++)   
122     {
123         poBandBlock[j] = new UINT[Width*leftH];
124         if (poBandBlock[j]==NULL)
125         {
126             cout<<"poBandBlock[j] = new UINT[Width*leftH]; failed!"<<endl;
127             system("pause");
128             return -1;
129         }
130         poBandBlock_out[j] = new UINT[Width*leftH];
131         if (poBandBlock_out[j]==NULL)
132         {
133             cout<<"poBandBlock_out[j] = new byte[Width*leftH]; failed!"<<endl;
134             system("pause");
135             return -1;
136         }
137     }
138     for (int k=0; k<nBands; k++ )
139     {
140         poBand[k]->RasterIO(GF_Read, 0,foreH,Width,leftH, poBandBlock[k], Width,leftH, gdt, 0, 0);
141     }
142     //poBandBlock像素操作////
143     for (int i=0;i<nBands;i++)
144     {
145         for (int j=0;j<leftH;j++)
146         {
147             for (int k=0;k<Width;k++)
148             {
149                 if (poBandBlock[i][j*Width+k]==0)
150                 {
151                     poBandBlock_out[i][j*Width+k]=0;
152                 }
153                 else if (poBandBlock[i][j*Width+k]>=1023)
154                 {
155                     poBandBlock_out[i][j*Width+k]=1023;
156                 }
157                 else
158                 {
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<nBands;j++)   
169     {
170         delete poBandBlock[j];    
171         poBandBlock[j]=NULL;
172 
173         delete poBandBlock_out[j];    
174         poBandBlock_out[j]=NULL;            
175     }
176     
177 
178     ////////对指定的某个区域数据重新赋值//////////////////////////////////////////////////////////////////
179     for (int j = 0; j < nBands; j++)
180     {
181         poBandBlock[j] = new UINT[1207 * 1215];
182         if (poBandBlock[j] == NULL)
183         {
184             cout << "poBandBlock[j] = new UINT[Width*leftH]; failed!" << endl;
185             system("pause");
186             return -1;
187         }
188         poBandBlock_out[j] = new UINT[1207 * 1215];
189         if (poBandBlock_out[j] == NULL)
190         {
191             cout << "poBandBlock_out[j] = new byte[Width*leftH]; failed!" << endl;
192             system("pause");
193             return -1;
194         }
195     }
196     for (int k = 0; k < nBands; k++)
197     {
198         poBand[k]->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 }

相关阅读:
HTML基础
 JavaScript基础
 J2EE技术（六）——JSP
2011.2—至今年度总结
 J2EE技术（二）——JDBC
DRP项目总结（一）——项目简介
 CSS基础
 J2EE技术（五）——EJB
J2EE技术（三）——JMS
J2EE技术（四）——JavaMail
原文地址：https://www.cnblogs.com/lyggqm/p/6737862.html