Как подсчитать растровые пиксели с определенным значением внутри заданного прямоугольника в C ++ с помощью GDAL?
В 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++) {
в
for (y=row0; y<=row1; y++) {
for (i=col0; i<=col1; i++) {
но теперь, с другим набором координат, r.report дает 249707, а моя сумма дает 250157 (разница = +450).
Решение
Это похоже на ошибку округления, вам нужно округлить индексы столбцов и строк. Проще рассуждать, если думать о каждом растровом пикселе как о прямоугольной области. Например, первый пиксель вашего растра простирается от (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];
}
}
Другие решения
Других решений пока нет …