C++. Перевод в пиксели

[pirateuse]

Здравствуйте.
Пришлось столкнуться с этим монстром - GDAL...
Вот есть у меня массив - изображение. Знаю, что изображение в проекции UTM (зона, гемисфера). Знаю геогр.координаты углов изображения, центра.
Есть у меня вектора координат точек (широта,долгота) в проекции WGS84. Надо поставить точки на изображение.
Пример того, что есть:

Код: Выделить всё

void Points2Img(uchar* img, int h, int w, GeoBox box, vector<double> &Lat, vector<double> &Lon)
{
	CPLSetConfigOption("GDAL_DATA", "..\\..\\gdal-data");
	OGRSpatialReference oSourceSRS, oTargetSRS;
	OGRCoordinateTransformation *poCT;

	oSourceSRS.importFromEPSG(4326); //WGS84 точек
	oTargetSRS.importFromEPSG(32628);//UTM изображения
	poCT = OGRCreateCoordinateTransformation(&oSourceSRS, &oTargetSRS);

	//перевод геогр.координат верхнего левого и нижнего правого углов в координаты UTM
	poCT->Transform(1, &box.left_top.x, &box.left_top.y);
	poCT->Transform(1, &box.rght_bot.x, &box.rght_bot.y);

	//сколько единиц в одном пикселе
	double unit_in_pixY = abs(box.left_top.y - box.rght_bot.y) / (h - 1);
	double unit_in_pixX = abs(box.rght_bot.x - box.left_top.x) / (w - 1);

	//перепроецируем координаты точек из WGS84 в заданный UTM
	poCT->Transform(Lat.size(), &Lon[0], &Lat[0]);

	//проставляем точки на изображении
	int y, x;
	for (size_t pp = 0; pp < Lat.size(); ++pp)
	{
		y = round(double(h-1) - (Lat[pp] - box.rght_bot.y) / unit_in_pixY); //ось Y на изображение направлена на "юг"
		x = round((Lon[pp] - box.left_top.x) / unit_in_pixX);
		if(проверка y,x)
			img[y*h + x] = 255;
	}
}

Перевод в пиксели производится в цикле из соображения, что при перепроецировании точек в UTM, координаты переходят в прямоугольную систему координат.
Вопрос в том, так ли это? У меня возникли сомнения.

Начал копать в другом направлении:

Код: Выделить всё

GDALDataset *pSrcDataset = (GDALDataset *)GDALOpen("tst.bmp", GA_ReadOnly);
OGRSpatialReference oSRS_proj;
oSRS_proj.importFromEPSG(32628);
char *pszSRS_WKT = 0;
oSRS_proj.exportToWkt(&pszSRS_WKT);
CPLErr err = pSrcDataset->SetProjection(pszSRS_WKT);

double adfGeoTransform[6] = {0, 0, 0, 0, 0, 0};
pSrcDataset->GetGeoTransform(adfGeoTransform);
double adfInvGeoTransform[6] = {0, 0, 0, 0, 0, 0};
GDALInvGeoTransform(adfGeoTransform, adfInvGeoTransform);
//......................
GDALApplyGeoTransform(adfInvGeoTransform, Lon[i],Lat[i], &new_x, &new_y);

Но возникает проблема: pSrcDataset->GetGeoTransform(adfGeoTransform) возвращает {0, 1, 0, 0, 0, 1}.
Вывод: pSrcDataset->SetProjection(pszSRS_WKT) не приводит к формированию коэффициентов.
В итоге не понятно как получить adfGeoTransform...

[Александр Мурый]

По сути, вам надо растеризовать точки и совместить два растра, так? "Пропатчить" изображение значениями точек?

[pirateuse]

По сути, да. Перенести точки на растр с учетом того, что растр и точки в разных проекциях (UTM и WGS84).

С другой стороны, неплохо было бы также уметь переводить из пиксельных координат растра в географические. Но одно следует из другого...

[Донецков]

А это именно нужно сделать в GDAL, нельзя ли их совместить в том же QGIS?

[pirateuse]

GDAL в любом случае используется, не хотелось бы подключать какие-то дополнительные библиотеки, которые опять же содержат те же самые модули для решения тех же задач... Да и разбираться надо будет дополнительно, как бы там тоже проблемы не возникли.
В общем, да, надо разобраться до конца с GDAL.

[Александр Мурый]

Почему бы не использовать возможности имеющегося API? Например, есть функции GDALRasterizeGeometries и GDALRasterizeLayers. Про перепроецирование растров написано здесь.
Я попытался бы сделать так:
- перепроецируем точки в СК растра;
- "прожигаем" растр значениями точек с помощью GDALRasterizeLayers. Причём пункт №1, возможно, можно будет пропустить, т.к. в этой функции есть параметр pfnTransformer для трасформации векторов в СК растра (как я понял).

[pirateuse]

Александр Мурый, спасибо за наводку. Посмотрю что за функции. Но мне так видятся две проблемы:
1) мой растр, загруженный с помощью

Код: Выделить всё

GDALDataset *pSrcDataset = (GDALDataset *)GDALOpen("tst.bmp", GA_ReadOnly);

не имеет географической привязки и прочей информации. Проекцию я задаю с помощью SetProjection, для привязки, наверное, SetGCPs... Надо пробовать...
2) хочется получить вектора с координатами точек отдельно, а не наносить их сразу на растр.

Все-таки правда ли, что координаты перепроецированных точек из градусов (из WGS84) в метры (в UTM) находятся в прямоугольной системе координат? Ведь тогда эту систему можно легко связать с прямоугольной системой координат изображения (координаты углов я знаю).
Либо мои точки после перепроекции не расположены в прямоугольной системе координат и тогда мне не получится обойтись без GetGeoTransform, GDALInvGeoTransform... В этом случае у меня такие данные: картинка, проекция картинки, геогр.координаты углов. При этом эти данные лежат отдельно друг от друга, надо ручками все засовывать в GDALDataset, о чем уже упомянул...

[Дмитрий Барышников]

Честно хочется помочь, но я запутался в ваших вопросах.
Вот выцепил первый:

Но возникает проблема: pSrcDataset->GetGeoTransform(adfGeoTransform) возвращает {0, 1, 0, 0, 0, 1}.
Вывод: pSrcDataset->SetProjection(pszSRS_WKT) не приводит к формированию коэффициентов.
В итоге не понятно как получить adfGeoTransform...

Ответ: ваш исходный растр не имеет географической привязки. В этом легко убедиться при помощи команды gdalinfo.
Как только привязка появится, GetGeoTransform будет возвращать иные значения.

[pirateuse]

Дмитрий Барышников, я и сам запутался в этом. Начал я с задачи перевода координат точек из одной проекции в другую, для этого подошел OGR, задача была решена достаточно быстро.
Но далее потребовалось нанести эти точки на растр, т.е. перевод в пиксельные координаты. Фактически географическая привязка растра мне известна. Наверно, если как-то задать нужные параметры в OGRSpatialReference, то задача будет решена, но как - не догадываюсь.
Поэтому пришлось обратиться к классу GDALDataset... Но в этом случае необходимо открывать растр средствами GDAL и вручную задавать параметры. Как раз по поводу географической привязки смотрю в сторону функции SetGCPs...

АП:
GDALOpen(...) - успех
SetProjection(...) - успех (GetProjectionRef() возвращает параметры проекции)
SetGCPs(...) - успех (GetGCPs() и GetGCPProjection() возвращают данные)
НО: GetProjectionRef() теперь возвращает пустую строку (так ли это должно быть?)
GetGeoTransform(...) опять возвращает {0, 1, 0, 0, 0, 1}
Да, все верно... Растр задается либо аффинными преобразованиями, либо через GCP. Одно другое исключает. Но как мне тогда задать географическую привязку и при этом получить параметры GeoTransform для перевода в пиксели?
АП: Хм, в GDAL нет механизма преобразования на основе опорных точек (GCPs). Значит путь геопривязки по координатам углов - тупиковый. Да и вообще похоже, что не получится получить параметры через GetGeoTransform(...), если загружаемый растр изначально не имеет привязки.

[Александр Мурый]

АП: Хм, в GDAL нет механизма преобразования на основе опорных точек (GCPs).

Откуда такая дезинформация? См. статью.

[pirateuse]

Откуда такая дезинформация? См. статью.

Нашел здесь: http://gdal.org/gdal_datamodel.html

The GDAL data model does not imply a transformation mechanism that must be generated from the GCPs ... this is left to the application. However 1st to 5th order polynomials are common.

Сейчас посмотрю вашу ссылку... Возможно, я неправильно выразился. Т.е. геопривязку задать можно через опорные точки, потом можно выполнять какие-то операции преобразования этого растра. Но как решать мою задачу перевода в пиксельные координаты - не понятно, так как GetGeoTransform не возвращает необходимые значения.

Да, в обсуждении приведенной статьи прозвучал вопрос

Очень классная статья! Спасибо большое.
Единственное, что мне осталось непонятным - как после привязки карты по контрольным точкам узнать географические координаты(долгота, широта) произвольной точки растрового изображения, зная ее пиксельные координаты.
Заранее спасибо.

Естественно, "просто загрузить в любую GIS, и посмотреть координаты курсора" (часть ответа от Максим Дубинин) мне не подходит. Надо это сделать "программно"...

[Дмитрий Барышников]

Давайте скорректируем вашу задачу, дабы избежать путаницы.
1) Для преобразования из географических координат в пиксельные OGRSpatialReference не подходит. Он вообще не может ничего преобразовывать. Скорее OGRCoordinateTransofrmation - да и то, только для преобразования между различными системами координат.
2) Для преобразования из пиксельных координат в географические (имеется в виду проекционные и географические) и обратно можно использовать GDALApplyGeoTransform или пересчитывать самостоятельно через матрицу.
3) Для использования GDALApplyGeoTransform нужно что бы исходное изображение было корректно привязано (gdalinfo должен выдавать информацию о привязке). В принципе привязку модно делать и программно.
4) Если известны координаты углов изображения совершенно не нужно использовать GCP. GCP нужно использовать если у вас есть координаты точек на местности и вы из можете идентифицировать на изображении. При этом следует понимать, что привязка может и не получиться - например, когда уравнение (полином) не имеет решения.
5) Для перехода от координат углов нужно только:
координаты верхнего левого угла + перейти от размеров изображения в географических координатах и количество пикселов по осям Х и Y из которых вы перейдете к размеру пиксела по X и Y.

Итого:
double dfULX, dfULY, dfURX, dfURY, dfLLX, dfLLR, dfLRX, dfLRY;// corners coordinates
double dfSizeX, dfSizY;
double dfPixSizeX = ((dfURX - dfULX) / 2 ) / dfSizeX; //assume dfURX - dfULX == dfLRX - dfLLX
double dfPixSizeY = ((dаLLY - dfULY) / 2 ) / dfSizeY; //assume dfURY - dfULY == dfLRY - dfLLY
double adfGeoTransform[6] = {dfULX, dfPixSizeX , 0, dfULY, 0, dfPixSizeY };

[pirateuse]

Дмитрий Барышников, спасибо! да, согласен по пунктам. Дельные уточнения. Таки не зря я все эти исследования провел, что-то теперь знаю)
Но надо все-таки остановиться на п.5.
1) Простите дурака, но не понял откуда деление на 2 берется: ((dfURX - dfULX) / 2 ) / dfSizeX;
Вроде надо взять интервал (dfURX - dfULX) и поделить на кол-во пикселей (dfSizeX). Почему еще на 2?
2) Из комментария //assume dfURX - dfULX == dfLRX - dfLLX напрашивается вывод, что система координат прямоугольная. Иными словами, если растр в проекции UTM, то координаты углов и других точек надо получить именно в этой проекции, т.е. единицы будут метры? Я же не смогу по этим формулам посчитать параметры, если координаты углов заданы в градусах (широта,долгота)?

[Дмитрий Барышников]

1) Да, это ошибка, на 2 делить не надо. Но зато нужно брать абсолютное значение от разности, а то размер может быть отрицательным.
2) Нет, такой вывод не напрашивается. Четко написано, предполагается что размер изображения в верхней части (upper) равен размеру в нижней части (low). Единица измерения не важна. Если размеры не совпадают - значит есть поворот - тогда лучше использовать GCP

[pirateuse]

Но зато нужно брать абсолютное значение от разности, а то размер может быть отрицательным.

Думаю, в данном случае не надо брать абсолютное значение, иначе GDALApplyGeoTransform() будет выдавать неверные значения, используя формулы:

Код: Выделить всё

//transformation between pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space
Xp = adfGeoTransform[0] + P*adfGeoTransform[1] + L*adfGeoTransform[2];
Yp = adfGeoTransform[3] + P*adfGeoTransform[4] + L*adfGeoTransform[5];

Где adfGeoTransform[6] = {dfULX, dfPixSizeX , 0, dfULY, 0, dfPixSizeY}, как было написано выше. Т.е. если потерять знак "-" у 1-го и/или 5-го коэффициента, то и выходные координаты будут противоположны истинным относительно координат углов.

И еще раз спасибо за помощь! Много вещей прояснилось.