摘要：
因为我不熟悉全球地图，所以我不怎么修改高程，而且我似乎没有这个功能。我只是手动修改了高程图tiff。由于我一直使用osg，所以我使用osgDB直接读取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 }