本文介绍如何使用GDAL/OGR 库对shapefile文件进行简单的操作,包括读取、创建、修改等.在GDAL官网上有读写shp文件的介绍,主要是针对点要素的操作例子:点击打开链接
1.读取shp文件
void ExtracInfo::ReadShapeFile()
{
OGRRegisterAll();
OGRDataSource *poDS;
poDS = OGRSFDriverRegistrar::Open("F:/test.shp", FALSE);
if (poDS == NULL)
{
// printf("Open failed \n");
}
OGRLayer *poLayer;
poLayer = poDS->GetLayerByName("test");
int layercount = poDS->GetLayerCount();
// printf("layer count is : %d\n", layercount);
OGRFeature *poFeature;
poLayer->ResetReading();
while((poFeature = poLayer->GetNextFeature()) != NULL)
{
//printf("hello!\n");
//get field
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
int iField;
int fieldcount = poFDefn->GetFieldCount();
// printf("Field count is : %d\n", poFDefn->GetFieldCount());
int fieldv1;
double fieldv2;
string feildv3;
for (iField = 0; iField < poFDefn->GetFieldCount(); iField ++)
{
OGRFieldDefn * poFieldDefn = poFDefn->GetFieldDefn(iField);
if (poFieldDefn->GetType() == OFTInteger)
{
fieldv1 = poFeature->GetFieldAsInteger(iField);
// printf("%d ", poFeature->GetFieldAsInteger(iField));
}
else if (poFieldDefn->GetType() == OFTReal)
{
fieldv2 = poFeature->GetFieldAsDouble(iField);
// printf("%3f ", poFeature->GetFieldAsDouble(iField));
}
else if (poFieldDefn->GetType() == OFTString)
{
feildv3 = poFeature->GetFieldAsString(iField);
// printf("%s ", poFeature->GetFieldAsString(iField));
}
// else
// printf("%s", poFeature->GetFieldAsString(iField));
}
//以上操作与gdal官网一样
//get Geometry获取polygon的信息
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
double x, y;double area;
if (poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon)
{
OGRPolygon *poPolygon = (OGRPolygon * )poGeometry;
area = poPolygon->get_Area();
OGRPoint point;
OGRLinearRing *pOGRLinearRing = poPolygon->getExteriorRing();
pOGRLinearRing->getPoint(0,&point);//得到第一个点的坐标
x = point.getX();
y = point.getY();
// printf("%.6f, %.6f\n", poPoint->getX(), poPoint->getY());
}
else
// printf("no point geometry\n");
OGRFeature::DestroyFeature(poFeature);
}
OGRDataSource::DestroyDataSource(poDS);
}
2.shp文件添加字段
有时需要对已经存在的shp文件进行添加字段操作。代码如下:
OGRDataSource *poDS;
poDS = OGRSFDriverRegistrar::Open("F:/test.shp", TRUE);//第二个参数表示只读还是可读写
if (poDS == NULL)
{
// printf("Open failed \n");
}
OGRLayer *poLayer;
poLayer = poDS->GetLayerByName("test");
int layercount = poDS->GetLayerCount();
OGRFieldDefn poField("NewField", OFTReal);
// poField.SetWidth(1); //设置的是格式化输出时的位数
poField.SetPrecision(8);//设置精度,除了real类型,其他的一般为0
if (poLayer->CreateField(&poField) != OGRERR_NONE)
{
// printf("Creating Name Field Failed\n");
// exit(1);
}
OGRFeature *poFeature;
poLayer->ResetReading();//设置从第0个要素开始读
int featureNum = 0;
while((poFeature = poLayer->GetNextFeature()) != NULL)
{
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
int fieldcount = poFDefn->GetFieldCount();
poFeature->SetField(fieldcount - 1, featureNum);
poLayer->SetFeature(poFeature);
OGRFeature::DestroyFeature(poFeature);
featureNum ++;
}
OGRDataSource::DestroyDataSource(poDS);
刚开始老是添加不成功,最后发现是open函数里第二个参数,如果默认设置为FALSE,则表示只读;需要把它设置为TRUE才可读写。
3.创建shp文件
下面主要说说创建多边形shp文件。创建的方法:使用的WKT格式的字符串来进行创建。也可以使用其他的方式进行创建(参考大神博客:点击打开链接)。但是它是wkt方法的C#、python代码,c++代码对应的函数没查到。所以用另一种方法创建如下:
const char *pszDriverName = "ESRI Shapefile";
OGRSFDriver *poDriver;
OGRRegisterAll();
poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverName);
if (poDriver == NULL)
{
// printf("%s driver not available \n", pszDriverName);
}
OGRDataSource *poDS;
poDS = poDriver->CreateDataSource("F:/zjq/testProject/ogrtest/polygon.shp", NULL);
if (poDS == NULL)
{
// printf("Create shapefile failed\n");
// exit(1);
}
OGRLayer *poLayer;
poLayer = poDS->CreateLayer("polygon1", NULL, wkbPolygon, NULL);
if (poLayer == NULL)
{
// printf("Layer create failed\n");
// exit(1);
}
//下面创建属性表
OGRFieldDefn poFieldID("ID", OFTInteger);//创建ID字段
poLayer->CreateField(&poFieldID);
OGRFieldDefn poFieldname("Name", OFTString);//创建Name字段
poFieldname.SetWidth(32);
poLayer->CreateField(&poFieldname);
double x = 0.0, y = 0.0;
for (int i = 0; i < 10; i++)
{
string szname = "a_" + QString::number(i).toStdString();
OGRFeature *poFeature;
poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
poFeature->SetField(0, i);
poFeature->SetField(1, szname.c_str());
OGRLinearRing ring;
ring.addPoint(0 + 10 * i, 0 + 10 *i);
ring.addPoint(0 +10 * i, 10 + 10 *i);
ring.addPoint(10 + 10 *i, 10 + 10 *i);
ring.addPoint(10 + 10 *i, 0 + 10*i);
ring.closeRings();//首尾点重合形成闭合环
OGRPolygon polygon;
polygon.addRing(&ring);
poFeature->SetGeometry(&polygon);
if ( poLayer->CreateFeature( poFeature) != OGRERR_NONE)
{
// printf("Failed to create feature in shapefile \n");
// exit(1);
}
OGRFeature::DestroyFeature(poFeature);
}
OGRDataSource::DestroyDataSource(poDS);
以上创建代码的结果如下图所示:(QGIS显示)