学习gdal的过程

代码片段1

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


using namespace std;


    int main()
    {
        printf("strat!\n");
        //注册文件格式
        GDALAllRegister();
        const char* file = "D:\\img1_1_32.tif";
        GDALDataset *ds = (GDALDataset*) GDALOpen (file, GA_ReadOnly);
        //A* p则使用：p->play(); 左边是结构指针。
        //A p 则使用：p.paly(); 左边是结构变量。
        int Xsize = ds->GetRasterXSize();
        int Ysize = ds->GetRasterYSize();
        printf("Xsize =%d\n" , Xsize);
        printf("影像长为%d,宽为%d,波段数为%d\n", ds->GetRasterXSize(), ds->GetRasterYSize(), ds->GetRasterCount());

        //读取影像第一波段
        GDALRasterBand* band = ds->GetRasterBand(1);
        int min,max;
        double minmax[2];
        //char a;
        //a = &min;
        minmax[0] = band->GetMinimum(&min);
        minmax[1] = band->GetMaximum(&max);
        printf("\n最小值是：%fd\n最大值是%f\n", minmax[0],minmax[1]);

    }

代码片段2

#include <iostream>
#include "gdal_priv.h"  // GDAL核心库头文件
#include "cpl_conv.h"   // 对于字符串处理函数

int main(int argc, char* argv[]) {
    if (argc != 2) {
        std::cerr << "Usage: " << argv[0] << " <raster-file>" << std::endl;
        return 1;
    }

    // 注册驱动
    GDALAllRegister();

    // 打开栅格文件
    GDALDataset *poDataset = static_cast<GDALDataset*>(GDALOpen(argv[1], GA_ReadOnly));
    if (poDataset == nullptr) {
        std::cerr << "Error: Cannot open raster file " << argv[1] << std::endl;
        return 1;
    }

    // 获取空间参考
    const  OGRSpatialReference *spatialRef = poDataset->GetSpatialRef();
    if (spatialRef != nullptr) {
        char *pszWKT = nullptr;
        spatialRef->exportToWkt(&pszWKT);  # 这里不理解
        std::cout << "WKT of the raster file:" << std::endl;
        std::cout << pszWKT << std::endl;

        // 释放分配的内存
        CPLFree(pszWKT);
    } else {
        std::cout << "No spatial reference system found for this raster." << std::endl;
    }

    // 关闭数据集
    GDALClose(poDataset);

    return 0;
}

分块计算NDVI

#include <iostream>
#include "gdal_priv.h"
#include <vector>
#include <cmath>
#include <stdexcept>
#include <chrono> // 引入计时库

const int BLOCK_SIZE = 256; // 分块大小

void calculateNDVI(float* nirData, float* redData, float* ndviData, int blockWidth, int blockHeight) {
    for (int y = 0; y < blockHeight; y++) {
        for (int x = 0; x < blockWidth; x++) {
            int index = y * blockWidth + x;
            float nir = nirData[index];
            float red = redData[index];

            // 计算NDVI
            if (nir + red != 0) {
                ndviData[index] = (nir - red) / (nir + red);
            } else {
                ndviData[index] = -1.0f; // 无效值
            }
        }
    }
}

int main() {
    // 初始化计时
    auto startTotal = std::chrono::high_resolution_clock::now(); // 总开始时间

    // 初始化GDAL
    GDALAllRegister();

    // 打开输入影像
    auto startOpen = std::chrono::high_resolution_clock::now(); // 打开文件开始时间
    GDALDataset* poDataset = (GDALDataset*)GDALOpen("D:\\temp\\GF1C_miluo2.tif", GA_ReadOnly);
    if (!poDataset) {
        std::cerr << "无法打开影像文件！" << std::endl;
        return 1;
    }
    auto endOpen = std::chrono::high_resolution_clock::now(); // 打开文件结束时间
    std::chrono::duration<double> elapsedOpen = endOpen - startOpen;
    std::cout << "打开文件耗时: " << elapsedOpen.count() << " 秒" << std::endl;

    // 获取影像信息
    int width = poDataset->GetRasterXSize();
    int height = poDataset->GetRasterYSize();
    int bandCount = poDataset->GetRasterCount();

    if (bandCount < 4) {
        std::cerr << "影像波段数不足，至少需要4个波段！" << std::endl;
        GDALClose(poDataset);
        return 1;
    }

    std::cout << "影像尺寸: " << width << "x" << height << std::endl;
    std::cout << "波段数: " << bandCount << std::endl;

    // 创建输出影像（单波段，存储NDVI）
    auto startCreate = std::chrono::high_resolution_clock::now(); // 创建输出文件开始时间
    GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* poDstDataset = poDriver->Create("ndvi_output.tif", width, height, 1, GDT_Float32, nullptr);

    // 设置NDVI影像的地理参考和投影信息
    double geoTransform[6];
    if (poDataset->GetGeoTransform(geoTransform) == CE_None) {
        poDstDataset->SetGeoTransform(geoTransform);
    }
    const char* projection = poDataset->GetProjectionRef();
    if (projection && strlen(projection) > 0) {
        poDstDataset->SetProjection(projection);
    }
    auto endCreate = std::chrono::high_resolution_clock::now(); // 创建输出文件结束时间
    std::chrono::duration<double> elapsedCreate = endCreate - startCreate;
    std::cout << "创建输出文件耗时: " << elapsedCreate.count() << " 秒" << std::endl;

    // 分块处理
    auto startProcess = std::chrono::high_resolution_clock::now(); // 分块处理开始时间
    for (int y = 0; y < height; y += BLOCK_SIZE) {
        for (int x = 0; x < width; x += BLOCK_SIZE) {
            // 计算当前块的实际大小
            int blockWidth = std::min(BLOCK_SIZE, width - x);
            int blockHeight = std::min(BLOCK_SIZE, height - y);

            // 分配内存
            std::vector<float> nirData(blockWidth * blockHeight);
            std::vector<float> redData(blockWidth * blockHeight);
            std::vector<float> ndviData(blockWidth * blockHeight);

            // 读取NIR波段（第4波段）
            GDALRasterBand* nirBand = poDataset->GetRasterBand(4);
            nirBand->RasterIO(GF_Read, x, y, blockWidth, blockHeight,
                              nirData.data(), blockWidth, blockHeight, GDT_Float32, 0, 0);

            // 读取Red波段（第3波段）
            GDALRasterBand* redBand = poDataset->GetRasterBand(3);
            redBand->RasterIO(GF_Read, x, y, blockWidth, blockHeight,
                              redData.data(), blockWidth, blockHeight, GDT_Float32, 0, 0);

            // 计算NDVI
            calculateNDVI(nirData.data(), redData.data(), ndviData.data(), blockWidth, blockHeight);

            // 写入NDVI数据
            GDALRasterBand* ndviBand = poDstDataset->GetRasterBand(1);
            ndviBand->RasterIO(GF_Write, x, y, blockWidth, blockHeight,
                               ndviData.data(), blockWidth, blockHeight, GDT_Float32, 0, 0);
        }
        std::cout << "已处理行: " << y << std::endl;
    }
    auto endProcess = std::chrono::high_resolution_clock::now(); // 分块处理结束时间
    std::chrono::duration<double> elapsedProcess = endProcess - startProcess;
    std::cout << "分块处理耗时: " << elapsedProcess.count() << " 秒" << std::endl;

    // 关闭数据集
    auto startClose = std::chrono::high_resolution_clock::now(); // 关闭文件开始时间
    GDALClose(poDataset);
    GDALClose(poDstDataset);
    auto endClose = std::chrono::high_resolution_clock::now(); // 关闭文件结束时间
    std::chrono::duration<double> elapsedClose = endClose - startClose;
    std::cout << "关闭文件耗时: " << elapsedClose.count() << " 秒" << std::endl;

    // 总耗时
    auto endTotal = std::chrono::high_resolution_clock::now(); // 总结束时间
    std::chrono::duration<double> elapsedTotal = endTotal - startTotal;
    std::cout << "总耗时: " << elapsedTotal.count() << " 秒" << std::endl;

    std::cout << "NDVI计算完成！结果已保存到 ndvi_output.tif" << std::endl;
    return 0;
}

对应的python代码

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2025/2/11 14:38 
# @File : block_compute_ndvi.py 
'''

'''
import osgeo.gdal as gdal
import numpy as np
from time import perf_counter

BLOCK_SIZE = 512  # 分块大小

def calculate_ndvi(nir_data, red_data):
    """计算NDVI"""
    ndvi = np.zeros_like(nir_data, dtype=np.float32)
    ndvi = (nir_data - red_data) / (nir_data + red_data)
    ndvi[np.isnan(ndvi)] = -1.0  # 处理无效值
    return ndvi

def main():
    start_total = perf_counter()  # 总开始时间

    # 打开输入影像
    start_open = perf_counter()  # 打开文件开始时间
    dataset = gdal.Open("D:\\temp\\GF1C_miluo2.tif", gdal.GA_ReadOnly)
    if not dataset:
        print("无法打开影像文件！")
        return 1
    end_open = perf_counter()  # 打开文件结束时间
    print(f"打开文件耗时: {end_open - start_open:.2f} 秒")

    # 获取影像信息
    width = dataset.RasterXSize
    height = dataset.RasterYSize
    band_count = dataset.RasterCount

    if band_count < 4:
        print("影像波段数不足，至少需要4个波段！")
        dataset = None
        return 1

    print(f"影像尺寸: {width}x{height}")
    print(f"波段数: {band_count}")

    # 创建输出影像（单波段，存储NDVI）
    start_create = perf_counter()  # 创建输出文件开始时间
    driver = gdal.GetDriverByName("GTiff")
    out_dataset = driver.Create("ndvi_output.tif", width, height, 1, gdal.GDT_Float32)

    # 设置NDVI影像的地理参考和投影信息
    geo_transform = dataset.GetGeoTransform()
    if geo_transform:
        out_dataset.SetGeoTransform(geo_transform)
    projection = dataset.GetProjection()
    if projection:
        out_dataset.SetProjection(projection)
    end_create = perf_counter()  # 创建输出文件结束时间
    print(f"创建输出文件耗时: {end_create - start_create:.2f} 秒")

    # 分块处理
    start_process = perf_counter()  # 分块处理开始时间
    for y in range(0, height, BLOCK_SIZE):
        for x in range(0, width, BLOCK_SIZE):
            block_width = min(BLOCK_SIZE, width - x)
            block_height = min(BLOCK_SIZE, height - y)

            # 读取NIR波段（第4波段）
            nir_band = dataset.GetRasterBand(4)
            nir_data = nir_band.ReadAsArray(x, y, block_width, block_height).astype(np.float32)

            # 读取Red波段（第3波段）
            red_band = dataset.GetRasterBand(3)
            red_data = red_band.ReadAsArray(x, y, block_width, block_height).astype(np.float32)

            # 计算NDVI
            ndvi_data = calculate_ndvi(nir_data, red_data)

            # 写入NDVI数据
            out_band = out_dataset.GetRasterBand(1)
            out_band.WriteArray(ndvi_data, x, y)

        print(f"已处理行: {y}")

    end_process = perf_counter()  # 分块处理结束时间
    print(f"分块处理耗时: {end_process - start_process:.2f} 秒")

    # 关闭数据集
    start_close = perf_counter()  # 关闭文件开始时间
    dataset = None  # 关闭输入影像
    out_dataset = None  # 关闭输出影像
    end_close = perf_counter()  # 关闭文件结束时间
    print(f"关闭文件耗时: {end_close - start_close:.2f} 秒")

    end_total = perf_counter()  # 总结束时间
    print(f"总耗时: {end_total - start_total:.2f} 秒")

    print("NDVI计算完成！结果已保存到 ndvi_output.tif")

if __name__ == "__main__":
    main()

获取矢量文件的空间参考信息

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

int main()
{
    // 注册所有 GDAL 驱动
    GDALAllRegister();

    // 指定要打开的 Shapefile 文件
    const char* pszFilename = "F:\\vector\\yili_vector.shp";
    GDALDataset* poDS = (GDALDataset*) GDALOpenEx(pszFilename, GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (poDS == NULL)
    {
        std::cerr << "打开文件失败：" << pszFilename << std::endl;
        return 1;
    }

    // 获取第一个图层
    OGRLayer* poLayer = poDS->GetLayer(0);
    if (poLayer == NULL)
    {
        std::cerr << "无法获取图层！" << std::endl;
        GDALClose(poDS);
        return 1;
    }

    // 获取图层的空间参考对象
    OGRSpatialReference* poSRS = poLayer->GetSpatialRef();
    if (poSRS == NULL)
    {
        std::cerr << "图层中没有空间参考信息。" << std::endl;
    }
    else
    {
        char* pszSRS_WKT = NULL;
        poSRS->exportToWkt(&pszSRS_WKT);
        std::cout << "空间参考信息：" << std::endl;
        std::cout << pszSRS_WKT << std::endl;
        CPLFree(pszSRS_WKT);
    }

    // 关闭数据集
    GDALClose(poDS);
    return 0;
}

将输入矢量数据转换到 WGS84 坐标系统，并输出到新文件

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

// 将输入矢量数据转换到 WGS84 坐标系统，并输出到新文件
bool transformVectorToWGS84(const std::string& inputFile,
                            const std::string& outputFile) {
    // 注册 GDAL 和 OGR 驱动
    GDALAllRegister();
    OGRRegisterAll();

    // 打开输入数据集（矢量数据）
    GDALDataset* poSrcDS = (GDALDataset*)GDALOpenEx(
            inputFile.c_str(),
            GDAL_OF_VECTOR,
            nullptr,
            nullptr,
            nullptr);
    if (poSrcDS == nullptr) {
        std::cerr << "无法打开输入文件: " << inputFile << std::endl;
        return false;
    }

    // 获取 ESRI Shapefile 驱动
    GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");
    if (poDriver == nullptr) {
        std::cerr << "无法加载 ESRI Shapefile 驱动" << std::endl;
        GDALClose(poSrcDS);
        return false;
    }

    // 创建输出数据集（新 Shapefile）
    GDALDataset* poDstDS = poDriver->Create(
            outputFile.c_str(),
            0, 0, 0,
            GDT_Unknown,
            nullptr);
    if (poDstDS == nullptr) {
        std::cerr << "无法创建输出文件: " << outputFile << std::endl;
        GDALClose(poSrcDS);
        return false;
    }

    // 定义目标空间参考（WGS84）
    OGRSpatialReference oTargetSRS;
    oTargetSRS.SetWellKnownGeogCS("WGS84");
    oTargetSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); // 强制经度-纬度顺序

    // 遍历输入数据集中的每个图层
    for (int iLayer = 0; iLayer < poSrcDS->GetLayerCount(); iLayer++) {
        OGRLayer* poSrcLayer = poSrcDS->GetLayer(iLayer);
        if (poSrcLayer == nullptr) {
            std::cerr << "11" << std::endl;
            continue;
        }

        // 获取输入图层的空间参考
        const OGRSpatialReference* poSrcSRS = poSrcLayer->GetSpatialRef();
        if (poSrcSRS == nullptr) {
            std::cerr << "图层中没有空间参考信息。" << std::endl;
            continue;
        }

        // 创建坐标转换对象：从输入空间参考转换到 WGS84
        OGRCoordinateTransformation* poCT = OGRCreateCoordinateTransformation(
                const_cast<OGRSpatialReference*>(poSrcSRS), &oTargetSRS);
        if (poCT == nullptr) {
            std::cerr << "无法创建坐标转换" << std::endl;
            continue;
        }

        // 创建输出图层，使用输入图层的名称、几何类型，并指定目标空间参考
        OGRLayer* poDstLayer = poDstDS->CreateLayer(
                poSrcLayer->GetName(),
                &oTargetSRS,
                poSrcLayer->GetGeomType(),
                nullptr);
        if (poDstLayer == nullptr) {
            std::cerr << "111" << std::endl;
            OGRCoordinateTransformation::DestroyCT(poCT);
            continue;
        }

        // 复制字段定义
        OGRFeatureDefn* poFDefn = poSrcLayer->GetLayerDefn();
        for (int i = 0; i < poFDefn->GetFieldCount(); i++) {
            poDstLayer->CreateField(poFDefn->GetFieldDefn(i));
        }

        // 遍历输入图层中的每个要素，并进行坐标转换
        OGRFeature* poSrcFeature = nullptr;
        poSrcLayer->ResetReading();
        while ((poSrcFeature = poSrcLayer->GetNextFeature()) != nullptr) {
            OGRGeometry* poGeometry = poSrcFeature->GetGeometryRef();
            if (poGeometry != nullptr) {
                // 执行几何对象的坐标转换
                if (poGeometry->transform(poCT) != OGRERR_NONE) {
                    std::cerr << "11" << std::endl;
                }
            }

            // 创建新的要素，复制属性信息，并添加到输出图层
            OGRFeature* poDstFeature = new OGRFeature(poDstLayer->GetLayerDefn());
            poDstFeature->SetFrom(poSrcFeature);
            if (poDstLayer->CreateFeature(poDstFeature) != OGRERR_NONE) {
                std::cerr << "1" << std::endl;
            }

//            OGRFeature::DestroyFeature(poDstFeature);
//            OGRFeature::DestroyFeature(poSrcFeature);
        }

        // 销毁坐标转换对象
        OGRCoordinateTransformation::DestroyCT(poCT);
    }

    // 关闭数据集，释放资源
    GDALClose(poSrcDS);
    GDALClose(poDstDS);

    return true;
}

// 示例：转换指定输入 Shapefile 到 WGS84 并保存为新的 Shapefile
int main() {
    std::string inputPath = "F:\\vector\\yili_vector.shp";
    std::string outputPath = "F:\\vector\\yili_vector_wgs84_cp.shp";

    if (transformVectorToWGS84(inputPath, outputPath)) {
        std::cout << "11" << std::endl;
    } else {
        std::cerr << "1111" << std::endl;
    }

    return 0;
}