学习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()