Skip to content

Latest commit

 

History

History
411 lines (353 loc) · 13.4 KB

ogr-read-write.rst

File metadata and controls

411 lines (353 loc) · 13.4 KB
Raw

rst

矢量数据读写

数据读取

与栅格的 gdalread 流程类似，但矢量数据定位明显，属性隐含，其几何形状和属性需要分开读取，GDAL有两个版本，读取接口不同，但是流程基本相同，以下是获取流程：

  • 注册驱动
  • 打开数据集
  • 打开图层
  • 获取要素类定义
  • 遍历要素
  • 获取几何形状以及属性信息

1.x 版本完整代码示例：

#include "ogrsf_frmts.h"
int main()
{
    //注册驱动
    OGRRegisterAll();
    //打开数据集
    OGRDataSource       *poDS;
    poDS = OGRSFDriverRegistrar::Open( "point.shp", FALSE );
    if( poDS == NULL )
    {
        printf( "Open failed.\n" );
        exit( 1 );
    }
    //打开图层
    OGRLayer  *poLayer;
    poLayer = poDS->GetLayerByName( "point" );
    OGRFeature *poFeature;
    //遍历要素
    poLayer->ResetReading();
    while( (poFeature = poLayer->GetNextFeature()) != NULL )
    {
        //获取要素类定义
        OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
        int iField;
        //获取属性信息
        for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
        {
            OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
            if( poFieldDefn->GetType() == OFTInteger )
                printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
            else if( poFieldDefn->GetType() == OFTReal )
                printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
            else if( poFieldDefn->GetType() == OFTString )
                printf( "%s,", poFeature->GetFieldAsString(iField) );
            else
                printf( "%s,", poFeature->GetFieldAsString(iField) );
        }
        //获取几何形状
        OGRGeometry *poGeometry;
        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 );
    }
    OGRDataSource::DestroyDataSource( poDS );
}

2.x 版本中最主要修改为 GDALDataset 代替OGRDataSource 其余接口大致不变.打开方式略有区别

/**
* @param pszFilename 文件名
* @param nOpenFlags  GDAL_OF_RASTER/GDAL_OF_VECTOR  数据类型
                     GDAL_OF_READONLY/GDAL_OF_UPDATE 打开方式
                     GDAL_OF_SHARED 共享模式
                     GDAL_OF_VERBOSE_ERROR 返回详细错误
* @param papszAllowedDrivers NULL或以NULL结尾的列表,指定驱动
* @param papszOpenOptions NULL或以NULL结尾的字符串列表,根据驱动设置
* @param papszSiblingFiles NULL或以NULL结尾的文件列表
*/
GDALDatasetH GDALOpenEx (   const char *    pszFilename,
unsigned int    nOpenFlags,
const char *const *     papszAllowedDrivers,
const char *const *     papszOpenOptions,
const char *const *     papszSiblingFiles 
)

Attention

  • 打开前建议添加 CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); 设置,保证windows下文件名中文路径正常打开
  • 打开前添加 CPLSetConfigOption("SHAPE_ENCODING", ""); 读出来编码和系统编码保持一致

2.x 版本完整代码示例：

#include "ogrsf_frmts.h"
int main()
{
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    CPLSetConfigOption("SHAPE_ENCODING", "");
    //主要修改
    GDALAllRegister();
    GDALDataset       *poDS;
    poDS = (GDALDataset*) GDALOpenEx( "point.shp", GDAL_OF_VECTOR, NULL, NULL, NULL );
    if( poDS == NULL )
    {
        printf( "Open failed.\n" );
        exit( 1 );
    }

    //后面与1.x版本一致
    OGRLayer  *poLayer;
    poLayer = poDS->GetLayerByName( "point" );
    //寻找特定字段
    OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();

    int nSIDIndex = poFDefn->GetFieldIndex("SS_ID");
    if (nSIDIndex < 0)
    {
        printf( "can't find SS_ID field.\n" );
    }

    OGRFeature *poFeature;
    poLayer->ResetReading();
    while( (poFeature = poLayer->GetNextFeature()) != NULL )
    {
        //获取固定字段值
        if(nSIDIndex >= 0)
        {
            std::string strSSID = poFeature->GetFieldAsString(nSIDIndex);
        }
        OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
        int iField;
        for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
        {
            OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
            if( poFieldDefn->GetType() == OFTInteger )
                printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
            else if( poFieldDefn->GetType() == OFTInteger64 )
                printf( CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64( iField ) );
            else if( poFieldDefn->GetType() == OFTReal )
                printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
            else if( poFieldDefn->GetType() == OFTString )
                printf( "%s,", poFeature->GetFieldAsString(iField) );
            else
                printf( "%s,", poFeature->GetFieldAsString(iField) );
        }
        OGRGeometry *poGeometry;
        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 );
}

Attention

  • 有很多矢量驱动是只读的，无法像栅格那样直接修改和创建文件。
  • 如果需要确定，可以使用 ogrinfo --formats 确定驱动的权限。

数据写入

数据写入与读取类似，流程如下：

  • 注册驱动
  • 打开或创建数据集
  • 创建图层
  • 创建要素类定义
  • 新建几何形状
  • 新建要素
  • 关闭数据集

1.x完整代码示例（ C++）：

#include "ogrsf_frmts.h"
int main()
{
    const char *pszDriverName = "ESRI Shapefile";
    OGRSFDriver *poDriver;
    OGRRegisterAll();
    poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(
                pszDriverName );
    if( poDriver == NULL )
    {
        printf( "%s driver not available.\n", pszDriverName );
        exit( 1 );
    }
    OGRDataSource *poDS;
    poDS = poDriver->CreateDataSource( "point_out.shp", 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 );
    }
    OGRDataSource::DestroyDataSource( poDS );
}

2.x 完整代码示例（ 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 );
}

数据修改

GDAL可以添加、删除、修改属性信息和要素，下面简单介绍下如何添加/编辑属性字段(2.x)：

GDALDataset* dst = (GDALDataset*)GDALOpenEx(shpName, GDAL_OF_UPDATE, NULL, NULL, NULL);
OGRLayer  *poLayer = dst->GetLayer(0);
OGRFieldDefn oField("imgNo", OFTString);  // Add fields
oField.SetWidth(128);
if (poLayer->CreateField(&oField, TRUE) != OGRERR_NONE)
{
printf("Creating Name field failed.\n");
}
OGRFeature *poFeature;

poLayer->ResetReading();
while ((poFeature = poLayer->GetNextFeature()) != NULL)
{
    poFeature->SetField("imgNo", "123456");  // Set fields
    poLayer->SetFeature(poFeature);
    OGRFeature::DestroyFeature(poFeature);
}
GDALClose(dst);

删除其实也类似,需要注意打开时一定要加上 GDAL_OF_UPDATE 设置

CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
CPLSetConfigOption("SHAPE_ENCODING", "");
GDALAllRegister();
GDALDataset *poDS = (GDALDataset*)GDALOpenEx(strIn.c_str(),\
 GDAL_OF_VECTOR| GDAL_OF_UPDATE, NULL, NULL, NULL);

OGRLayer* pLayer = poDS->GetLayer(0);
if (!pLayer->TestCapability(OLCDeleteFeature))
{
    LOGI( "Layer does not support delete feature capability");
}
pLayer->DeleteFeature(1);
pLayer->DeleteFeature(2);
pLayer->DeleteFeature(3);
pLayer->DeleteFeature(4);
pLayer->DeleteFeature(5);
pLayer->SyncToDisk();
std::stringstream sql;
//GDAL版本小于2.2时,shapefile需要这一步才能真正写入
sql << "REPACK " << pLayer->GetName();
poDS->ExecuteSQL(sql.str().c_str(), NULL, NULL);
GDALClose(poDS);

Attention

  • GDAL版本小于2.2时,处理shp文件时记得最后要执行 REPACK 表名 ,否则只是临时标记,不会真正删除。
  • GDAL 2.2之后的版本在调用 FlushCache() 或 SyncToDisk() 时会自动执行 REPACK

创建空间索引

有时候合并和处理比较慢,可以尝试先建立空间索引再处理,建空间索引方式跟执行sql一样,调用 GDALDataset::ExecuteSQL 方法

  • 删除空间索引 DROP SPATIAL INDEX ON tablename
  • 创建空间索引 CREATE SPATIAL INDEX ON tablename [DEPTH N] N从1-12,指示索引树深度
GDALDataset *poDS = (GDALDataset*)GDALOpenEx(strIn.c_str(),\
 GDAL_OF_VECTOR| GDAL_OF_SHARED, NULL, NULL, NULL);
 OGRLayer* pLayer = poDS->->GetLayer(0);
std::string strCreateIndex = "CREATE SPATIAL INDEX ON ";
strCreateIndex += pLayer->GetName();
//执行创建空间索引的返回是NULL,不需要判断这里
OGRLayer *idxBase = dsBase->ExecuteSQL(strCreateBase.ToPlatform().c_str(), NULL, NULL);
...

Attention

  • 如果表名有空格,在sql里用转义的双引号括起来即可