电子海图 GDAL解析

C++调用GDAL实现栅格数据的读取

#include <iostream>
#include <gdal_priv.h>
#include <gdal_alg_priv.h>
#include <gdal.h>
using namespace std;
 
int main() {
	//注册所有的驱动
	GDALAllRegister();
	//设置支持中文路径和文件名
	//CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); 
	
	//打开文件
	string file_path_name = "H:/DEM/DEM12.5.tif";
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(file_path_name.c_str(), GA_ReadOnly);
	if (poDataset == NULL)
	{
		cout << "指定的文件不能打开!" <<endl;
		return 0;
	}
 
	// 获取波段信息
	GDALRasterBand* poBand = poDataset->GetRasterBand(1);
	int nXSize = poDataset->GetRasterXSize();
	int nYSize = poDataset->GetRasterYSize();
 
	cout << "width = " << nXSize << endl;
	cout << "height = " << nYSize << endl;
 
	// 获取投影参考系
	string projection = poDataset->GetProjectionRef();
	cout << "Projection = " << projection << endl;
 
	// 获取投影坐标范围
	double geoTransform[6];
	poDataset->GetGeoTransform(geoTransform);
	double minX = geoTransform[0];
	double minY = geoTransform[3] + nXSize * geoTransform[4] + nYSize * geoTransform[5];
	double maxX = geoTransform[0] + nXSize * geoTransform[1] + nYSize * geoTransform[2];
	double maxY = geoTransform[3];
	cout.setf(ios::fixed);
	cout.precision(8);
	cout << "X-axis minimum:" << minX << ",\nY-axis minimum:" << minY << ",\nX-axis maximum:" << maxX << ",\nY-axis maximum:" << maxY << endl;
 
    // 读取数据
	int* pafScanline = new int[nXSize];
	for (int i = 0; i < nYSize; i++) {
		poBand->RasterIO(GF_Read, 0, i, nXSize,1, pafScanline, nXSize, 1, GDT_Int16, 0, 0);
	}
 
	delete[] pafScanline;
 
	//关闭栅格文件
	GDALClose(poDataset);
 
	return 0;
}

---

GDAL api读写shp文件的方法

这里注意新版本的GDAL用GDALDataset代替OGRDataSource，头文件里也有描述，说OGRDataSource将被抛弃。

#include "GDAL/gdal.h"
#include "GDAL/gdal_priv.h"
#include "GDAL/ogrsf_frmts.h"
void testWriteShp()
{
	GDALAllRegister();

	const char *pFileName = "test.shp";
	const char *pszDriverName = "ESRI Shapefile";
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName);
	if (poDriver == NULL)
		return;

	GDALDataset *poDs = poDriver->Create(pFileName, 0, 0, 0, GDT_Unknown, NULL); // 需要close
	if (poDs == NULL)
		return;

	// wkbPolygon表示绘制多边形；wkbLinearRing表示绘制线环
	OGRLayer *poLayer = poDs->CreateLayer("ring", NULL, wkbMultiLineString, NULL); // 不删
	if (poLayer == NULL)
	{
		return;
	}

	OGRFieldDefn idField("id", OFTInteger64);
	if (poLayer->CreateField(&idField) != OGRERR_NONE)
	{
		return;
	}

	OGRFieldDefn nameField("Name", OFTString);
	nameField.SetWidth(32);
	if (poLayer->CreateField(&nameField) != OGRERR_NONE)
	{
		return;
	}

	// 设置字段并加入几何图形
	OGRFeature *poFeature;

	// 第1个图形和属性
	poFeature = OGRFeature::CreateFeature(poLayer->GetLayerDefn());
	poFeature->SetField("id", (GIntBig)1);
	poFeature->SetField("Name", "ring");
	// 绘制线环，内部不标记颜色.线的绘制原则是，从第一个点开始，逐次绘制直到最后一个点。
	// 对于闭合的线环，由于首点和尾点是一个点。所以4边形绘制需要5个点，首点被重复一次
	OGRLineString line;
	line.setNumPoints(5);
	line.setPoint(0, 0, 1200, 0); //第一个输入变量是序号，第二、三、四分别是线XYZ坐标
	line.setPoint(1, 0, 0, 0);
	line.setPoint(2, 1200, 0, 0);
	line.setPoint(3, 1200, 1200, 0);
	line.setPoint(4, 0, 1200, 0);
	poFeature->SetGeometry(&line);
	if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
	{
		return;
	}
	OGRFeature::DestroyFeature(poFeature);

	// 第2个图形和属性
	poFeature = OGRFeature::CreateFeature(poLayer->GetLayerDefn());
	poFeature->SetField("id", (GIntBig)2);
	poFeature->SetField("Name", "poly");
	// 绘制多边形。多边形内部亦被颜色标记
	OGRPolygon poly;
	std::vector<OGRPoint> vec;
	vec.push_back(OGRPoint(0, 200));
	vec.push_back(OGRPoint(200, 200));
	vec.push_back(OGRPoint(200, 0));
	vec.push_back(OGRPoint(0, 0));
	OGRLinearRing ring;
	for (int i = 0; i < 4; i++)
	{
		ring.addPoint(vec[i].getX(), vec[i].getY());
	}
	ring.closeRings();
	poly.addRing(&ring);
	poFeature->SetGeometry(&poly);
	if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
	{
		return;
	}
	OGRFeature::DestroyFeature(poFeature);

	// 第3个图形和属性
	poFeature = OGRFeature::CreateFeature(poLayer->GetLayerDefn());
	poFeature->SetField("id", (GIntBig)3);
	poFeature->SetField("Name", "multiLine");
	// 多条线
	OGRMultiLineString multiLine;
	OGRLineString line1, line2;
	line1.setNumPoints(5);
	line1.setPoint(0, 0, 1000, 0); //第一个输入变量是序号，第二、三、四分别是线XYZ坐标
	line1.setPoint(1, 0, 0, 0);
	line1.setPoint(2, 1000, 0, 0);
	line1.setPoint(3, 1000, 1000, 0);
	line1.setPoint(4, 0, 1000, 0);
	line2.setNumPoints(5);
	line2.setPoint(0, 500, 1500, 0); //第一个输入变量是序号，第二、三、四分别是线XYZ坐标
	line2.setPoint(1, 500, 500, 0);
	line2.setPoint(2, 1500, 500, 0);
	line2.setPoint(3, 1500, 1500, 0);
	line2.setPoint(4, 500, 1500, 0);
	multiLine.addGeometry((const OGRGeometry *)&line1);
	multiLine.addGeometry((const OGRGeometry *)&line2);
	poFeature->SetGeometry(&multiLine);
	if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
	{
		return;
	}
	OGRFeature::DestroyFeature(poFeature);

	GDALClose(poDs);

	OGRCleanupAll();
}

void testReadShp()
{
	OGRRegisterAll();

	GDALDataset *poDS = (GDALDataset *)GDALOpenEx("test.shp", GDAL_OF_READONLY, NULL, NULL, NULL);
	if (poDS == NULL)
	{
		printf("Open failed.\n");
		return;
	}

	OGRLayer *poLayer = poDS->GetLayer(0);

	OGRFeature *poFeature;

	poLayer->ResetReading();
	while ((poFeature = poLayer->GetNextFeature()) != NULL)
	{
		OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
		int iField;

		for (iField = 0; iField < poFDefn->GetFieldCount(); iField++)
		{
			OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
			if (poFieldDefn->GetType() == OFTInteger)
				printf("%d, ", poFeature->GetFieldAsInteger(iField));
			else if (poFieldDefn->GetType() == OFTInteger64)
				printf("%lld, ", poFeature->GetFieldAsInteger64(iField));
			else if (poFieldDefn->GetType() == OFTReal)
				printf("%.3f, ", poFeature->GetFieldAsDouble(iField));
			else if (poFieldDefn->GetType() == OFTString)
				printf("%s, ", poFeature->GetFieldAsString(iField));
			else
				printf("%s, ", poFeature->GetFieldAsString(iField));
		}
		printf("\n");

		OGRGeometry *poGeometry = poFeature->GetGeometryRef();
		if (poGeometry != NULL)
		{
			switch (poGeometry->getGeometryType())
			{
			case wkbLineString:
			{
				OGRLineString *line = (OGRLineString *)poGeometry;
				for (int i = 0; i < line->getNumPoints(); i++)
				{
					printf("wkbLineString %d: x=%g y=%g z=%g\n", i, line->getX(i), line->getY(i), line->getZ(i));
				}
				break;
			}
			case wkbPolygon:
			{
				OGRPolygon *poly = (OGRPolygon *)poGeometry;
				OGRLinearRing *ring = (OGRLinearRing *)poly->getExteriorRingCurve();
				for (int i = 0; i < ring->getNumPoints(); i++)
				{
					printf("wkbPolygon %d: x=%g y=%g z=%g\n", i, ring->getX(i), ring->getY(i), ring->getZ(i));
				}
				break;
			}
			case wkbMultiLineString:
			{
				OGRMultiLineString *multiLine = (OGRMultiLineString *)poGeometry;
				for (int n = 0; n < multiLine->getNumGeometries(); n++)
				{
					OGRLineString *line = (OGRLineString *)multiLine->getGeometryRef(n);
					for (int i = 0; i < line->getNumPoints(); i++)
					{
						printf("wkbMultiLineString %d-%d: x=%g y=%g z=%g\n", n, i, line->getX(i), line->getY(i), line->getZ(i));
					}
				}

				break;
			}
			default:
				printf("no point geometry\n");
				break;
			}
		}

		OGRFeature::DestroyFeature(poFeature);
	}

	GDALClose(poDS);
}

int main(int argc, char const *argv[])
{
	testWriteShp();
	testReadShp();

	return 0;
}