Как подсчитать растровые пиксели с определенным значением внутри заданного прямоугольника в C ++ с помощью GDAL? - PullRequest
Новая
       42

Как подсчитать растровые пиксели с определенным значением внутри заданного прямоугольника в C ++ с помощью GDAL?

0 голосов
/ 05 декабря 2018

В C ++ есть функция для подсчета всех пикселей с определенным значением внутри данного прямоугольника (xmin, ymin, xmax, ymax), используя библиотеку GDAL?Или мне придется читать каждый блок и подсчитывать их все попиксельно?Я искал, но нашел только несколько скриптов Python, которые делают это.

Я сделал это для себя (это логический растр - 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++) {

to 

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

но теперь, с другим набором координат, r.report дает 249707, в то время как моя сумма дает 250157 (разница = +450).

1 Ответ

0 голосов
/ 18 декабря 2018

Похоже на ошибку округления, вам нужно округлить индексы столбцов и строк.Проще рассуждать, если думать о каждом растровом пикселе как о прямоугольной области.Например, первый пиксель вашего растра охватывает от (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];
    }
}
...