c++ - 如何使用 GDAL 在 C++ 中计算给定矩形内具有特定值的光栅像素？

Question

在 C++ 中，是否有一个函数可以使用 GDAL 库计算给定矩形(xmin、ymin、xmax、ymax)内具有特定值的所有像素？或者我必须读取每个 block 并逐个像素地计算它们？我进行了搜索，但只找到了一些执行此操作的 Python 脚本。

我自己做了(它是一个 bool 栅格 - 0/1 - 我正在计算“1”像素)，但结果与 GRASS GIS r.report 函数(通过 QGIS)中给出的结果不同。

long long openTIF(string ftif,double x0,double y0,double x1,double y1) {
    long long sum = 0;
    GDALDataset *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset*)GDALOpen(ftif.c_str(),GA_ReadOnly);
    if (poDataset == NULL) {
        cout << "Error reading raster '" << ftif << "'\n";
        exit(1);
    }
    int tx=poDataset->GetRasterXSize(),
        ty=poDataset->GetRasterYSize();
    double geoTransf[6], rx0=0, ry0=0, rx1=0, ry1=0;
    if (poDataset->GetGeoTransform(geoTransf) == CE_None) {
        rx0 = geoTransf[0];
        ry0 = geoTransf[3]-ty*geoTransf[1];
        rx1 = geoTransf[0]-tx*geoTransf[5];
        ry1 = geoTransf[3];
    } else {
        exit(2);
    }
    int nBlockXSize,nBlockYSize;
    GDALRasterBand *poBand;
    poBand = poDataset->GetRasterBand(1);
    poBand->GetBlockSize(&nBlockXSize,&nBlockYSize);
    uint32_t i;
    int y,
        col0 = round((rx0-x0)/geoTransf[5]),
        col1 = round((rx0-x1)/geoTransf[5]),
        row0 = round((ry1-y1)/geoTransf[1]),
        row1 = round((ry1-y0)/geoTransf[1]);
    CPLErr error;
    uint8_t *data = (uint8_t*)CPLMalloc(nBlockXSize*nBlockYSize);
    GDALDataType type = poBand->GetRasterDataType();
    if (type == GDT_Byte) {
        cout << "Byte" << endl;
    }
    for (y=row0; y<row1; y++) {
        error = poBand->ReadBlock(0,y,data);
        if (error > 0) {
            cout << "Error reading data block.\n";
            exit(3);
        }
        for (i=col0; i<col1; i++) {
            sum += (uint8_t)data[i];
        }
    }
    CPLFree(data);
    return sum;
}

对于这个给定的栅格和给定的坐标(精度 %.8f)，GRASS GIS r.report 函数报告 28096011 个值为 1 的像素，而我的总和为 28094709(差异 = -1302)。使用其他坐标 r.report 给出 5605458，我的总和给出 5604757(差异 = -701)。知道会发生什么吗？

编辑:因为我的总和总是小于 GRASS GIS r.report，所以我考虑包括最后一行和最后一列，更改行

    for (y=row0; y<row1; y++) {
        for (i=col0; i<col1; i++) {

到

    for (y=row0; y<=row1; y++) {
        for (i=col0; i<=col1; i++) {

但是现在，使用另一组坐标，r.report 给出 249707，而我的总和给出 250157(差 = +450)。

Answer 1

这看起来像是一个舍入错误，您需要向下舍入您的列和行索引。如果您将每个光栅像素都视为一个矩形区域，则更容易推理。例如，栅格的第一个像素从 (x, y) 到 (x + width, y + height)。

rx = geoTransf[0];
rx_size = geoTransf[1]; // you got geoTransf[1] and geoTransf[5] swapped!
ry = geoTransf[3];    
ry_size = geoTransf[5];
// assert, that geoTrans[2] and geoTrans[4] are zero

col_first = floor((x0 - rx) / rx_size));
col_last = floor((x1 - rx) / rx_size));
row_first = floor((y0 - ry) / ry_size));
row_last = floor((y1 - ry) / ry_size));

for (int y = row_first; y <= row_last; ++y) {
    poBand->ReadBlock(col_first, y, data);
    for (int x = 0; x <= col_last - col_first; ++x) {
         data[x];
    }
}