//总体来说这个过程就是构建数据源->构建层->构建要素->构建形状->关闭数据源。 //要包含的GDAL头文件 #include <gdal_priv.h> #include <ogrsf_frmts.h> #include <iostream> using namespace std; #pragma comment(lib,"gdal_i.lib") bool Creatshape(const char* pszFileName ,int line,int row); #include <tchar.h> //_TCHAR* 类型在该头文件里 int _tmain(int argc, _TCHAR* argv[]) { const char *pszFileName="C:\Users\Public\Pictures\Sample Pictures\srtm_51_03.tif"; Creatshape(pszFileName,7,9); system("pause"); return 0; } /************************************************************************/ /*创建过程:构建数据源->构建层->构建要素->构建形状->关闭数据源 /* 参数pszFileName 为输入的文件名 参数linenum为划分的行数 参数rows为划分的列数*/ /************************************************************************/ bool Creatshape(const char* pszFileName ,int linenum,int rows) { //获取影像信息 GDALDataset *poDataSet; GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO"); poDataSet=(GDALDataset*)GDALOpen(pszFileName,GA_ReadOnly); //打开数据集 if (poDataSet==NULL) { "Failed to open this dataset!"; //代开失败的话则给出提示并退出 exit(1); } double Trans[6];//坐标转换参数数组,作为GetGeoTransform()函数的参数 int width,height; //影像的像行列数; width=poDataSet->GetRasterXSize(); //获取影像列数,为后续划分网格做准备 height=poDataSet->GetRasterYSize();//获取影像行数,为后续划分网格做准备 poDataSet->GetGeoTransform(Trans); for (int i=0;i<6;i++) { cout<<Trans[i]<<endl; //循环输出仿射变换参数; } //注册shape文件驱动 const char* pszDriverName="ESRI Shapefile"; OGRSFDriver *poDriver; OGRRegisterAll(); poDriver=OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverName); if (poDriver==NULL) { printf("%s driver is not available!",pszDriverName); exit(1); } //创建shape文件; OGRDataSource *poDS; //创建一个叫myshapefile的目录,存放生成的文件; //如果名字有.shp后缀,则直接在当前目录下生成文件; poDS=poDriver->CreateDataSource("myshapefile.shp",NULL); if (poDS==NULL) { printf("Create my shape file failed!"); exit(1); } //创建输出图层; OGRLayer *poLayer;; const char *prj=poDataSet->GetProjectionRef(); //获取栅格影像的空间参考信息 cout<<"栅格数据空间参考信息为: "<<prj<<endl<<endl; OGRSpatialReference oSRS; oSRS.SetFromUserInput(prj); //将获得的空间参考信息字符串做为文本一次性赋给矢量数据的OGRSpatialReference对象; poLayer=poDS->CreateLayer(pszFileName,&oSRS, wkbUnknown, NULL); if (poLayer==NULL) { printf("Creat layer failed!"); exit(1); } //创建层数据的属性fields; OGRFieldDefn oField("Point",OFTString); oField.SetWidth(10); if (poLayer->CreateField(&oField)!=OGRERR_NONE) { printf("Create Point Field Failed!"); exit(1); } //创建features,写入feature到磁盘; OGRFeature *poFeature; poFeature=OGRFeature::CreateFeature(poLayer->GetLayerDefn()); //绘制外边框 OGRLineString Line; OGRPoint Point1(Trans[0],Trans[3]); OGRPoint Point2(Trans[0]+width*Trans[1],Trans[3]); OGRPoint Point3(Trans[0]+width*Trans[1],Trans[3]+width*Trans[4]+height*Trans[5]); OGRPoint Point4(Trans[0],Trans[3]+width*Trans[4]+height*Trans[5]);//四个角点的地理坐标,通过行列号计算地理坐标,是通过六参数得到的; Line.addPoint(&Point1); Line.addPoint(&Point2); Line.addPoint(&Point3); Line.addPoint(&Point4); Line.addPoint(&Point1); //水平方向加线; OGRLineString SubHline[50]; OGRPoint PointLeft[50],PointRight[50]; for (int i=1;i<linenum;i++) { PointLeft[i].setX(Trans[0]); //设置左边框上要加线的起点X坐标 PointLeft[i].setY((Point4.getY()-Trans[3])/linenum*i+Trans[3]);//设置左边框上要加线的起点Y坐标 PointRight[i].setX(Point2.getX());//设置右边框上要加线的起点X坐标 PointRight[i].setY((Point4.getY()-Trans[3])/linenum*i+Trans[3]);//设置右边框上要加线的起点Y坐标 } for (int i=1;i<linenum;i++) { SubHline[i].addPoint(&PointLeft[i]); //左边框上加点; SubHline[i].addPoint(&PointRight[i]);//右边框上加点; if (i<linenum-1) { SubHline[i].addPoint(&PointRight[i+1]); //从右边框的上一点转到下一点,避免交叉斜线的出现; } Line.addSubLineString(&SubHline[i]);//将SubHline数组中的每一个线做为子线段添加到Line对象中; } Line.addPoint(&Point2); //垂直方向加线 OGRLineString SubVline[50]; OGRPoint PointUp[50],PointDown[50]; for (int j=1;j<rows;j++) { //添加上边框上的点; PointUp[j].setX((Point2.getX()-Trans[0])/rows*j+Trans[0]); PointUp[j].setY(Trans[3]); //添加下边框上的点; PointDown[j].setX((Point2.getX()-Trans[0])/rows*j+Trans[0]); PointDown[j].setY(Point4.getY()); } for (int j=1;j<rows;j++) { SubVline[j].addPoint(&PointUp[j]); SubVline[j].addPoint(&PointDown[j]); if (j<rows-1) { SubVline[j].addPoint(&PointDown[j+1]); //从下边框的前一点转到后一点,避免交叉斜线的出现; } Line.addSubLineString(&SubVline[j]); } //poFeature->SetGeometryDirectly(&Line);//问题出在这啊!!用SetGeometryDirectly报错,要用SetGeometry poFeature->SetGeometry(&Line); if (poLayer->CreateFeature(poFeature)!=OGRERR_NONE) { printf("Failed create feature in shapefile!"); exit(1); } OGRFeature::DestroyFeature(poFeature); OGRDataSource::DestroyDataSource(poDS); printf("创建矢量数据成功! "); system("pause"); return TRUE; }
用ENVI打开生成的shp文件:
注意:在计算四个角点的坐标时,用到了六参数,参考博客:http://blog.csdn.net/ivanljf/article/details/9226463