gdal工具函数

参考

void setProjection(const std::string& sourceFile)
{
	GDALAllRegister();
	CPLSetConfigOption("SHAPE_ENCODING", "");			//解决中文乱码问题
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");	//解决中文路径崩溃

	char*				srs_str;
	OGRSpatialReference ogrSrs;
	ogrSrs.SetProjCS("UTM 50N /WGS84");
	ogrSrs.SetWellKnownGeogCS("WGS84");
	ogrSrs.SetUTM(50, true);
	ogrSrs.exportToPrettyWkt(&srs_str);

	// 创建输出数据集，格式为GeoTiff格式
	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("Gtiff");

	GDALDataset* pDataset = (GDALDataset*)(GDALOpen(sourceFile.c_str(), GA_Update));
	if (!pDataset)
	{
		std::cout << "读取tif数据失败！" << std::endl;
		return;
	}

	pDataset->SetProjection(srs_str);
	return;
}

//像素到投影坐标
osg::Vec3 pixel2Geo(GDALDataset* pDataset, int x, int y)
{
	double geoTransform[6];
	CPLErr err = pDataset->GetGeoTransform(geoTransform);

	osg::Matrixd m = osg::Matrixd::identity();
	m.set(geoTransform[1], geoTransform[4], 0.0, 0.0,
		geoTransform[2], geoTransform[5], 0.0, 0.0,
		0.0, 0.0, 1.0, 0.0,
		geoTransform[0], geoTransform[3], 0.0, 1.0);

	return osg::Vec3d(x, y, 0.0) * m;
}


osg::Vec3d transformation(const osg::Vec3d& pt, OGRSpatialReference* src, OGRSpatialReference* tagert)
{
	osg::Vec3d ret = pt;
	OGRCoordinateTransformation* trans = OGRCreateCoordinateTransformation(src, tagert);
	trans->Transform(1, &ret[0], &ret[1], &ret[2]);

	return ret;
}

void processDomGray(const std::string& sourceFile)
{
	GDALAllRegister();
	CPLSetConfigOption("SHAPE_ENCODING", "");			//解决中文乱码问题
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");	//解决中文路径崩溃

	// 创建输出数据集，格式为GeoTiff格式
	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("Gtiff");

	GDALDataset* pDataset = (GDALDataset*)(GDALOpen(sourceFile.c_str(), GA_Update));
	if (!pDataset)
	{
		std::cout << "读取影像数据失败！" << std::endl;
		return;
	}

	int numBands = pDataset->GetRasterCount();

	for (int i = 1; i <= numBands; i++)
	{
		GDALRasterBand* band = pDataset->GetRasterBand(i);
		band->SetColorInterpretation( GDALColorInterp(i + 2));
	}
}

void convertTifFromWorld2Projection(const std::string& sourceFile, const std::string& desFile)
{
	GDALAllRegister();
	CPLSetConfigOption("SHAPE_ENCODING", "");			//解决中文乱码问题
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");	//解决中文路径崩溃

	char*				srs_str;
	OGRSpatialReference ogrSrs;
	ogrSrs.SetProjCS("UTM 50N /WGS84");
	ogrSrs.SetWellKnownGeogCS("WGS84");
	ogrSrs.SetUTM(50, true);
	ogrSrs.exportToPrettyWkt(&srs_str);

	// 创建输出数据集，格式为GeoTiff格式
	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("Gtiff");

	GDALDataset* pDataset = ( GDALDataset* )(GDALOpen(sourceFile.c_str(), GA_Update));
	if (!pDataset)
	{
		std::cout << "读取tif数据失败！" << std::endl;
		return;
	}

	const int pixwidth  = pDataset->GetRasterXSize();
	const int pixheight = pDataset->GetRasterYSize();

	GDALDataset* pNewDataset = pDriver->Create(desFile.c_str(), pixwidth, pixheight, 1, GDT_Float32, NULL);
	pNewDataset->SetProjection(srs_str);

	double geoTransform[6];
	CPLErr err = pDataset->GetGeoTransform(geoTransform);

	// 设置六参数
	float* pData = new float[pixwidth * pixheight];
	memset(pData, 0, sizeof(float) * pixwidth * pixheight);

	pNewDataset->SetGeoTransform(geoTransform);

	int index = 0;
	auto band = selectband(pDataset);

	std::vector<osg::Vec3> PTs;
	for (int j = 0; j < pixheight; j++)
	{
		for (int i = 0; i < pixwidth; i++)
		{
			float res = -999;

			band->RasterIO(GF_Read, i, j, 1, 1, ( void* )(&res), 1, 1, GDT_Float32, 0, 0);
			pData[index++] = res;
		}
	}

	// 写入tif
	pNewDataset->RasterIO(GF_Write, 0, 0, pixwidth, pixheight, pData, pixwidth, pixheight, GDT_Float32, 1, 0, 0, 0, 0);
}