• C++、GDAL创建shapefile,并向矢量文件中添加网格


    //总体来说这个过程就是构建数据源->构建层->构建要素->构建形状->关闭数据源。
    //要包含的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

  • 相关阅读:
    菜鸟学存储:网络存储IP SAN与IB SAN
    读xml高手
    预先加载图片
    xred520
    最简单准确的硬盘整数分区设置操作方法
    Google 每天处理约 20000TB 的数据
    IE 8 无法正常使用网站后台编辑器问题
    常用的JS技术1
    adodb stream 使用说明
    [Tools] JDGUI(Java Decompiler)
  • 原文地址:https://www.cnblogs.com/aukle/p/3225769.html
Copyright © 2020-2023  润新知