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;
}