python GDAL开发源码转c++源码（tif文件的解析）

公司考虑python方面的开发人手偏少，需要后续代码维护，另一方面c++运行效率较高，被分配了一个python源码转c++源码的活，给了一份python源码文件，经过大约两周多的时间，主体功能基本转化完成。

python源码主要功能是通过GDAL的相关API对tif文件进行文件解析，生成相应目标文件方便其他程序使用。

这里以个人理解记录下开发过程中python和c++调用GDAL之间的库函数对应关系，以及开发和调试时遇见的问题，中间有表述错误或不当的欢迎交流，一起学习。

GDAL相关的帮助文档及环境下载

开发帮助文档：

GDAL — GDAL 文档 (osgeo.cn)

三方库下载：

GISInternals Support Site

GDALAPI中python和c++使用时对照关系

0.简单的函数对照

说明	python	c++
导入文件	from osgeo import gdal from osgeo import osr	#include "gdal.h" #include "ogr_spatialref.h"  #include "gdalwarper.h"  #include "gdal_priv.h"
gdal的数据初始化注册操作	gdal.AllRegister()	GDALAllRegister();
gdal驱动获取，参数可以选择其他类型	gdal.GetDriverByName( 'MEM' )	(GDALDriver*)GDALGetDriverByName("MEM");
根据输入文件获取gdal数据集	gdal.Open(self.input, gdal.GA_ReadOnly)	(GDALDataset*)GDALOpen(m_input.c_str(), GA_ReadOnly);
一个参考空间坐标系的声明	self.in_srs = osr.SpatialReference()	OGRSpatialReference in_srs;
获得wkt坐标系	self.in_srs = osr.SpatialReference() in_srs_wkt=self.in_srs.ExportToWkt()	OGRSpatialReference in_srs ； char *in_srs_wkt=NULL; in_srs .exportToWkt(&in_srs_wkt);
获得wkt坐标系	self.in_ds = gdal.Open(self.input, gdal.GA_ReadOnly) self.in_srs_wkt = self.in_ds.GetProjection()	GDALDataset* in_ds= (GDALDataset*)GDALOpen(m_input.c_str(), GA_ReadOnly); in_srs_wkt = in_ds->GetProjectionRef();
获得 proj4格式 坐标系	one_sproj4=in_srs.ExportToProj4()	OGRSpatialReference in_srs char* one_sproj4; in_srs.exportToProj4(&one_sproj4);
	OGRSpatialReference in_srs，srs4326； osr.CoordinateTransformation(self.in_srs, srs4326)	OGRSpatialReference in_srs，srs4326； OGRCreateCoordinateTransformation(&in_srs,&srs4326);
打印显示进度条函数	gdal.TermProgress_nocb(complete)	GDALTermProgress(complete,NULL,NULL)

1.GetNoDataValue。获取gdal数据集的nodata值

也不太清除nodata是个啥东西，就这样表述了。

self.in_ds = gdal.Open(self.input, gdal.GA_ReadOnly)
for i in range(1, self.in_ds.RasterCount+1):
  if self.in_ds.GetRasterBand(i).GetNoDataValue() != None:
    #直接添加nodata列表
    self.in_nodata.append( self.in_ds.GetRasterBand(i).GetNoDataValue() )

GDALDataset* poDSRef = (GDALDataset*)GDALOpen(m_input.c_str(), GA_ReadOnly);
int nCount = poDSRef->GetRasterCount();
vector<double> m_vNoData;
for (int i = 1; i < nCount+1; ++i)
{
	int bSuccess(0); 
	double dNoData;
	if (poDSRef->GetRasterBand(i) != NULL) {
        //获取nodata值
		dNoData = poDSRef->GetRasterBand(i)->GetNoDataValue(&bSuccess);
	}
	if (bSuccess)
	{//当存在nodata值时
		m_vNoData.push_back(dNoData);
	}
}

2.ImportFromEPSG。导入一个坐标系

self.out_srs = osr.SpatialReference()
#两种坐标系，墨卡托、大地
if self.options.profile == 'mercator':
   self.out_srs.ImportFromEPSG(900913)
elif self.options.profile == 'geodetic':
   self.out_srs.ImportFromEPSG(4326)

/*
需要预先设置空间坐标系数据值的所在位置，否则下面的importFromEPSG函数会报错找不到，可以通过函数返回值进行查看是否获取成功。
该数据文件会跟gdal库文件一同下载
*/
OGRSpatialReference m_out_srs；
CPLSetConfigOption("GDAL_DATA", "F:\\gdal2.3.2_vs2019\\data");
if (m_options.m_profile == "mercator") {
	int test=m_out_srs.importFromEPSG(900913);//900913 3857
}
else if(m_options.m_profile == "geodetic"){
	int test = m_out_srs.importFromEPSG(4326);
}

所需数据文件集所在位置：

3.GetGeoTransform。

GDALDataset in_ds
if (in_ds.GetGeoTransform() == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) ：
    pass

double trsf[6] = { 0 };
GDALDataset*poDSRef
poDSRef->GetGeoTransform(trsf); 
if (trsf[0] == 0.0 && trsf[1] == 1.0 
    && trsf[2] == 0.0 && trsf[3] == 0.0 
    && trsf[4] == 0.0 && trsf[5] == 1.0 )
{
	printf("test);
}

4.AutoCreateWarpedVRT.

GDALDatasetin_srs;
const char* in_srs_wkt;
OGRSpatialReference out_srs;
GDALDataset self.out_ds = gdal.AutoCreateWarpedVRT( self.in_ds, self.in_srs_wkt, self.out_srs.ExportToWkt() )
                    

/*
需要声明下proj.db的文件路径，要不GDALAutoCreateWarpedVRT函数报错，返回值为NULL
*/
const char* newprojs[] = { "F:\\release-1928-x64-gdal-3-2-2-dll\\bin\\proj7\\share",NULL };
OSRSetPROJSearchPaths(newprojs);

GDALDataset * in_ds;
const char* in_srs_wkt;
const char* out_srs_wkt;
GDALDataset *m_out_ds = (GDALDataset*)GDALAutoCreateWarpedVRT(in_ds, in_srs_wkt, out_srs_wkt, GRA_NearestNeighbour , 0,NULL);

5.CreateCopy复制数据集

#把数据集out_ds中数据复制到tempfilename文件中
GDALDataset out_ds
self.out_ds.GetDriver().CreateCopy(tempfilename, self.out_ds)

GDALDataset *m_out_ds；
//将m_out_ds数据复制到文件original_backup_file
GDALDataset* copyTmp= m_out_ds->GetDriver()->CreateCopy(original_backup_file.c_str(), m_out_ds,false,NULL,0,NULL);
//数据集复制完成后，函数返回指针会占用生成的数据文件（delete后数据好像才正式写到目标文件中），对生成文件的读写操作前，需要释放该指针
delete copyTmp;

6.读写操作。python(ReadRaster、WriteRaster)，c++(RasterIO)

python	c++
ReadRaster(..)	RasterIO(GF_Read,...)
WriteRaster(...)	RasterIO(GF_Write,...)

data = ds.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize, band_list=list(range(1,self.dataBandsCount+1)))



alpha = self.alphaband.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize)

int m_dataBandsCount;
int* p_band_list = (int*)malloc(sizeof(int) * m_dataBandsCount);
for (int i = 1; i < m_dataBandsCount + 1; ++i) {
	*(p_band_list + i - 1) = i;
}
double* outDsRasterIO = (double*)CPLMalloc(sizeof(double)* wxsize * wysize);
GDALDataset *m_out_ds;
if (m_out_ds != NULL) {
/*
功能：读m_out_ds 数据并写到缓冲区outDsRasterIO
这个函数与python代码比较，python代码调用时不定的少很多参数，根据使用场景进行修改定吧。
描述：
（1）GF_Read：决定读还是写
（2）rx, ry, rxsize, rysize：读数据的位置信息
（3）outDsRasterIO：数据读到的缓冲区
（4）wxsize, wysize：数据缓冲区的大小，一般与outDsRasterIO开辟的大小有关，正常来说outDsRasterIO定了这个参数就跟着定了。python有地方用的None，需要相应修改。
（5）m_tiledata：数据操作的格式，如int还是float。这里python中有时不写，需要自己根据情况决定数据操作格式。像我这里工程整体基本用的int16，我就用它来做默认值了。
（6）m_dataBandsCount, p_band_list：需要读的波段数据的数量以及其指针列表。
（7） 0, 0, 0：最后三个参数，按默认值0传参，需要的话可以研究下去改变参数
*/
	m_out_ds->RasterIO(GF_Read, rx, ry, rxsize, rysize, outDsRasterIO, wxsize, wysize, m_tiledata, m_dataBandsCount, p_band_list, 0, 0, 0);
}

//读缓冲区outDsRasterIO数据并写到m_out_ds
//这里的参数调用与上面的类似处理就行
m_out_ds->RasterIO(GF_Write, wx, wy, wxsize, wysize, outDsRasterIO, wxsize, wysize, m_tiledata, m_dataBandsCount, p_band_list, 0, 0, 0);


double* outDsAlphabandIO = NULL;
outDsAlphabandIO = (double*)CPLMalloc(sizeof(double) * wxsize * wysize);
GDALRasterBand *m_alphaband = m_out_ds->GetRasterBand(1)->GetMaskBand();
if (m_alphaband != NULL) {
//这里的参数调用与上面的类似处理就行
	m_alphaband->RasterIO(GF_Read, rx, ry, rxsize, rysize, outDsAlphabandIO, wxsize, wysize, m_alphaband->GetRasterDataType(), 0, 0, 0);
}

 +++++++++++++++++++++++++++++++++++++++

20231010补充

项目的后续数据测试过程中出现问题了，得稍微深入了解下RasterIO的函数作用了，由于以前没做过相关的开发，理解起来有些困难，粘个大佬的文章链接，方便随时学习。深入解析GDAL库的RasterIO()函数-CSDN博客

 ++++++++++++++++++++++++++++++++++++++++

20231025补充

关于参数tileData（数据的格式），又碰见问题了，在读灰度图像文件的时候（其格式为float类型），一直和python的切片结果不一样，python中切出来的图片能够看到层次的灰色，但自己转的c代码不行看不到，转的源代码如下：

//开辟数据缓冲内存
double* outDsRasterIO = (double*)malloc(sizeof(double)* wxsize * wysize* m_dataBandsCount);
//进行数据读取
//m_tileData=Byte
//m_out_ds数据类型为float
int test= m_out_ds->RasterIO(GF_Read, rx, ry, rxsize, rysize, outDsRasterIO, wxsize, wysize, m_tileData, m_dataBandsCount, p_band_list,0, 0, 0);
					

 最初理解，以double类型开辟空间，空间够大什么类型的数据都能够支持，省的自己再根据每个数据集类型进行不同的初始化。问题就出在m_tileData这个参数这里，一直获取到错误数据，后面经过测试以及又仔细看了下上面链接中对该参数的描述，特别当读数据时，该参数类型要与开辟的存储空间数据类型一致，这里就是tileData与outDsRasterIO类型不一致导致数据异常。然后做了以下尝试，数据都正确。

//1.开辟float类型空间并使用float类型调用IO
float* outDsRasterIO = (float*)malloc(sizeof(float)* wxsize * wysize* m_dataBandsCount);
int test= m_out_ds->RasterIO(GF_Read, rx, ry, rxsize, rysize, outDsRasterIO, wxsize, wysize, GDT_Float32, m_dataBandsCount, p_band_list,0, 0, 0);
//2.开辟double类型空间并使用double类型调用IO
double* outDsRasterIO = (double*)malloc(sizeof(double)* wxsize * wysize* m_dataBandsCount);
int test= m_out_ds->RasterIO(GF_Read, rx, ry, rxsize, rysize, outDsRasterIO, wxsize, wysize, GDT_Float64, m_dataBandsCount, p_band_list,0, 0, 0);
					

 不过根据描述，当m_out_ds数据类型为double时，然后用float初始化空间可能也会产生异常，自己没做测试。

20231101补充

在跑数据进行测试的时候，数据比较大，28g，在基于墨卡托、大地两种坐标系格式下，读取文件的时候很慢，断点测试的时候，是RasterIO函数每次读写操作影响的速度。然后就想开辟线程以为多几条路一起跑代码就会快很多，但结果是更慢了，蒙圈.......

然后在网上搜索有没有前辈们走过gdal多线程切片的路子，看到了下面的一些知识，又学废了，粘过来记录一下，跟以往理解的开辟多线程都能提高效率有点出入，结论就是受制与硬盘的IO速度，使用多线程加速RasterIO的速度是不可行的：

（1） 读写最好还是不要多线程，硬盘读写的速度有限，单线程时已经满负荷了，多线程又会增加线程之间的切换，会增加时间。如果想增加读写速度，应该增加硬盘，做raid

（2）首先是硬盘的写入是串行的，CPU的计算才是并行的，如果你偏重计算那么多线程能提高，要不怎么叫做并行计算呢；如果侧重存储，除非数据量达到足以体现优势的程度，否则加上线程之间切换的损耗当然会效率更加地下。

（3）这个是按照算法来说的，目前来说大多数的算法都是很快的，瓶颈都在磁盘的IO上，我们针对大多数的算法都进行过测试，基本一半以上的时间都耗费在磁盘的IO上。比如我处理一个影像，处理数据用了1分钟，写入图像用了2分钟，那你把你的算法优化的很牛逼，10秒中搞定，你的效率提高了多少，但是如果我多线程写入的话，我效率提高一倍，也就是写入图像用了1分钟，那这个效率明显比你优化你的算法来的实惠。这个东西还是要针对算法来说的。

（4）磁盘IO单线程顺序写时最快的，如果多线程写，磁盘的磁头要不断重新寻址，所以写入速度反而会慢。

python函数与c++的对应

1.python中的fromstring

/*
这里实现了跟python代码一样的效果，但刚刚写这个文章的时候又看不懂了，但目前数据测试没有问题的。!!后面数据测试过程中有问题再改吧。
*/
/*
因为目前只做了int16的数据转化，一下以int16举例
这是python的函数
fromstring(data, int16)//data=char[2]
if a >= -1000 :
   a = (a + 1000) * 5
else:
   a = a - a 
   a = (a + 1000) * 5
dst.write(a)
*/
data=char[2];
unsigned char inTmp[2] = { data[0],data[1] };
INT16 tmp=(inTmp[1] << 8) | inTmp[0];
if (tmp >= -1000) {
	tmp = (tmp + 1000) * 5;
}
else {
	tmp -= tmp;
	tmp = (tmp + 1000) * 5;
}
unsigned char ret[2];
ret[1] = tmp >> 8;
ret[0] = tmp;
dst.write(ret, sizeof(char) * 2);

还有就是上面流程转化完之后，过程中是有产生一些问题，到现在也还不太理解问题所在（尽管现在先曲线救国了），在我另外提的问答里，c++文件二进制读写过程，INT16转char操作后产生的数据异常_编程语言-CSDN问答

这里有大佬解释比较详细：

python 语言， 详解fromstring 函数 -------------从字符串到Ascii 码的转换-CSDN博客

过程中遇见c++的一些库函数使用

0.简单的c++函数

说明	函数
依据输入的创建文件路径，创建文件夹	#include <windows.h> CreateDirectory(stringToLPCWSTR(m_output), NULL)

1.tmpnam.临时文件的建立和删除

/*
代码运行时提示tmpnam函数是不安全的，此定义可以屏蔽此提示。
vs建议使用tmpnam_s,但使用过程遇见了些问题，也就没用了。
*/
#define _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
//生成临时文件，返回生成的临时文件名
std::string original_backup_file = std::tmpnam(NULL);
//删除临时文件
_unlink(fileName.c_str());

2.读文件所有内容

#include<fstream>
//读一个文件的所有内容，fileName文件名；readType  ios::in | ios::binary等文件方式
std::string self_readFileContent(const std::string& fileName,  ios_base::openmode readType=ios::in) {
	std::ifstream fInput(fileName, readType);
	std::stringstream buffer;
	buffer << fInput.rdbuf();
	std::string fContent(buffer.str());
	fInput.close();
	return fContent;
}

3.boost::replace_all_copy，字符串替换

#include <boost/algorithm/string.hpp>

//将文本fContent中的lowStr替换为newStr，fContent被修改
boost::replace_all(fContent, lowStr, newStr);
//将文本fContent中的lowStr替换为newStr，返回新的字符串，fContent不被修改
string newStr=boost::replace_all_copy(fContent, lowStr, newStr);

4.元组的使用

#include <tuple>
tuple<double, double, double, double> m_swne 
    = tuple<double, double, double, double>(0，0，0，0);
//数据使用
get<0>(m_swne), get<1>(m_swne), get<2>(m_swne), get<3>(m_swne)

5.double转string并设置精度

#include<sstream>
#include <iomanip>
//转化数据，保留精度
string doubleToString(double val, int fixed) {
	std::stringstream  oss;
	oss << setprecision(fixed) << val;
	return oss.str();
}

6.判断文件是否存在

#include <sys/stat.h>
//判断字符串表示的文件或文件夹路径是否是否存在.0:不存在，1：存在文件夹，2：存在文件
int judgeFilePath(const string filePath) {
	int res = 0;
	struct stat s;
	if (stat(filePath.c_str(), &s) == 0) {
		if (s.st_mode & S_IFDIR) {
			res = 1;
			//std::cout << "it's a directory" << std::endl;
		}
		else if (s.st_mode & S_IFREG) {
			res = 2;
			//std::cout << "it's a file" << std::endl;
		}
		else {
			//std::cout << "not file not directory" << std::endl;
		}
	}
	else {
		//std::cout << "error, doesn't exist" << std::endl;
	}
	return res;
}

7.string转LPCWSTR

LPCWSTR stringToLPCWSTR(std::string orig)
{
	size_t origsize = orig.length() + 1;
	const size_t newsize = 100;
	size_t convertedChars = 0;
	wchar_t* wcstring = (wchar_t*)malloc(sizeof(wchar_t) * (orig.length() - 1));
	mbstowcs_s(&convertedChars, wcstring, origsize, orig.c_str(), _TRUNCATE);

	return wcstring;
}