ORG/GDAL利用csv文件生成shapefile点文件

       输入csv文件,确保前三列为ID/经纬度,后面的列数及列的类型不限制,生成对应的点文件,默认设置地理坐标系为WGS84,投影坐标为通用墨卡托49N。

       

#include <QtCore>
#include "ogrsf_frmts.h"
#include "gdal.h"
#include "gdal_priv.h"
#include "cpl_string.h" 
#include <string>
#include <iostream>
#include <strstream>
#include <time.h>
using namespace std;

struct DataStructure
{
	int nFID;
	double dlng;
	double dlat;
	QStringList vattr;
};
QStringList vAttr;

QList<DataStructure> getCsvTableData(const char* strFileName)
{
	QFile infile(strFileName);
	if (!infile.open(QIODevice::ReadOnly))
	{
		cout << strFileName << " File open Error!!!\n ";
	}
	QTextStream _in(&infile);
	//读取表头。在此之前必须把前三行设置为FID, lon, lat;
	
	vAttr.clear();

	QString smsg = _in.readLine();
	QStringList slist = smsg.split(",");
	//提取附属属性;
	if (slist.size() > 2)
	{
		for(int i = 2; i < slist.size(); i ++)
		{
			vAttr.append(slist[i]);
		}
	}

	QList <DataStructure> vDatas;
	vDatas.clear();
	while(! _in.atEnd())
	{
		smsg = _in.readLine();
		slist = smsg.split(",");
		DataStructure mydata;
		//mydata.nFID = slist[0].trimmed().toInt();
		mydata.dlng = slist[1].trimmed().toDouble();//经度
		mydata.dlat = slist[0].trimmed().toDouble();//纬度
		for (int i = 0; i < vAttr.size(); i ++)
		{
			mydata.vattr.append(slist[i + 2]);
		}
		vDatas.append(mydata);
	}
	return vDatas;

}

void convertToShp(QList<DataStructure> vData, char *outshp)
{   
	//使属性表字段支持中文
	CPLSetConfigOption("SHAPE_ENCODING","");
	OGRRegisterAll();//注册所有的驱动
	//创建ESRI shp文件
	char *pszDriverName = "ESRI Shapefile";
	//调用对Shape文件读写的Driver
	OGRSFDriver *poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverName);
	if (poDriver == NULL)
	{
		cout<<pszDriverName<<"驱动不可用!"<<endl;
		return;
	}
	//创建数据源
	OGRDataSource *poDs = poDriver->CreateDataSource(outshp, NULL);
	if (poDs == NULL)
	{
		cout<<"DataSource Creation Error"<<endl;
		return;
	}
	//创建图层Layer
	string outShapName = outshp;
	string layerName = outShapName.substr(0, outShapName.length()-4);
	//layerName.c_str()表示将string转为char *类型
	//参数说明:新图层名称,坐标系,图层的几何类型,创建选项,与驱动有关
	OGRSpatialReference poSpatialRef;
	poSpatialRef.SetProjCS( "UTM 49(WGS84) in northern hemisphere." );
	poSpatialRef.SetWellKnownGeogCS("WGS84");
	poSpatialRef.SetUTM( 49, TRUE );
	
	OGRLayer *poLayer = poDs->CreateLayer(layerName.c_str(), &poSpatialRef, wkbPoint, NULL);
	if (poLayer == NULL)
	{
		cout<<"Layer Creation Failed"<<endl;
		OGRDataSource::DestroyDataSource(poDs);
		return;
	}
	//下面创建属性表
	//先创建一个“ID”整型属性
	OGRFieldDefn oFieldId("nFID", OFTInteger);
	oFieldId.SetWidth(5);
	poLayer->CreateField(&oFieldId);
	//再创建一个"X"double属性
	OGRFieldDefn oFieldX("lon", OFTString);
	oFieldX.SetWidth(30);
	poLayer->CreateField(&oFieldX);
	//再创建一个"Y"double属性
	OGRFieldDefn oFieldY("lat", OFTString);
	oFieldY.SetWidth(30);
	poLayer->CreateField(&oFieldY);
	//创建其他属性
	for (int i =0; i < vAttr.size(); i ++)
	{
		OGRFieldDefn oField(vAttr[i].toStdString().c_str(), OFTString);
		oField.SetWidth(20);
		poLayer->CreateField(&oField);

	}
	
	//创建一个feature
	for (int i =0; i < vData.size(); i ++)
	{
		OGRFeature *poFeature;  
		poFeature = OGRFeature::CreateFeature(poLayer->GetLayerDefn());//GetLayerDefn()获取当前图层的属性表结构
		//给属性表中我们刚创建的列赋值;

		poFeature->SetField("nFID", i);
		poFeature->SetField("lon", vData[i].dlng);
		poFeature->SetField("lat", vData[i].dlat);
		for (int j = 0; j < vAttr.size(); j ++)
		{
			poFeature->SetField(vAttr[j].toStdString().c_str(), vData[i].vattr[j].toStdString().c_str());
		}
		//创建一个OGRPoint对象
		OGRPoint point;
		point.setX(vData[i].dlng);
		point.setY(vData[i].dlat);
		//point.setZ(0);

		poFeature->SetGeometry(&point);

		if(poLayer->CreateFeature(poFeature) != OGRERR_NONE )
		{
			printf( "Failed to create feature in shapefile.\n" );
			exit( 1 );
		}
		OGRFeature::DestroyFeature(poFeature);
	OGRDataSource::DestroyDataSource(poDs);     
}



int main(int argc, char *argv[])
{
	QTextCodec* codec =QTextCodec::codecForLocale();
	QTextCodec::setCodecForCStrings(codec);
	QTextCodec::setCodecForTr(codec);

	time_t t1 = time(NULL);

	QList<DataStructure> vData = getCsvTableData("../2016_POIs.csv");
	convertToShp(vData, "../poi_pro.shp");
	time_t t2 = time(NULL);
	cout << t2-t1<<endl;
	return 0;
}

       以87万个POI为例,在i3, 4G内存的笔记本电脑上运行,用时44s,相比于ArcGIS的导出shp文件还是很快的。


                 

       如上属性表,csv表内的其余属性也会自动加入,只要宽度不超过20.

转载自:https://blog.csdn.net/liuph_/article/details/52262405

You may also like...

退出移动版