前言
gdal库是什么？gdal的官网上是这么写的：GDAL 是一个用于栅格和矢量地理空间数据格式的转换器库，由开源地理空间基金会根据 MIT 风格的开源许可证发布。作为一个库，它向调用应用程序提供所有支持的格式的单个光栅抽象数据模型和单个向量抽象数据模型。它还带有各种有用的命令行实用程序，用于数据转换和处理。
事实上gdal已经成为了gis行业的一种标准和基石，arcgis，qgis等gis软件无不在底层使用了gdal库（可能对gdal库进行了一定修改）。gdal库本身也支持各种语言的调用，留出了众多api接口。
作为一个地信（遥感）专业学生，在我的学习经历中，第一次接触gdal库是在大二的一次gis二次开发大赛中，我需要去对栅格图像进行一个TDVI指数的提取，需要在python环境中使用gdal库去读取一副tiff格式的卫星影像并输出一个tiff结果的影像，并保持原有的坐标系信息不丢失，因此不能使用opencv等库把卫星图像当普通图像打开（会丢失坐标信息）。不得不说使用python调包真的很方便，那个时候根本就没想过为什么要用c++去实现。但是在实现算法后想打包成exe的时候才深刻体会到了python的痛点，使用pyinstaller打包时会将整个python环境附带第三方库全部打包，结果就是生成的exe体积极大，运行极慢，尤其当使用到numpy和pandas等科学计算库时，运行时要耗费几分钟去初始化。这就不得不提到那个老生常谈的问题：python和c++选哪个语言好？
Python vs C++
我觉得关于哪个语言问题好的问题其实是很流氓的，每个语言有它的优点和短处，根据实际的情况和场景选择合适的语言去做就好了。
对于一个地学而非计算机的学生，我觉得本科时就把python用好用熟就已经很好了，毕竟本科阶段不会接触到什么实际项目，很少涉及到打包，部署等工程问题。很多时候就是面向数据的，负责分析的，输出一个矢量或栅格或统计的结果，写个小论文/报告就完事了，另一方面本科时间有限，学习c++这种复杂深奥的语言有时候真的来不及（对大部分学生来说，光是c++基础可能就要花一年以上时间来初步掌握，之后更要花数年时间来实践，踩坑）。再加上这几年深度学习实在是太火了，不管你是做什么的都得懂点深度学习，因此本科时候还是得学点python，调调pytorch/paddle/tf等框架对深度学习有个基本的了解，现在保研考研面试的时候老师也太喜欢问有没有深度学习经验了，个人认为深度学习真的是本科时性价比最高的学习方向，当然等研究生之后还要不要卷这个方向另说（反正我是不卷了，否则也不会来钻gdal这个方向）。
回到语言本身，我们先说说python的优缺点
优点：
- 学习成本低，语法糖多，上手快
- 写码速度快，python五行代码写完的c++可能要几十行，有个想法能快速编程实现
- 包/库的配置极为方便（pip install/conda install），且种类丰富，基本所有领域都有库可以直接调用，配合anaconda配置环境非常方便
缺点：
- 运行效率低，封装的太好很难去底层优化。（numpy和pandas等科学库的优化相当好，但是不是所有运算都可以用这些库的api来实现，一旦不能通过矩阵运算来实现而要进行循环遍历时就效率感人了。还有就是python对多线程也不友好）
- 打包是个难题，上面已经说过了
- 对大型项目不够友好（这方面我也没经验，对于地学生来说也影响不大）
接下来说说c++，个人觉得其实c++就凭两点特性就可以活下去：
1运行效率高
2对c的兼容，只要系统底层还是c写的，c++一时半会就死不了
除此之外在gis桌面软件开发领域不管是基于mfc框架还是qt框架，都还是c++实现。
至于c++的缺点，那就不用多说了，对于地学生来说光是编译个库就是个大麻烦，至于指针，地址这种东西也不好理解，更不要谈每三年更新一次版本里的新特性，从智能指针，模板到后面的auto，类型推导，完美转发，一堆说了脑袋都大的知识点。但真当你有那个需求的时候，c++就成了唯一选择，所以纵使有什么缺点，你也得捏着鼻子用（说不定用着用着就爱上这门语言了呢？）。
GDAL编译及导入
对于python，gdal库直接pip install大概率失败，最好通过whl文件来安装。这方面教程很多，不细说，打开cmd，到python目录下（有环境变量就不用），输入pip install，然后把whl文件拖进来就会自动填入whl文件的地址，按回车安装就完事了。导入的话gdal库是在osgeo库中，导入语句是from osgeo import gdal
c++编译和导入
参考这位博主
VS2019C++编译GDAL3.3.2+SQLite3+PROJ6+GEOS3.7.3+HDF4+HDF5（保姆级教程）_gdal2.4.2 sqlite3_Jack丶Wang的博客-CSDN博客
需要注意的是如果同时编译debug和release版本，需要编译好一个版本使用clean命令进行清理，在编译另一个版本。另外分清前置依赖库的debug和release版本。
导入或如果报错缺少什么dll，就去相应的库的编译文件中找，再拖到编译出来的gdal库的bin文件夹中就好。
vs2019中创建gdal项目，并设置gitee仓库
vs2019使用Gitee_gitee vs2019_故里2130的博客-CSDN博客
仓库地址gislxz/GDAL_Learning - 码云 - 开源中国 (gitee.com)
GDAL坐标系与投影
坐标系基础知识
大地水准面-地球椭圆体-基准面-地理坐标系-投影坐标系
OGR空间参考简介
定义地理坐标系
void CreateMyCRS() {
OGRSpatialReference oSRS;
oSRS.SetGeogCS("My geographic CRS",
"World Geodetic System 1984",
"My WGS84 Spheroid",
SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
"Greenwich", 0.0,
"degree", atof(SRS_UA_DEGREE_CONV));
}在上面的代码中，"My geographic CRS"为定义的地理坐标系的名称，"My WGS84 Spheroid"为定义的椭球体的名字，"Greenwich"和"degree"分别表示格林尼治子午线及其单位，这些主要是用来说明这个地理坐标系的一些信息，是展现给用户的而不是什么关键词（key）。然而，"World Geodetic System 1984"是一个定义大地基准的关键词，注意：这里的大地基准必须是一个有效的大地基准，不能随便指定。而上面代码中的SRS_WGS84_SEMIMAJOR、SRS_WGS84_INVFLATTENING和SRS_UA DEGREE_CONV都是OGR定义的一些常数，分别表示WGS84椭球体的长半轴（6378137 米）、WGS84 椭球体扁率的倒数（298.257223563）以及度和弧度的换算比例(0.0174532925199433)。
在我的gdal版本中报错，提示SRS_UA_DEGREE_CONV为const char *而形参类型是double，SRS_UA_DEGREE_CONV是geos库中宏定义的，用atof转成double后不报错。
看了一下原函数的形参，最后一个参数是弧度转角度的比值，对应的就是这个SRS_UA_DEGREE_CONV，不知是什么情况。
OGRSpatialReference内置了对一些众所周知的CRS的支持， 包括“NAD27”、“NAD83”、“WGS72”和“WGS84”，可以是 在对 OGRSpatialReference：：SetWellKnownGeogCS（） 的单个调用中定义。
oSRS.SetWellKnownGeogCS( "WGS84" );如果 EPSG 数据库可用，可用 GCS 代码编号设置。
oSRS.SetWellKnownGeogCS( "EPSG:4326" );为了方便地和其他库进行交互，OGRSpatialReference提供了与OpenGIS的WKT格式的相互转换。OGRSpatialReference 可以使用一个WKT来进行初始化，或者将里面的信息导出为WKT格式。
char *pszWKT = NULL;
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.exportToWkt( &pszWKT );
printf( "%s\n", pszWKT );
CPLFree(pszWKT);gdal3.0+版本可以设置参数
char *pszWKT = nullptr;
oSRS.SetWellKnownGeogCS( "WGS84" );
const char* apszOptions[] = { "FORMAT=WKT2_2018", "MULTILINE=YES", nullptr };
oSRS.exportToWkt( &pszWKT, apszOptions );
printf( "%s\n", pszWKT );
CPLFree(pszWKT);定义投影坐标系
void CreateMyUtm() {
//创建空间参考对象
OGRSpatialReference oSRS;
//设置地理坐标系和投影坐标系
oSRS.SetProjCS("UTM 17 (WGS84) in northern hemisphere.");
//设置地理坐标系和投影坐标系
oSRS.SetWellKnownGeogCS("WGS84");
oSRS.SetUTM(17, TRUE);
//转成WKT并打印
char* pszWKT = NULL;
const char* apszOptions[] = { "FORMAT=WKT2_2018", "MULTILINE=YES", nullptr };
oSRS.exportToWkt(&pszWKT, apszOptions);
std::cout << pszWKT << std::endl;
CPLFree(pszWKT);
}解析空间参考信息
接下来我们主要说明如何从一个OGRSpatialReference对象中获取空间参考的各种信息，比如可以使用 OGRSpatialReference：：IsProjected（）和 OGRSpatialReference：：IsGeographic（）方法来判断该空间参考是地理坐标系还是投影坐标系。函数OGRSpatialReference：GetSemiMajor（）、OGRSpatialReference：GetSemiMinor（）和 OGRSpatialReference：GetInvFlattening（）用来获取参考椭球体的信息，分别是参考椭球体的长半轴、短半轴以及扁率的倒数。使用OGRSpatial Reference：GetAttrValue（）方法可以用来获取PROJCS、GEOGCS、DATUM、SPHEROID和PROJECTION的名称字符串。使用OGRSpatialReference：GetProjParm（）方法可以获取投影参数信息。使用OGRSpatialReference：GetLinearUnits（）方法可以获取长度单位类型并将其转换为米。
const char *pszProjection = poSRS->GetAttrValue("PROJECTION");
if( pszProjection == NULL )
{
if( poSRS->IsGeographic() )
sprintf( szProj4+strlen(szProj4), "+proj=longlat " );
else
sprintf( szProj4+strlen(szProj4), "unknown " );
}
else if( EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
{
sprintf( szProj4+strlen(szProj4),
"+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ",
poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0) );
}坐标变换
OGRSpatialReference oSourceSRS, oTargetSRS;
OGRCoordinateTransformation *poCT;
double x, y;
oSourceSRS.importFromEPSG( atoi(papszArgv[i+1]) );
oTargetSRS.importFromEPSG( atoi(papszArgv[i+2]) );
poCT = OGRCreateCoordinateTransformation( &oSourceSRS,
&oTargetSRS );
x = atof( papszArgv[i+3] );
y = atof( papszArgv[i+4] );
if( poCT == NULL || !poCT->Transform( 1, &x, &y ) )
printf( "Transformation failed.\n" );
else
{
printf( "(%f,%f) -> (%f,%f)\n",
atof( papszArgv[i+3] ),
atof( papszArgv[i+4] ),
x, y );
}OGR库
OGR体系结构
OGR中所有的类都是以OpenGIS提供的简单要素的API为蓝本、使用C++进行实现的，主要由以下七部分内容组成。
● Geometry（ogr geometry.h）
几何对象∶几何对象类（OGRGeometry等）封装了OpenGIS模型矢量数据以及一些几何操作，并且提供常用的二进制格式和文本格式之间的转换。一个几何对象包含一个空间参考（既投影）。
● Spatial Reference（ogr_spatialref.h）
空间参考∶一个OGRSpatialReference定义了投影和水准面等信息。
● Feature（ogr_feature.h）
要素∶OGRFeature类里面定义了一个完整的要素，由一个几何对象和一系列的属性组成。
● Feature Class Definition （ogr_feature.h）
要素类定义∶OGRFeatureDefn类体现了一组相关要素（通常是整个图层）的格式（字段定义的集合）。
● Layer（ogrsf_frmts.h）
图层∶OGRLayer是一个抽象基类，用来描述在一个OGRData Source中的一个图层中的所有要素。
● Data Source（ogrsf_frmts.h）
数据源∶OGRDataSource是一个抽象基类，用于表示一个文件或者数据库，里面包含一个或者多个OGRLayer对象。
● Drivers（ogrsffrmts.h）
驱动∶OGRSFDriver表示一个指定格式的转换器，打开文件返回一个 OGRDataSource 对象，所有支持的驱动使用 OGRSFDriverRegistrar 来进行管理。
其中Geometry子类比较多，分别对应不同的矢量数据类型
Geometry（几何对象）
矢量数据操作
基本操作
下面以一个shp文件为例测试如何操作一个矢量数据
数据说明：该数据为重庆市县级行政区划数据，shapefile格式，要素为多边形，
注意点：
1 gdal应该是默认utf-8编码，注意vs本身GB编码带来的问题
2 注意数据源的释放，主要是GDALDataset和OGRFeature，但是GDAL新版本结合c++11的特性有不用释放的写法
3 打开数据时有参数可以输入，详见gdal文档中GDALOpenEX函数说明，主要是第二个参数可以用逻辑或运算符输入多个，可以同时设定文件类型（GDAL_OF_VECTOR）以及打开方式（GDAL_OF_READONLY (exclusive)or GDAL_OF_UPDATE or GDAL_OF_SHARED），这关系到能否对数据进行修改。
gdal.h: Raster C API — GDAL documentation
//.h
#include <iostream>
#include "ogrsf_frmts.h"
void OpenESRIShp(const char*);
//.cpp
void OpenESRIShp(const char* FileDIR) {
//支持中文路径
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
//属性表字符编码
//CPLSetConfigOption("SHAPE_ENCODING", "");
//注册所有驱动
OGRRegisterAll();
//打开数据
//指针打开
//auto *poDS = (GDALDataset*)GDALOpenEx(FileDIR, GDAL_OF_VECTOR, NULL, NULL, NULL);
//智能指针打开
GDALDatasetUniquePtr poDS(GDALDataset::Open(FileDIR, GDAL_OF_VECTOR));
if (poDS == NULL)
{
std::cout << "Open failed" << std::endl;
exit(1);
}
//获取图层个数
int layerNum = poDS->GetLayerCount();
std::cout << layerNum << std::endl;
//打开第一个图层（Layer）
auto poLayer = poDS->GetLayer(0);
std::cout << poLayer->GetFeatureCount() << std::endl;
//打开第一个要素（feature）
//auto poFeature = poLayer->GetFeature(0);
//std::cout << poFeature->GetFID() << std::endl;
//打开属性表，输出表头
OGRFeatureDefn* poFDefn = poLayer->GetLayerDefn();
std::cout << poFDefn->GetFieldCount() << std::endl;
for (int FieldIndex = 0; FieldIndex < poFDefn->GetFieldCount(); FieldIndex++) {
auto&& layerField0 = poFDefn->GetFieldDefn(FieldIndex);
std::cout << layerField0->GetNameRef() << " ";
}
//获得每个要素的属性表
//c++ 11 写法
//万能引用，oField实际类型为左值引用迭代器 const OGRFeature::FiledValue&
//也可以不判断每个oField的类型，直接GetFieldAsString
//printf懒得改成cout了
for (const auto& poFeature : poLayer) {
for (auto&& oField : *poFeature)
{
if (oField.IsUnset())
{
printf("(unset),");
continue;
}
if (oField.IsNull())
{
printf("(null),");
continue;
}
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:
// GetString() returns a C string
printf("%s,", oField.GetString());
break;
default:
// Note: we use GetAsString() and not GetString(), since
// the later assumes the field type to be OFTString while the
// former will do a conversion from the original type to string.
printf("%s,", oField.GetAsString());
break;
}
}
std::cout << std::endl;
}
//旧写法
//OGRFeatureDefn* poFDefn = poLayer->GetLayerDefn();
//for (int iField = 0; iField < poFDefn->GetFieldCount(); iField++)
//{
// if (!poFeature->IsFieldSet(iField))
// {
// printf("(unset),");
// continue;
// }
// if (poFeature->IsFieldNull(iField))
// {
// printf("(null),");
// continue;
// }
// 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;
// }
//}
//读取几何图形（Geometry）
//读取图层第一个要素（feature）
auto poFeature = poLayer->GetFeature(0);
OGRGeometry* poGeometry;
poGeometry = poFeature->GetGeometryRef();
//如果为点
if (poGeometry != NULL
&& wkbFlatten(poGeometry->getGeometryType()) == wkbPoint)
{
OGRPoint* poPoint = poGeometry->toPoint();
printf("%.3f,%3.f\n", poPoint->getX(), poPoint->getY());
}
//如果为多边形，则输出多边形所有点
else if(poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon)
{
OGRPolygon* poPolygon = poGeometry->toPolygon();
//看看多边形有几条边，一般都是一条边从头到底，即外环（有空洞的情况不能这么写）
int CurveNum = 0;
for (auto poCurve : *poPolygon) {
CurveNum++;
for (auto poPoint : *poCurve) {
std::cout << "x=" << poPoint.getX() << " ,";
std::cout << "y=" << poPoint.getY() << std::endl;
}
std::cout << CurveNum << std::endl;
}
}
else {
std::cout << poGeometry->getGeometryName() << std::endl;
std::cout << poGeometry->getGeometryType() << std::endl;
std::cout << "几何图案读取失败" << std::endl;
}
//保险起见，销毁对象，gdal的文档搞不明白新版本这个对象要不要释放
OGRFeature::DestroyFeature(poFeature);
//gdal3 c++ 11 不需要释放资源
//GDALClose(poDS);
}
//参考代码
// gdal2.3+ c++11
//{
//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)
// {
// if (oField.IsUnset())
// {
// printf("(unset),");
// continue;
// }
// if (oField.IsNull())
// {
// printf("(null),");
// continue;
// }
// 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:
// // GetString() returns a C string
// printf("%s,", oField.GetString());
// break;
// default:
// // Note: we use GetAsString() and not GetString(), since
// // the later assumes the field type to be OFTString while the
// // former will do a conversion from the original type to string.
// 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;
//}
//旧版本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++)
// {
// if (!poFeature->IsFieldSet(iField))
// {
// printf("(unset),");
// continue;
// }
// if (poFeature->IsFieldNull(iField))
// {
// printf("(null),");
// continue;
// }
// 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);
//}特别说明：编码问题
vs2019的编码问题非常头疼，首先vs和win系统都是gbk编码，尤其是win系统的gbk编码，意味着如果你强行设置全局utf-8编码后会导致根本读取不了中文文件，所以中文路径读取的时候一定是gbk编码
//设置文件名不是utf-8编码
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");其次gdal内部完全是utf-8编码的，这点我们也改变不了，但由于win系统是gbk编码，因此cmd控制台也只支持gbk，权衡之下我们只能放弃控制台，接受控制台调试的时候看到乱码，才能保证涉及到中文功能的正常运行。我们只要保证最后输出到文件的结果是正常的就行了。
另外c++ 11之后支持string强制utf-8编码，凭借这点我们可以对属性表的字符串进行比较了，后面实例时详细说。
string s1 = u8"市辖区";在网上查阅了相当多的资料，无非是以下几种：
1 源码文件保存为utf-8
2 设置全局环境utf-8
system("chcp 65001");3 通过命令行设置输入为utf-8，执行时为gbk（项目->属性->c/c++->命令行）
/source-charset:utf-8 /execution-charset:gbk 还没完全搞明白到底怎么最合适，以后再更
坐标系相关
待更新
准备单独写一篇文章将栅格和矢量数据的坐标系相关操作统一归纳。
属性相关
基本操作中已经交代了属性表的读取方法，下面记录一些编辑操作。新建一个类用静态成员函数去实现各类功能。首先注意的是文件打开的属性，上文提到了GDAL_OF_READONLY (exclusive)or GDAL_OF_UPDATE or GDAL_OF_SHARED这三个选项，如果不设置的话默认是readonly，再修改属性表时就会报错。
//CRS_defn.h
#pragma once
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include <iostream>
class CRS_defn {
public:
//添加属性列
static void addDefn(OGRLayer* , OGRFieldDefn*);
//给属性表某列赋值
template <typename F>
static void assignDefn(OGRLayer* , int& , F&&);
//删除指定属性表列
static void deleteDefn(OGRLayer*, int);
};属性表中插入列
api就是OGRLayer::createField，参数输入OGRFieldDefn*。OGRFieldDefn类中设置属性表列的名称和数据类型，这里列名设置成testDefn，数据类型为int。
注意代码规范，设置有效性检查。
//CRS_defn.cpp
void CRS_defn::addDefn(OGRLayer* tarLayer , OGRFieldDefn* newField) {
if (tarLayer == NULL || newField == NULL) {
std::cout << "输入无效" << std::endl;
return;
}
if (tarLayer->CreateField(newField) == OGRERR_NONE) {
std::cout << "新增属性列成功" << std::endl;
}
}
//int main()
//全局utf-8编码
system("chcp 65001");
// 增加属性列
GDALAllRegister();
const char* chongqingDIR = "C:/Users/23005/Desktop/chongqing/chongqing.shp";
GDALDatasetUniquePtr poDS(GDALDataset::Open(chongqingDIR, GDAL_OF_VECTOR || GDAL_OF_UPDATE));
if (poDS == NULL)
{
printf("Open failed.\n");
exit(1);
}
OGRLayer* poLayer = poDS->GetLayer(0);
OGRFieldDefn* testFieldDefn = new OGRFieldDefn("testDefn",OFTInteger);
CRS_defn::addDefn(poLayer,testFieldDefn);
delete testFieldDefn;属性列删除
属性列删除比较简单，调用OGRLayer::deleteDefn即可
//CRS_defn.cpp
void CRS_defn::deleteDefn(OGRLayer* tarLayer, int fieldIndex)
{
if (tarLayer == NULL || fieldIndex < 0 || fieldIndex >= tarLayer->GetLayerDefn()->GetFieldCount()) {
std::cout << "输入无效" << std::endl;
return;
}
tarLayer->DeleteField(fieldIndex);
}
//int main()
//获得属性表最后一列的序号
int lastFieldIndex = poLayer->GetLayerDefn()->GetFieldCount() - 1;
//删除
CRS_defn::deleteDefn(poLayer, lastFieldIndex);属性表值修改
属性表值修改复杂一点，需要对图层（Layer）循环，对图层中每个要素（Feature）进行赋值操作，api是OGRFeature::SetField()，这个函数有很多重载，对应不同数据类型，第一个参数是需要修改属性值的列号，第二个是写入的值，根据写入值的类型自动匹配相应函数，可以与读取属性值的函数相比较，是一一对应的，不过读取值的话就需要指定读取值的类型，因此有数个函数对应不同类型。
我们首先尝试给新增列简单赋值一个数字，结果运行后没有效果，查了很久的文档，得知修改feature的数据后需要使用setFeature来更新数据，同样的，新建一个矢量数据在设置完相应数据后也需要用createFeature来激活。只是对feature层面的修改需要调用setFeature，像上面直接在Layer层面增减属性列是不用的。
但是OGRLayer::setFeature这个函数的输入是OGRFeature*，而我们用的新版本的写法遍历，迭代器得到的是OGRFeature的Unique_Ptr智能指针，因此还需要用get函数获得里面的裸指针传入。
还有个api是OGRLayer::SyncToDisk()，是写入磁盘，但我测试不需要这个函数也能对shp文件进行修改了。
void CRS_defn::assignDefn(OGRLayer* tarLayer, int& fieldIndex)
{
if (tarLayer == NULL || fieldIndex < 0 || fieldIndex >= tarLayer->GetLayerDefn()->GetFieldCount()) {
std::cout << "输入无效" << std::endl;
return;
}
for (const auto& oFeature : *tarLayer) {
oFeature->SetField(fieldIndex, 1);
tarLayer->SetFeature(oFeature.get());
}
}下面再深入的进行一些修改，显然在实际应用中不可能就这么简单的给属性赋一个常值，肯定是要将需要计算的结果进行赋值。上面说过给属性表赋值需要遍历图层中的要素，实际上大部分情况下我们也需要对每个要素进行单独计算再赋值。比如说我们现在需要对示例数据中类型这一项进行一个编码处理，类型中一共有市辖区，县，自治县三种，我们分别以0，1，2来映射并写入testDefn。
我们一步一步来，首先最简单的写法开始。最基础的写法就是在循环内部进行映射处理，示意如下
if (tarLayer == NULL || fieldIndex < 0 || fieldIndex >= tarLayer->GetLayerDefn()->GetFieldCount()) {
std::cout << "输入无效" << std::endl;
return;
}
for (const auto& oFeature : *tarLayer) {
//在这里进行映射计算,分三类对result赋值
int result;
switch(){...}
oFeature->SetField(fieldIndex, result);
tarLayer->SetFeature(oFeature.get());
}
}这种写法没有问题，但问题在于如果我们需要其他赋值策略的时候，每次都要新写一个这样的函数，事实上我们每次需要写的只是循环中间计算的那块，上下都是重复的，用软件工程的话来说就是不符合低内聚高耦合。
实际上我们可以把赋值策略写成一个函数，然后作为函数指针传进去，在循环时调用就行了。我们当然可以把函数指针作为参数来传入我们写的函数，但是又有一个问题是函数指针的返回类型是要给定的，比如我们当前这个要做的这个映射，返回的是0，1，2，int类型，函数指针就得是int*，且形参类型和个数也要确定，这还是不符合我们的需求。
//函数指针
Int (*p)(int x, int y);//*p的括号很重要
Int max(int x , int y){
return (x>y)? x : y;
};
p=max;
p(3,4);
//将指向函数的指针变量作为函数参数
int wm ( int x, int y, (*midfunc)(int,int)){
Int result=midfunc(x,y);
Return result;
};
Wm(3,5,max);//函数名本身就是指针
Wm(3,5,p);//也行
这里就需要用到模板来进行类型推导以便让我们传入一个不受类型限制的函数指针，这样我们就可以实现一个通用的效果。
//注意模板成员函数的实现必须写在头文件里
//CRS_defn.h
template <typename F>
void CRS_defn::assignDefn(OGRLayer* tarLayer, int& fieldIndex, F assignRule)
{
if (tarLayer == NULL || fieldIndex < 0 || fieldIndex >= tarLayer->GetLayerDefn()->GetFieldCount() || assignRule == NULL) {
std::cout << "输入无效" << std::endl;
return;
}
for (const auto& oFeature : *tarLayer) {
oFeature->SetField(fieldIndex, assignRule());
tarLayer->SetFeature(oFeature.get());
}
//tarLayer->SyncToDisk();
}
//先传入一个简单的函数
int func1() {
return 1;
}
//主函数中进行调用
//修改属性列的值
int lastFieldIndex = poLayer->GetLayerDefn()->GetFieldCount() - 1;
CRS_defn::assignDefn(poLayer, lastFieldIndex, func1);如果功能比较简单，像这种赋以一个常值的需求，我们还可以用传入lambda函数来实现。
//主函数中
int lastFieldIndex = poLayer->GetLayerDefn()->GetFieldCount() - 1;
CRS_defn::assignDefn(poLayer, lastFieldIndex, [] {return 1; });接下来我们来实现这个映射功能。本来想用static enum加switch来写的，然而不会，算了，直接if判断好了。
这边就涉及到了上面说的编码问题，多次尝试后发现这样写是ok的
//CRS_defn.h
template <typename F>
void CRS_defn::assignDefn(OGRLayer* tarLayer, int& fieldIndex, F assignRule)
{
if (tarLayer == NULL || fieldIndex < 0 || fieldIndex >= tarLayer->GetLayerDefn()->GetFieldCount() || assignRule == NULL) {
std::cout << "输入无效" << std::endl;
return;
}
for (const auto& oFeature : *tarLayer) {
oFeature->SetField(fieldIndex, assignRule(oFeature.get()));
tarLayer->SetFeature(oFeature.get());
}
//tarLayer->SyncToDisk();
}
//计算函数
int func1(OGRFeature* poFeature) {
string xianType = poFeature->GetFieldAsString(6);
if (xianType == u8"市辖区")
return 0;
else if (xianType == u8"县")
return 1;
else if (xianType == u8"自治县")
return 2;
else
return -1;
}
//主函数
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
const char* chongqingDIR = "C:/Users/23005/Desktop/chongqing/重庆.shp";
GDALDatasetUniquePtr poDS(GDALDataset::Open(chongqingDIR, GDAL_OF_VECTOR || GDAL_OF_UPDATE));
if (poDS == NULL)
{
printf("Open failed.\n");
exit(1);
}
CPLSetConfigOption("SHAPE_ENCODING", "");
OGRLayer* poLayer = poDS->GetLayer(0);
//修改属性列的值
int lastFieldIndex = poLayer->GetLayerDefn()->GetFieldCount() - 1;
CRS_defn::assignDefn(poLayer, lastFieldIndex, func1);下面测试一下能否将中文正常写入属性表，果然不加u8前缀结果就是乱码
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
const char* chongqingDIR = "C:/Users/23005/Desktop/chongqing/重庆.shp";
GDALDatasetUniquePtr poDS(GDALDataset::Open(chongqingDIR, GDAL_OF_VECTOR || GDAL_OF_UPDATE));
if (poDS == NULL)
{
printf("Open failed.\n");
exit(1);
}
OGRLayer* poLayer = poDS->GetLayer(0);
//新建string类型属性表列并插入
OGRFieldDefn* testFieldDefn = new OGRFieldDefn("testStringDefn", OFTString);
CRS_defn::addDefn(poLayer, testFieldDefn);
delete testFieldDefn;
//修改属性列的值，注意utf-8编码
int lastFieldIndex = poLayer->GetLayerDefn()->GetFieldCount() - 1;
CRS_defn::assignDefn(poLayer, lastFieldIndex, [](OGRFeature* poFeature) {return u8"你好1"; });实际开发时这个函数的返回值应该是一个信号，用来传递是否成功修改，而不是现在的void
要素相关
新建CRS_feature类实现要素相关操作
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include <iostream>
class CRS_feature {
public:
static void addFeature(OGRLayer*);
static void deleteFeature(OGRLayer*, long);
static void moveFeature(OGRLayer*, long, const std::vector<float>&);
};删除要素
删除要素很简单，调用api就一句话
//根据FID删除指定要素
void CRS_feature::deleteFeature(OGRLayer* tarLayer, long FID) {
tarLayer->DeleteFeature(FID);
}修改要素
修改要素中集合体的坐标。一开始用auto循环发现并不能对源数据生效，只能通过指针来操作。
这里简单说明一下polygon的格式。OGRPolygon大概是我们最常碰到的要素类型，上面矢量基本操作的代码对polygon的操作太简单粗暴了。实际上polygon可能由于存在空洞问题而有数条环（OGRLinearRing），如果polygon是简单多边形，那就只有一条外环。至于setPoint这个函数应该还是复制了传入点的坐标，不用担心相关释放问题。
//平移要素（简单多边形）
void CRS_feature::moveFeature(OGRLayer* tarLayer, long FID, const std::vector<float>& offset) {
OGRFeature* tarFeature = tarLayer->GetFeature(FID);
OGRGeometry* tarGeometry = tarFeature->GetGeometryRef();
OGRPolygon* poPolygon = tarGeometry->toPolygon();
OGRLinearRing* tarRing = poPolygon->getExteriorRing();
for (int i = 0; i < tarRing->getNumPoints(); i++) {
OGRPoint* originPoint = new OGRPoint();
tarRing->getPoint(i, originPoint);
originPoint->setX(originPoint->getX() + offset[0]);
originPoint->setY(originPoint->getY() + offset[1]);
tarRing->setPoint(i, originPoint);
delete originPoint;
}
tarLayer->SetFeature(tarFeature);
tarLayer->SyncToDisk();
OGRFeature::DestroyFeature(tarFeature);
}插入要素
新增要素，需要准备好新要素的集合体和属性，绑定后用图层加入。
//新增要素
void CRS_feature::addFeature(OGRLayer* tarLayer) {
//拿图层中第一个要素（万州区）平移一下，作为新要素的几何体
OGRFeature* copyFeature = tarLayer->GetFeature(0);
OGRPolygon* newPolygon = (OGRPolygon*)OGRGeometryFactory::createGeometry(wkbPolygon);
OGRLinearRing* newLinearRing = (OGRLinearRing*)OGRGeometryFactory::createGeometry(wkbLinearRing);
OGRPoint newPoint;
auto copyGeometry = copyFeature->GetGeometryRef();
auto copyPolygon = copyGeometry->toPolygon();
auto copyRing = copyPolygon->getExteriorRing();
for (auto copyPoint : copyRing) {
newPoint.setX(copyPoint.getX() + 1.0);
newPoint.setY(copyPoint.getY() + 1.0);
newLinearRing->addPoint(&newPoint);
}
//将环封闭
newLinearRing->closeRings();
//将环赋给多边形作为外环
newPolygon->addRing(newLinearRing);
//根据原图层属性新建一个元素
OGRFeature* newFeature = OGRFeature::CreateFeature(tarLayer->GetLayerDefn());
//填入属性（填个name，其他懒得填）
newFeature->SetField(1, u8"万州区平移");
//绑定几何体
newFeature->SetGeometry(newPolygon);
//将要素加入图层
tarLayer->CreateFeature(newFeature);
//释放相关指针
OGRFeature::DestroyFeature(newFeature);
std::cout << "插入要素成功" << std::endl;
}选择要素
按空间范围选择要素，可以直接使用SetSpatialFilter设置空间过滤器，也可以用SetSpatialFilterRect传入一个矩形范围。
//按四至范围的四分之一（左下角）选择元素
void CRS_feature::selectBySpatialFilter(OGRLayer* tarLayer) {
OGREnvelope tarLayerExtent;
tarLayer->GetExtent(&tarLayerExtent);
tarLayerExtent.MaxX = tarLayerExtent.MaxX - (tarLayerExtent.MaxX - tarLayerExtent.MinX) / 2;
tarLayerExtent.MaxY = tarLayerExtent.MaxY - (tarLayerExtent.MaxY - tarLayerExtent.MinY) / 2;
tarLayer->SetSpatialFilterRect(tarLayerExtent.MinX, tarLayerExtent.MinY, tarLayerExtent.MaxX, tarLayerExtent.MaxY);
//重置选择
//tarLayer->SetSpatialFilter(NULL);
//读取所有要素的所有属性
for (const auto& poFeature : tarLayer) {
for (auto&& oField : *poFeature)
{
if (oField.IsUnset())
{
printf("(unset),");
continue;
}
if (oField.IsNull())
{
printf("(null),");
continue;
}
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:
// GetString() returns a C string
printf("%s,", oField.GetString());
break;
default:
// Note: we use GetAsString() and not GetString(), since
// the later assumes the field type to be OFTString while the
// former will do a conversion from the original type to string.
printf("%s,", oField.GetAsString());
break;
}
}
std::cout << std::endl;
}
}按属性选择见下一章SQL说明
SQL说明
使用SQL进行属性查询/筛选
OGR SQL dialect — GDAL documentation
SQL查询可以从数据源开始（完整SQL语句），也可以从图层查询（WHERE子句）
注意GDALDatasetUniquePtr传入一定是引用，另外中文的话要用转义符\前后加上引号，并且结果要回收。
//按属性选择要素(输入数据源)
void CRS_feature::selectByAttribute(GDALDatasetUniquePtr& tarDS) {
string tarLayerName = tarDS->GetLayer(0)->GetName();
string SQLCommand = "SELECT * FROM \"" + tarLayerName + "\" WHERE testDefn = 2";
OGRLayer* tarLayer = tarDS->ExecuteSQL(SQLCommand.c_str(), NULL, NULL);
//读取所有要素的所有属性
for (const auto& poFeature : tarLayer) {
for (auto&& oField : *poFeature)
{
if (oField.IsUnset())
{
printf("(unset),");
continue;
}
if (oField.IsNull())
{
printf("(null),");
continue;
}
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:
// GetString() returns a C string
printf("%s,", oField.GetString());
break;
default:
// Note: we use GetAsString() and not GetString(), since
// the later assumes the field type to be OFTString while the
// former will do a conversion from the original type to string.
printf("%s,", oField.GetAsString());
break;
}
}
std::cout << std::endl;
}
tarDS->ReleaseResultSet(tarLayer);
}
//按属性选择要素(输入数据源)
void CRS_feature::selectByAttribute(OGRLayer* tarLayer) {
string SQLCommand = "testDefn = 2";
tarLayer->SetAttributeFilter(SQLCommand.c_str());
//读取所有要素的所有属性
for (const auto& poFeature : tarLayer) {
for (auto&& oField : *poFeature)
{
if (oField.IsUnset())
{
printf("(unset),");
continue;
}
if (oField.IsNull())
{
printf("(null),");
continue;
}
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:
// GetString() returns a C string
printf("%s,", oField.GetString());
break;
default:
// Note: we use GetAsString() and not GetString(), since
// the later assumes the field type to be OFTString while the
// former will do a conversion from the original type to string.
printf("%s,", oField.GetAsString());
break;
}
}
std::cout << std::endl;
}
}OGRSQL的内容还不少，等用到了再来更新。
矢量数据保存
待更新
Comments 3 条评论
博主 匿名
gdal，在conda环境中，直接conda install gdal，速度还是可以的，maybe是因为我换源了？
博主 匿名
@匿名 哈哈这我就不清楚了，很久不搞python了，可能现在新版本更好安装了吧
博主 匿名
@匿名 应该是我换了国内镜像的原因