矢量API教程
本文档旨在使用OGR C++类来从文件中读取和写入数据。强烈建议读者先复习 矢量数据模型 描述关键类及其在OGR中角色的文档。
它还包括C和Python中相应函数的代码片段。
从OGR中读取
为了演示使用OGR进行读取，我们将构造一个小实用程序，用于以逗号分隔格式将点层从OGR数据源转储到stdout。
最初，需要注册所需的所有格式驱动程序。这通常是通过调用 GDALAllRegister() 它注册所有内置到GDAL/OGR中的格式驱动程序。
在C++中：
#include "ogrsf_frmts.h"
int main()
{
GDALAllRegister();
在C中：
#include "gdal.h"
int main()
{
GDALAllRegister();
接下来我们需要打开输入OGR数据源。根据所使用的驱动程序，数据源可以是文件、rdbmse、充满文件的目录，甚至是远程web服务。但是，数据源名称始终是单个字符串。在这种情况下，我们被硬编码为打开一个特定的shapefile。第二个参数（向量的GDAL_）告诉 OGROpen() 方法，我们希望使用向量驱动程序且不需要更新访问。失败时返回NULL，并报告错误。
在C++中：
GDALDataset *poDS;
poDS = (GDALDataset*) GDALOpenEx( "point.shp", GDAL_OF_VECTOR, NULL, NULL, NULL );
if( poDS == NULL )
{
printf( "Open failed.\n" );
exit( 1 );
}
在C中：
GDALDatasetH hDS;
hDS = GDALOpenEx( "point.shp", GDAL_OF_VECTOR, NULL, NULL, NULL );
if( hDS == NULL )
{
printf( "Open failed.\n" );
exit( 1 );
}
GDALDataset可能有许多层与之关联。可以使用查询可用的层数 GDALDataset::GetLayerCount() 以及使用 GDALDataset::GetLayer() . 不过，我们只是按名称获取层。
在C++中：
OGRLayer *poLayer;
poLayer = poDS->GetLayerByName( "point" );
在C中：
OGRLayerH hLayer;
hLayer = GDALDatasetGetLayerByName( hDS, "point" );
现在我们要从层开始读取特性。在我们开始之前，我们可以给图层指定一个属性或空间过滤器来限制我们返回的特征集，但是现在我们对获取所有特征感兴趣。
用GDAL 2.3和C++ 11：
for( auto& poFeature: poLayer )
{
GDAL 2.3和C：
OGR_FOR_EACH_FEATURE_BEGIN(hFeature, hLayer)
{
如果使用旧的GDAL版本，虽然在这种情况下并不是绝对必要的，因为我们是从新的层开始的，但是调用 OGRLayer::ResetReading() 以确保我们从层的开始开始。我们使用OGRLayer：：GetNextFeature（）遍历层中的所有功能。当功能用完时，它将返回空值。
GDAL＜2.3和C++：
OGRFeature *poFeature;
poLayer->ResetReading();
while( (poFeature = poLayer->GetNextFeature()) != NULL )
{
当GDAL<2.3和C时：
OGRFeatureH hFeature;
OGR_L_ResetReading(hLayer);
while( (hFeature = OGR_L_GetNextFeature(hLayer)) != NULL )
{
为了转储特性的所有属性字段，获取 OGRFeatureDefn . 这是一个对象，与图层关联，包含所有字段的定义。我们循环遍历所有字段，并根据其类型获取和报告属性。
用GDAL 2.3和C++ 11：
for( auto&& oField: *poFeature )
{
switch( oField.GetType() )
{
case OFTInteger:
printf( "%d,", oField.GetInteger() );
break;
case OFTInteger64:
printf( CPL_FRMT_GIB ",", oField.GetInteger64() );
break;
case OFTReal:
printf( "%.3f,", oField.GetDouble() );
break;
case OFTString:
printf( "%s,", oField.GetString() );
break;
default:
printf( "%s,", oField.GetAsString() );
break;
}
}
GDAL＜2.3和C++：
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
for( int iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
switch( poFieldDefn->GetType() )
{
case OFTInteger:
printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
break;
case OFTInteger64:
printf( CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64( iField ) );
break;
case OFTReal:
printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
break;
case OFTString:
printf( "%s,", poFeature->GetFieldAsString(iField) );
break;
default:
printf( "%s,", poFeature->GetFieldAsString(iField) );
break;
}
}
在C中：
OGRFeatureDefnH hFDefn = OGR_L_GetLayerDefn(hLayer);
int iField;
for( iField = 0; iField < OGR_FD_GetFieldCount(hFDefn); iField++ )
{
OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, iField );
switch( OGR_Fld_GetType(hFieldDefn) )
{
case OFTInteger:
printf( "%d,", OGR_F_GetFieldAsInteger( hFeature, iField ) );
break;
case OFTInteger64:
printf( CPL_FRMT_GIB ",", OGR_F_GetFieldAsInteger64( hFeature, iField ) );
break;
case OFTReal:
printf( "%.3f,", OGR_F_GetFieldAsDouble( hFeature, iField) );
break;
case OFTString:
printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) );
break;
default:
printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) );
break;
}
}
与上面显式处理的字段类型相比，还有一些字段类型，但是可以使用 OGRFeature::GetFieldAsString() 方法。事实上，我们可以通过对所有类型使用GetFieldAsString（）来缩短上述过程。
接下来，我们要从特征中提取几何图形，并写出点几何图形x和y OGRGeometry 指针。然后我们确定具体的几何类型，如果它是一个点，我们将它转换为一个点并对其进行操作。如果它是其他东西，我们写占位符。
在C++中：
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
{
#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(2,3,0)
OGRPoint *poPoint = poGeometry->toPoint();
#else
OGRPoint *poPoint = (OGRPoint *) poGeometry;
#endif
printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() );
}
else
{
printf( "no point geometry\n" );
}
在C中：
OGRGeometryH hGeometry;
hGeometry = OGR_F_GetGeometryRef(hFeature);
if( hGeometry != NULL
&& wkbFlatten(OGR_G_GetGeometryType(hGeometry)) == wkbPoint )
{
printf( "%.3f,%3.f\n", OGR_G_GetX(hGeometry, 0), OGR_G_GetY(hGeometry, 0) );
}
else
{
printf( "no point geometry\n" );
}
这个 wkbFlatten() 宏用于将wkbPoint25D（具有z坐标的点）的类型转换为基本二维几何类型代码（wkbPoint）。对于每个二维几何图形类型，都有相应的2.5D类型代码。2D和2.5D几何案例由同一C++类处理，所以我们的代码将正确处理2D或3D案例。
可以将多个几何体字段关联到一个特征。
在C++中：
OGRGeometry *poGeometry;
int iGeomField;
int nGeomFieldCount;
nGeomFieldCount = poFeature->GetGeomFieldCount();
for(iGeomField = 0; iGeomField < nGeomFieldCount; iGeomField ++ )
{
poGeometry = poFeature->GetGeomFieldRef(iGeomField);
if( poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
{
#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(2,3,0)
OGRPoint *poPoint = poGeometry->toPoint();
#else
OGRPoint *poPoint = (OGRPoint *) poGeometry;
#endif
printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() );
}
else
{
printf( "no point geometry\n" );
}
}
在C中：
OGRGeometryH hGeometry;
int iGeomField;
int nGeomFieldCount;
nGeomFieldCount = OGR_F_GetGeomFieldCount(hFeature);
for(iGeomField = 0; iGeomField < nGeomFieldCount; iGeomField ++ )
{
hGeometry = OGR_F_GetGeomFieldRef(hFeature, iGeomField);
if( hGeometry != NULL
&& wkbFlatten(OGR_G_GetGeometryType(hGeometry)) == wkbPoint )
{
printf( "%.3f,%3.f\n", OGR_G_GetX(hGeometry, 0),
OGR_G_GetY(hGeometry, 0) );
}
else
{
printf( "no point geometry\n" );
}
}
在 Python 中：
nGeomFieldCount = feat.GetGeomFieldCount()
for iGeomField in range(nGeomFieldCount):
geom = feat.GetGeomFieldRef(iGeomField)
if geom is not None and geom.GetGeometryType() == ogr.wkbPoint:
print "%.3f, %.3f" % ( geom.GetX(), geom.GetY() )
else:
print "no point geometry\n"
注意 OGRFeature::GetGeometryRef() 和 OGRFeature::GetGeomFieldRef() 返回一个指向OGRFeature拥有的内部几何图形的指针。在这里，我们实际上并没有删除返回几何体。
使用GDAL 2.3和C++ 11，循环结束的特征只由关闭的卷曲括号终止。
}
在GDAL 2.3和C中，循环功能的终止方式如下。
}
OGR_FOR_EACH_FEATURE_END(hFeature)
当GDAL<2.3时 OGRLayer::GetNextFeature() 方法返回我们现在拥有的功能的副本。所以在使用结束时我们必须释放这个特性。我们可以“删除”它，但这可能会导致windows构建中的问题，其中GDAL DLL与主程序有不同的“堆”。为了安全起见，我们使用一个GDAL函数来删除这个特性。
在C++中：
OGRFeature::DestroyFeature( poFeature );
}
在C中：
OGR_F_Destroy( hFeature );
}
OGRLayer返回 GDALDataset::GetLayerByName() 也是对GDALDataset拥有的内部层的引用，因此我们不需要删除它。但我们确实需要删除数据源才能关闭输入文件。我们再次使用自定义的delete方法来避免特殊的win32堆问题。
在C/C++中：
GDALClose( poDS );
}
我们的程序看起来是这样的。
用GDAL 2.3和C++ 11：
#include "ogrsf_frmts.h"
int main()
{
GDALAllRegister();
GDALDatasetUniquePtr poDS(GDALDataset::Open( "point.shp", GDAL_OF_VECTOR));
if( poDS == nullptr )
{
printf( "Open failed.\n" );
exit( 1 );
}
for( const OGRLayer* poLayer: poDS->GetLayers() )
{
for( const auto& poFeature: *poLayer )
{
for( const auto& oField: *poFeature )
{
switch( oField.GetType() )
{
case OFTInteger:
printf( "%d,", oField.GetInteger() );
break;
case OFTInteger64:
printf( CPL_FRMT_GIB ",", oField.GetInteger64() );
break;
case OFTReal:
printf( "%.3f,", oField.GetDouble() );
break;
case OFTString:
printf( "%s,", oField.GetString() );
break;
default:
printf( "%s,", oField.GetAsString() );
break;
}
}
const OGRGeometry *poGeometry = poFeature->GetGeometryRef();
if( poGeometry != nullptr
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
{
const OGRPoint *poPoint = poGeometry->toPoint();
printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() );
}
else
{
printf( "no point geometry\n" );
}
}
}
return 0;
}
在C++中：
#include "ogrsf_frmts.h"
int main()
{
GDALAllRegister();
GDALDataset *poDS = static_cast<GDALDataset*>(
GDALOpenEx( "point.shp", GDAL_OF_VECTOR, NULL, NULL, NULL ));
if( poDS == NULL )
{
printf( "Open failed.\n" );
exit( 1 );
}
OGRLayer *poLayer = poDS->GetLayerByName( "point" );
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
poLayer->ResetReading();
OGRFeature *poFeature;
while( (poFeature = poLayer->GetNextFeature()) != NULL )
{
for( int iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
switch( poFieldDefn->GetType() )
{
case OFTInteger:
printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
break;
case OFTInteger64:
printf( CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64( iField ) );
break;
case OFTReal:
printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
break;
case OFTString:
printf( "%s,", poFeature->GetFieldAsString(iField) );
break;
default:
printf( "%s,", poFeature->GetFieldAsString(iField) );
break;
}
}
OGRGeometry *poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
{
OGRPoint *poPoint = (OGRPoint *) poGeometry;
printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() );
}
else
{
printf( "no point geometry\n" );
}
OGRFeature::DestroyFeature( poFeature );
}
GDALClose( poDS );
}
在C中：
#include "gdal.h"
int main()
{
GDALAllRegister();
GDALDatasetH hDS;
OGRLayerH hLayer;
OGRFeatureH hFeature;
OGRFeatureDefnH hFDefn;
hDS = GDALOpenEx( "point.shp", GDAL_OF_VECTOR, NULL, NULL, NULL );
if( hDS == NULL )
{
printf( "Open failed.\n" );
exit( 1 );
}
hLayer = GDALDatasetGetLayerByName( hDS, "point" );
hFDefn = OGR_L_GetLayerDefn(hLayer);
OGR_L_ResetReading(hLayer);
while( (hFeature = OGR_L_GetNextFeature(hLayer)) != NULL )
{
int iField;
OGRGeometryH hGeometry;
for( iField = 0; iField < OGR_FD_GetFieldCount(hFDefn); iField++ )
{
OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, iField );
switch( OGR_Fld_GetType(hFieldDefn) )
{
case OFTInteger:
printf( "%d,", OGR_F_GetFieldAsInteger( hFeature, iField ) );
break;
case OFTInteger64:
printf( CPL_FRMT_GIB ",", OGR_F_GetFieldAsInteger64( hFeature, iField ) );
break;
case OFTReal:
printf( "%.3f,", OGR_F_GetFieldAsDouble( hFeature, iField) );
break;
case OFTString:
printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) );
break;
default:
printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) );
break;
}
}
hGeometry = OGR_F_GetGeometryRef(hFeature);
if( hGeometry != NULL
&& wkbFlatten(OGR_G_GetGeometryType(hGeometry)) == wkbPoint )
{
printf( "%.3f,%3.f\n", OGR_G_GetX(hGeometry, 0), OGR_G_GetY(hGeometry, 0) );
}
else
{
printf( "no point geometry\n" );
}
OGR_F_Destroy( hFeature );
}
GDALClose( hDS );
}
在 Python 中：
import sys
from osgeo import gdal
ds = gdal.OpenEx( "point.shp", gdal.OF_VECTOR )
if ds is None:
print "Open failed.\n"
sys.exit( 1 )
lyr = ds.GetLayerByName( "point" )
lyr.ResetReading()
for feat in lyr:
feat_defn = lyr.GetLayerDefn()
for i in range(feat_defn.GetFieldCount()):
field_defn = feat_defn.GetFieldDefn(i)
# Tests below can be simplified with just :
# print feat.GetField(i)
if field_defn.GetType() == ogr.OFTInteger or field_defn.GetType() == ogr.OFTInteger64:
print "%d" % feat.GetFieldAsInteger64(i)
elif field_defn.GetType() == ogr.OFTReal:
print "%.3f" % feat.GetFieldAsDouble(i)
elif field_defn.GetType() == ogr.OFTString:
print "%s" % feat.GetFieldAsString(i)
else:
print "%s" % feat.GetFieldAsString(i)
geom = feat.GetGeometryRef()
if geom is not None and geom.GetGeometryType() == ogr.wkbPoint:
print "%.3f, %.3f" % ( geom.GetX(), geom.GetY() )
else:
print "no point geometry\n"
ds = None
写信给OGR
作为通过OGR编写的一个例子，我们将做与上面大致相反的事情。从输入文本中读取逗号分隔值的短程序将通过OGR写入点形状文件。
像往常一样，我们首先注册所有驱动程序，然后获取Shapefile驱动程序，因为我们需要它来创建输出文件。
在C++中：
#include "ogrsf_frmts.h"
int main()
{
const char *pszDriverName = "ESRI Shapefile";
GDALDriver *poDriver;
GDALAllRegister();
poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );
if( poDriver == NULL )
{
printf( "%s driver not available.\n", pszDriverName );
exit( 1 );
}
在C中：
#include "ogr_api.h"
int main()
{
const char *pszDriverName = "ESRI Shapefile";
GDALDriver *poDriver;
GDALAllRegister();
poDriver = (GDALDriver*) GDALGetDriverByName(pszDriverName );
if( poDriver == NULL )
{
printf( "%s driver not available.\n", pszDriverName );
exit( 1 );
}
接下来我们创建数据源。ESRI Shapefile驱动程序允许我们创建一个包含Shapefile的目录，或者创建一个Shapefile作为数据源。在本例中，我们将通过在名称中包含扩展名来显式创建单个文件。其他司机的行为也不同。第二、第三、第四和第五个参数与栅格尺寸有关（如果驱动程序具有栅格功能）。调用的最后一个参数是一个选项值列表，但在本例中我们只使用默认值。支持的选项的详细信息也是特定于格式的。
在C++中：
GDALDataset *poDS;
poDS = poDriver->Create( "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );
if( poDS == NULL )
{
printf( "Creation of output file failed.\n" );
exit( 1 );
}
在C中：
GDALDatasetH hDS;
hDS = GDALCreate( hDriver, "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );
if( hDS == NULL )
{
printf( "Creation of output file failed.\n" );
exit( 1 );
}
现在我们创建输出层。在这种情况下，由于数据源是单个文件，因此我们只能有一个层。我们传递wkbPoint来指定此层支持的几何体类型。在这种情况下，我们不会传递任何坐标系信息或其他特殊的图层创建选项。
在C++中：
OGRLayer *poLayer;
poLayer = poDS->CreateLayer( "point_out", NULL, wkbPoint, NULL );
if( poLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
在C中：
OGRLayerH hLayer;
hLayer = GDALDatasetCreateLayer( hDS, "point_out", NULL, wkbPoint, NULL );
if( hLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
既然层已经存在，我们需要创建任何应该出现在层上的属性字段。在写入任何功能之前，必须将字段添加到层中。要创建字段，我们初始化 OGRField 对象，其中包含有关字段的信息。在shapefile的情况下，字段宽度和精度在创建output.dbf文件时非常重要，因此我们专门设置它，尽管通常默认值是可以的。对于这个例子，我们只有一个属性，一个与x，y点相关联的名称字符串。
注意，我们传递到的模板OGRField OGRLayer::CreateField() 在内部复制。我们保留该物品的所有权。
在C++中：
OGRFieldDefn oField( "Name", OFTString );
oField.SetWidth(32);
if( poLayer->CreateField( &oField ) != OGRERR_NONE )
{
printf( "Creating Name field failed.\n" );
exit( 1 );
}
在C中：
OGRFieldDefnH hFieldDefn;
hFieldDefn = OGR_Fld_Create( "Name", OFTString );
OGR_Fld_SetWidth( hFieldDefn, 32);
if( OGR_L_CreateField( hLayer, hFieldDefn, TRUE ) != OGRERR_NONE )
{
printf( "Creating Name field failed.\n" );
exit( 1 );
}
OGR_Fld_Destroy(hFieldDefn);
下面的剪循环从stdin读取“x，y，name”形式的行并对它们进行解析。
在C++和C语言中：
double x, y;
char szName[33];
while( !feof(stdin)
&& fscanf( stdin, "%lf,%lf,%32s", &x, &y, szName ) == 3 )
{
要将要素写入磁盘，必须先创建本地OGRFeature，设置属性并附加几何图形，然后再尝试将其写入层。必须从与要写入的层关联的OGRFeatureDefn实例化此功能。
在C++中：
OGRFeature *poFeature;
poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
poFeature->SetField( "Name", szName );
在C中：
OGRFeatureH hFeature;
hFeature = OGR_F_Create( OGR_L_GetLayerDefn( hLayer ) );
OGR_F_SetFieldString( hFeature, OGR_F_GetFieldIndex(hFeature, "Name"), szName );
我们创建一个局部几何体对象，并将其副本（间接）指定给特征。这个 OGRFeature::SetGeometryDirectly() 不同于 OGRFeature::SetGeometry() 在这种情况下，直接方法将几何图形的所有权赋予特征。这通常更有效，因为它避免了几何体的超深对象副本。
在C++中：
OGRPoint pt;
pt.setX( x );
pt.setY( y );
poFeature->SetGeometry( &pt );
在C中：
OGRGeometryH hPt;
hPt = OGR_G_CreateGeometry(wkbPoint);
OGR_G_SetPoint_2D(hPt, 0, x, y);
OGR_F_SetGeometry( hFeature, hPt );
OGR_G_DestroyGeometry(hPt);
现在我们在文件中创建一个特性。这个 OGRLayer::CreateFeature() 不会拥有我们的功能，所以我们在完成后会清理它。
在C++中：
if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature in shapefile.\n" );
exit( 1 );
}
OGRFeature::DestroyFeature( poFeature );
}
在C中：
if( OGR_L_CreateFeature( hLayer, hFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature in shapefile.\n" );
exit( 1 );
}
OGR_F_Destroy( hFeature );
}
最后，我们需要关闭数据源，以确保以有序的方式写出头并恢复所有资源。
在C/C++中：
GDALClose( poDS );
}
一个块中的同一程序如下所示：
在C++中：
#include "ogrsf_frmts.h"
int main()
{
const char *pszDriverName = "ESRI Shapefile";
GDALDriver *poDriver;
GDALAllRegister();
poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );
if( poDriver == NULL )
{
printf( "%s driver not available.\n", pszDriverName );
exit( 1 );
}
GDALDataset *poDS;
poDS = poDriver->Create( "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );
if( poDS == NULL )
{
printf( "Creation of output file failed.\n" );
exit( 1 );
}
OGRLayer *poLayer;
poLayer = poDS->CreateLayer( "point_out", NULL, wkbPoint, NULL );
if( poLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
OGRFieldDefn oField( "Name", OFTString );
oField.SetWidth(32);
if( poLayer->CreateField( &oField ) != OGRERR_NONE )
{
printf( "Creating Name field failed.\n" );
exit( 1 );
}
double x, y;
char szName[33];
while( !feof(stdin)
&& fscanf( stdin, "%lf,%lf,%32s", &x, &y, szName ) == 3 )
{
OGRFeature *poFeature;
poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
poFeature->SetField( "Name", szName );
OGRPoint pt;
pt.setX( x );
pt.setY( y );
poFeature->SetGeometry( &pt );
if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature in shapefile.\n" );
exit( 1 );
}
OGRFeature::DestroyFeature( poFeature );
}
GDALClose( poDS );
}
在C中：
#include "gdal.h"
int main()
{
const char *pszDriverName = "ESRI Shapefile";
GDALDriverH hDriver;
GDALDatasetH hDS;
OGRLayerH hLayer;
OGRFieldDefnH hFieldDefn;
double x, y;
char szName[33];
GDALAllRegister();
hDriver = GDALGetDriverByName( pszDriverName );
if( hDriver == NULL )
{
printf( "%s driver not available.\n", pszDriverName );
exit( 1 );
}
hDS = GDALCreate( hDriver, "point_out.shp", 0, 0, 0, GDT_Unknown, NULL );
if( hDS == NULL )
{
printf( "Creation of output file failed.\n" );
exit( 1 );
}
hLayer = GDALDatasetCreateLayer( hDS, "point_out", NULL, wkbPoint, NULL );
if( hLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
hFieldDefn = OGR_Fld_Create( "Name", OFTString );
OGR_Fld_SetWidth( hFieldDefn, 32);
if( OGR_L_CreateField( hLayer, hFieldDefn, TRUE ) != OGRERR_NONE )
{
printf( "Creating Name field failed.\n" );
exit( 1 );
}
OGR_Fld_Destroy(hFieldDefn);
while( !feof(stdin)
&& fscanf( stdin, "%lf,%lf,%32s", &x, &y, szName ) == 3 )
{
OGRFeatureH hFeature;
OGRGeometryH hPt;
hFeature = OGR_F_Create( OGR_L_GetLayerDefn( hLayer ) );
OGR_F_SetFieldString( hFeature, OGR_F_GetFieldIndex(hFeature, "Name"), szName );
hPt = OGR_G_CreateGeometry(wkbPoint);
OGR_G_SetPoint_2D(hPt, 0, x, y);
OGR_F_SetGeometry( hFeature, hPt );
OGR_G_DestroyGeometry(hPt);
if( OGR_L_CreateFeature( hLayer, hFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature in shapefile.\n" );
exit( 1 );
}
OGR_F_Destroy( hFeature );
}
GDALClose( hDS );
}
在Python中：
import sys
from osgeo import gdal
from osgeo import ogr
import string
driverName = "ESRI Shapefile"
drv = gdal.GetDriverByName( driverName )
if drv is None:
print "%s driver not available.\n" % driverName
sys.exit( 1 )
ds = drv.Create( "point_out.shp", 0, 0, 0, gdal.GDT_Unknown )
if ds is None:
print "Creation of output file failed.\n"
sys.exit( 1 )
lyr = ds.CreateLayer( "point_out", None, ogr.wkbPoint )
if lyr is None:
print "Layer creation failed.\n"
sys.exit( 1 )
field_defn = ogr.FieldDefn( "Name", ogr.OFTString )
field_defn.SetWidth( 32 )
if lyr.CreateField ( field_defn ) != 0:
print "Creating Name field failed.\n"
sys.exit( 1 )
# Expected format of user input: x y name
linestring = raw_input()
linelist = string.split(linestring)
while len(linelist) == 3:
x = float(linelist[0])
y = float(linelist[1])
name = linelist[2]
feat = ogr.Feature( lyr.GetLayerDefn())
feat.SetField( "Name", name )
pt = ogr.Geometry(ogr.wkbPoint)
pt.SetPoint_2D(0, x, y)
feat.SetGeometry(pt)
if lyr.CreateFeature(feat) != 0:
print "Failed to create feature in shapefile.\n"
sys.exit( 1 )
feat.Destroy()
linestring = raw_input()
linelist = string.split(linestring)
ds = None
可以将多个几何体字段关联到一个特征。此功能仅适用于一些文件格式，如PostGIS。
要创建此类数据源，必须首先创建几何字段。空间参考系统对象可以关联到每个几何体字段。
在C++中：
OGRGeomFieldDefn oPointField( "PointField", wkbPoint );
OGRSpatialReference* poSRS = new OGRSpatialReference();
poSRS->importFromEPSG(4326);
oPointField.SetSpatialRef(poSRS);
poSRS->Release();
if( poLayer->CreateGeomField( &oPointField ) != OGRERR_NONE )
{
printf( "Creating field PointField failed.\n" );
exit( 1 );
}
OGRGeomFieldDefn oFieldPoint2( "PointField2", wkbPoint );
poSRS = new OGRSpatialReference();
poSRS->importFromEPSG(32631);
oPointField2.SetSpatialRef(poSRS);
poSRS->Release();
if( poLayer->CreateGeomField( &oPointField2 ) != OGRERR_NONE )
{
printf( "Creating field PointField2 failed.\n" );
exit( 1 );
}
在C中：
OGRGeomFieldDefnH hPointField;
OGRGeomFieldDefnH hPointField2;
OGRSpatialReferenceH hSRS;
hPointField = OGR_GFld_Create( "PointField", wkbPoint );
hSRS = OSRNewSpatialReference( NULL );
OSRImportFromEPSG(hSRS, 4326);
OGR_GFld_SetSpatialRef(hPointField, hSRS);
OSRRelease(hSRS);
if( OGR_L_CreateGeomField( hLayer, hPointField ) != OGRERR_NONE )
{
printf( "Creating field PointField failed.\n" );
exit( 1 );
}
OGR_GFld_Destroy( hPointField );
hPointField2 = OGR_GFld_Create( "PointField2", wkbPoint );
OSRImportFromEPSG(hSRS, 32631);
OGR_GFld_SetSpatialRef(hPointField2, hSRS);
OSRRelease(hSRS);
if( OGR_L_CreateGeomField( hLayer, hPointField2 ) != OGRERR_NONE )
{
printf( "Creating field PointField2 failed.\n" );
exit( 1 );
}
OGR_GFld_Destroy( hPointField2 );
要将功能写入磁盘，必须先创建本地OGRFeature、设置属性并附加几何体，然后再尝试将其写入层。必须从与要写入的层关联的OGRFeatureDefn实例化此功能。
在C++中：
OGRFeature *poFeature;
OGRGeometry *poGeometry;
char* pszWKT;
poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
pszWKT = (char*) "POINT (2 49)";
OGRGeometryFactory::createFromWkt( &pszWKT, NULL, &poGeometry );
poFeature->SetGeomFieldDirectly( "PointField", poGeometry );
pszWKT = (char*) "POINT (500000 4500000)";
OGRGeometryFactory::createFromWkt( &pszWKT, NULL, &poGeometry );
poFeature->SetGeomFieldDirectly( "PointField2", poGeometry );
if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature.\n" );
exit( 1 );
}
OGRFeature::DestroyFeature( poFeature );
在C中：
OGRFeatureH hFeature;
OGRGeometryH hGeometry;
char* pszWKT;
poFeature = OGR_F_Create( OGR_L_GetLayerDefn(hLayer) );
pszWKT = (char*) "POINT (2 49)";
OGR_G_CreateFromWkt( &pszWKT, NULL, &hGeometry );
OGR_F_SetGeomFieldDirectly( hFeature,
OGR_F_GetGeomFieldIndex(hFeature, "PointField"), hGeometry );
pszWKT = (char*) "POINT (500000 4500000)";
OGR_G_CreateFromWkt( &pszWKT, NULL, &hGeometry );
OGR_F_SetGeomFieldDirectly( hFeature,
OGR_F_GetGeomFieldIndex(hFeature, "PointField2"), hGeometry );
if( OGR_L_CreateFeature( hFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature.\n" );
exit( 1 );
}
OGR_F_Destroy( hFeature );
在Python中：
feat = ogr.Feature( lyr.GetLayerDefn() )
feat.SetGeomFieldDirectly( "PointField",
ogr.CreateGeometryFromWkt( "POINT (2 49)" ) )
feat.SetGeomFieldDirectly( "PointField2",
ogr.CreateGeometryFromWkt( "POINT (500000 4500000)" ) )
if lyr.CreateFeature( feat ) != 0:
print( "Failed to create feature.\n" );
sys.exit( 1 );