GDAL数据集写入空间坐标参考

1. 概述

可以通过GDAL给地理数据写入空间参考信息，不过要注意的是GDAL给矢量数据和栅格数据写入空间坐标参考的接口不太一样。

2. 栅格数据

实现代码如下：

#include <iostream>
#include <gdal_priv.h>
#include <string>

using namespace std;

int main()
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路径
	CPLSetConfigOption("SHAPE_ENCODING", "");  //解决中文乱码问题	
	CPLSetConfigOption("GDAL_DATA", "D:/Work/GDALBuild/gdal-2.4.2/install/data");

	GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTIFF"); //图像驱动
	char** ppszOptions = NULL;
	ppszOptions = CSLSetNameValue(ppszOptions, "BIGTIFF", "IF_NEEDED"); //配置图像信息
	const char* dstPath = "dst.tif";
	GDALDataset* dst = pDriver->Create(dstPath, 256, 256, 3, GDT_Byte, ppszOptions);
	if (dst == nullptr)
	{
		printf("Can't Write Image!");
		return false;
	}

	//空间参考
	OGRSpatialReference spatialReference;
	spatialReference.importFromEPSG(4326);				//wgs84地理坐标系
	char *pszWKT = nullptr;
	spatialReference.exportToWkt(&pszWKT);

	dst->SetProjection(pszWKT);

	CPLFree(pszWKT);
	pszWKT = nullptr;

	//坐标信息
	double padfTransform[6] = {
		114.0,		//左上角点坐标X
		0.000001,		//X方向的分辨率
		0,		//旋转系数，如果为0，就是标准的正北向图像
		34.0,			//左上角点坐标Y
		0,			//旋转系数，如果为0，就是标准的正北向图像
		0.000001,			//Y方向的分辨率
	};
	dst->SetGeoTransform(padfTransform);

	GDALClose(dst);
}

这里创建了一个wgs84地理坐标系空间参考的栅格数据，通过OGRSpatialReference类导出了描述空间参考的wkt字符串，写入到GDAL数据集中。

3. 矢量数据

实现代码如下：

#include <iostream>
#include <gdal_priv.h>
#include <ogrsf_frmts.h>

using namespace std;

int main()
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路径
	CPLSetConfigOption("SHAPE_ENCODING", "");  //解决中文乱码问题	
	CPLSetConfigOption("GDAL_DATA", "D:/Work/GDALBuild/gdal-2.4.2/install/data");

	//空间参考
	OGRSpatialReference spatialReference;
	spatialReference.importFromEPSG(4326);				//wgs84地理坐标系
	
	//创建
	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");
	if (!driver)
	{
		printf("Get Driver ESRI Shapefile Error！\n");
		return false;
	}

	GDALDataset* dataset = driver->Create("dst.shp", 0, 0, 0, GDT_Unknown, NULL);
	OGRLayer* poLayer = dataset->CreateLayer("houseType", &spatialReference, wkbPolygon, NULL);

	//创建属性字段
	{
		// 字符串
		OGRFieldDefn oField1("名称", OFTString);
		oField1.SetWidth(8);
		if (poLayer->CreateField(&oField1) != OGRERR_NONE) {
			printf("Creating Name field failed.\n"); return FALSE;
		}

		// 浮点数
		OGRFieldDefn oField2("面积", OFTReal);
		oField2.SetPrecision(3);
		if (poLayer->CreateField(&oField2) != OGRERR_NONE) {
			printf("Creating Name field failed.\n"); return FALSE;
		}

		// 整型
		OGRFieldDefn oField3("结点数", OFTInteger);
		if (poLayer->CreateField(&oField3) != OGRERR_NONE) {
			printf("Creating Name field failed.\n"); return FALSE;
		}
	}

	//创建特征
	OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn());

	OGRLinearRing ogrring;
	int pNum = 4;
	ogrring.setNumPoints(pNum);
	ogrring.setPoint(0, 114.0, 34.0, 0.0);
	ogrring.setPoint(1, 115.0, 34.0, 0.0);
	ogrring.setPoint(2, 115.0, 35.0, 0.0);
	ogrring.setPoint(3, 114.0, 35.0, 0.0);
	   
	OGRPolygon polygon;
	polygon.addRing(&ogrring);
	poFeature->SetGeometry(&polygon);

	poFeature->SetField("名称", "多边形");
	poFeature->SetField("面积", polygon.get_Area());
	poFeature->SetField("结点数", pNum);

	if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
	{
		printf("Failed to create feature in shapefile.\n");
		return false;
	}

	//释放
	GDALClose(dataset);
	dataset = nullptr;
}

与写入到栅格数据不同，空间参考信息写入到矢量数据是写入到GDAL数据集的图层类中的，并且直接传入OGRSpatialReference类即可。