使用ChatGPT写一个GDALDatasetRasterIOEx功能增强函数

GDALDatasetRasterIOEx函数是GDAL用来读取栅格数据文件像素数据的函数,默认是不支持参数传入的读取窗口范围超出栅格的实际范围的,但是很多时候,需要按照地理范围读取,就会有超出栅格范围的情况。
之前我常用的做法,就是创建一个 基于 MemVRT(虚拟栅格),VRT的地理范围取自己想要读取的范围即可,这种做法实际上是使用了GDAL自身的功能进行的一些简单的组合应用。
ChatGPT出来之后,简单尝试了一下,用它来写代码确实可以减少很大的工作量,也解决不少变量起名的问题。

下面是使用 ChatGPT 来编写一个函数,解决这个问题的记录。

写一个C++代码程序,封装一下GDALDatasetRasterIOEx 函数,使其支持访问范围超出影像大小,超出部分使用nodata值填充

以下是一个示例代码,封装了GDALDatasetRasterIOEx函数,支持访问范围超出影像大小,超出部分使用nodata值填充。在代码中,我们首先获取影像的nodata值和像素宽高等信息,然后再根据读取范围对应的实际影像坐标进行判断,判断是否超出影像范围。如果超出影像范围,则填充nodata值。

下面是 ChatGPT 写出来的代码,这个代码整体结构很清晰,变量命名也非常切意,不过还需要进行一下优化。

#include "gdal_priv.h"
#include "cpl_conv.h" // for CPLMalloc()

// 封装GDALDatasetRasterIOEx函数
CPLErr customGDALDatasetRasterIOEx(GDALDataset *pDataset, GDALRWFlag eRWFlag,
                                   int nXOff, int nYOff, int nXSize, int nYSize,
                                   void *pData, int nBufXSize, int nBufYSize,
                                   GDALDataType eBufType, int nBandCount,
                                   int *panBandMap, int nPixelSpace, int nLineSpace, int nBandSpace,
                                   const char* pszOrder, int nCompeletedPixels, int nCompeletedLines)
{
    if (pDataset == NULL) {
        return CE_Failure;
    }

    CPLErr eErr = CE_None;

    // 获取像素和行之间的字节偏移量
    int nPixelBytes = GDALGetDataTypeSize(eBufType) / 8;
    int nLineBytes = nPixelBytes * nBufXSize;

    // 获取影像的像素、行数和波段数等信息
    int nRasterXSize = pDataset->GetRasterXSize();
    int nRasterYSize = pDataset->GetRasterYSize();
    int nBands = pDataset->GetRasterCount();

    // 获取nodata值
    double dfNoDataValue;
    int bGotNoData = false;
    if (nBandCount == 1) {
        GDALRasterBand *pBand = pDataset->GetRasterBand(panBandMap[0]);
        dfNoDataValue = pBand->GetNoDataValue(&bGotNoData);
    }

    // 计算读取范围对应的实际影像坐标
    int nActualXOff = nXOff < 0 ? 0 : nXOff;
    int nActualYOff = nYOff < 0 ? 0 : nYOff;
    int nActualXEnd = nXOff + nXSize > nRasterXSize ? nRasterXSize : nXOff + nXSize;
    int nActualYEnd = nYOff + nYSize > nRasterYSize ? nRasterYSize : nYOff + nYSize;

    // 计算超出部分的长度和宽度
    int nExceedXStart = nXOff < 0 ? -nXOff : 0;
    int nExceedYStart = nYOff < 0 ? -nYOff : 0;
    int nExceedXEnd = nXOff + nXSize > nRasterXSize ? nXOff + nXSize - nRasterXSize : 0;
    int nExceedYEnd = nYOff + nYSize > nRasterYSize ? nYOff + nYSize - nRasterYSize : 0;

    // 实际要读取的像素宽度和高度
    int nActualXSize = nActualXEnd - nActualXOff;
    int nActualYSize = nActualYEnd - nActualYOff;

    // 创建一个临时数组存储影像数据
    char *pTempData = (char*)CPLMalloc(nPixelBytes * nBufXSize * nBufYSize * nBands);

    // 读取实际部分
    if (nActualXSize > 0 && nActualYSize > 0) {
        eErr = pDataset->RasterIO(eRWFlag, nActualXOff, nActualYOff, nActualXSize, nActualYSize,
                                   pTempData + (nExceedXStart + nExceedYStart * nBufXSize) * nPixelBytes * nBands,
                                   nBufXSize, nBufYSize, eBufType, nBandCount, panBandMap,
                                   nPixelSpace, nLineSpace, nBandSpace);
    }

    // 填充超出部分的nodata值
    for (int i = 0; i < nExceedYStart; i++) {
        char *pLineData = pTempData + i * nLineBytes;
        for (int j = 0; j < nBufXSize; j++) {
            for (int k = 0; k < nBands; k++) {
                if (bGotNoData) {
                    GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
                                  pLineData + (j + k * nBufXSize) * nPixelBytes, eBufType, 0,
                                  1);
                } else {
                    memset(pLineData + (j + k * nBufXSize) * nPixelBytes, 0,
                           nPixelBytes);
                }
            }
        }
    }

    for (int i = nExceedYStart; i < nBufYSize - nExceedYEnd; i++) {
        char *pLineData = pTempData + i * nLineBytes;
        for (int j = 0; j < nExceedXStart; j++) {
            for (int k = 0; k < nBands; k++) {
                if (bGotNoData) {
                    GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
                                  pLineData + (j + k * nBufXSize) * nPixelBytes, eBufType, 0,
                                  1);
                } else {
                    memset(pLineData + (j + k * nBufXSize) * nPixelBytes, 0,
                           nPixelBytes);
                }
            }
        }

        for (int j = nActualXOff; j < nActualXEnd; j++) {
            char *pPixelData = pLineData + (j - nActualXOff + nExceedXStart) * nPixelBytes * nBands;
            eErr = pDataset->RasterIO(eRWFlag, j, i, 1, 1,
                                       pPixelData, 1, 1, eBufType, nBandCount, panBandMap,
                                       nPixelSpace, nLineSpace, nBandSpace);
        }

        for (int j = nBufXSize - nExceedXEnd; j < nBufXSize; j++) {
            for (int k = 0; k < nBands; k++) {
                if (bGotNoData) {
                    GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
                                  pLineData + (j + k * nBufXSize) * nPixelBytes, eBufType, 0,
                                  1);
                } else {
                    memset(pLineData + (j + k * nBufXSize) * nPixelBytes, 0,
                           nPixelBytes);
                }
            }
        }
    }

    for (int i = nBufYSize - nExceedYEnd; i < nBufYSize; i++) {
        char *pLineData = pTempData + i * nLineBytes;
        for (int j = 0; j < nBufXSize; j++) {
            for (int k = 0; k < nBands; k++) {
                if (bGotNoData) {
                    GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
                                  pLineData + (j + k * nBufXSize) * nPixelBytes, eBufType, 0,
                                  1);
                } else {
                    memset(pLineData + (j + k * nBufXSize) * nPixelBytes, 0,
                           nPixelBytes);
                }
            }
        }
    }

    // 将临时数组中的数据写入到输出指针所指向的空间中
    memcpy(pData, pTempData, nPixelBytes * nBufXSize * nBufYSize * nBands);

    CPLFree(pTempData);

    return eErr;
}
posted @ 2023-04-10 14:49  乌合之众  阅读(125)  评论(1编辑  收藏  举报
clear