1. 程式人生 > 實用技巧 >GDAL資料集寫入空間座標參考

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類即可。