Re[3]: GDAL API C++ - C/C++. Прикладные вопросы

Здравствуйте, Serpuh, Вы писали:

S>А есть на форуме кто работает с этой картографической библиотекой? Нужно преобразовать tiff в geoTiff с известными GCP. Примеры есть урывками и полную pipeline закодить пока не получается.

GeoTiff — это "привязанный" TIFF с преобразованием и координатной системой?
http://www.gdal.org/warptut.html — вроде тут все есть.
Для исходного Tiff-а вызвать SetGCPs, потом создать новый Tiff и сделать Warp из исходного в новый.

Здравствуйте, nrwl, Вы писали:
N>GeoTiff — это "привязанный" TIFF с преобразованием и координатной системой?
Да, можно в QGIS закачивать и работать с ним, ортофото делаю из 3D модели полученную из фоток с дрона
3d моделька

	Скрытый текст

ортофото с нее

	Скрытый текст

Яндекс https://yandex.ru/maps/?ll=37.064733%2C55.784008&z=18&l=sat%2Cskl

Еще одна моделька

	Скрытый текст

N>http://www.gdal.org/warptut.html — вроде тут все есть.
Я сделал через CreateCopy(), но дьявол в деталях. Как указывать правильно проекции (закоментаренные geographicRef), как метры перевести в UTM, там вроде один в один, но есть некий масштабный коэфициент 0.9996. Поиском просто так не ищется, у GDAL help очень краткий.

string  pszFilename("Scene_dense_mesh_texture_orthomap.tif");
GDALDataset  *poSrcDS;
GDALAllRegister();
poSrcDS = (GDALDataset *)GDALOpen(pszFilename.c_str(), GA_ReadOnly);

.....

double adfGeoTransform[6] = { TopLeftX, WidthResolution, 0, TopLeftY, 0, HeightResolution };

const char *pszFormat = "GTiff";
GDALDriver *poDriver;
GDALDataset *dstDataset;
GDALAllRegister();
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
dstDataset = poDriver->CreateCopy("Scene_dense_mesh_texture_orthomap_gdal2.tif", poSrcDS, false, NULL, NULL, NULL);


// specify a equirectangular projection (also known as geographic or uprojected lat/lon) for the VRT
OGRSpatialReference geographicRef;
geographicRef.SetWellKnownGeogCS("WGS84"); // standard WGS84 coordinate system
geographicRef.SetUTM( 37, TRUE );

//geographicRef.SetProjection(SRS_PT_EQUIRECTANGULAR);
//geographicRef.SetProjCS("Equirectangular / WGS84");
// geographicRef.SetEquirectangular(0,0,0,0); // gives same results as above
//geographicRef.SetAngularUnits(SRS_UA_DEGREE,0.0174532925199433); // has no effect
//geographicRef.SetAngularUnits(SRS_UA_RADIAN,1.0); // has no effect 
char *pszDstWkt;

geographicRef.exportToWkt(&pszDstWkt);
dstDataset->SetProjection(pszDstWkt);
cout << pszDstWkt;

GDALSetGeoTransform(dstDataset, adfGeoTransform);

if (dstDataset != NULL) GDALClose((GDALDatasetH)dstDataset);        
if (poSrcDS != NULL) GDALClose((GDALDatasetH)dstDataset);

Здравствуйте, Serpuh, Вы писали:

S>Как указывать правильно проекции (закоментаренные geographicRef)

А что именно получить хочется? Если не надо какую-то специфическую координатную систему, то проще взять или обычный WGS84, или гугловую проекцию (EPSG:3857)
Ну и если GCP-щные точки есть, то должна бы быть и координатная система для них.

OGRSpatialReference ref;
ref.importFromEpsg(3857);
char *pszDstWkt;
ref.exportToWkt(&pszDstWkt);
dataset->SetProjection(pszDstWkt);

S>как метры перевести в UTM, там вроде один в один, но есть некий масштабный коэффициент 0.9996

UTM и так в метрах. И что за масштабный коэффициент?

S>double adfGeoTransform[6] = { TopLeftX, WidthResolution, 0, TopLeftY, 0, HeightResolution };

А параметры для SetGeoTransform — они откуда взялись?

Здравствуйте, nrwl, Вы писали:
N>А что именно получить хочется? Если не надо какую-то специфическую координатную систему, то проще взять или обычный WGS84, или гугловую проекцию (EPSG:3857)
Я пока так указал, geoTiff получился и размеры от экранной линейки в QGIS вроде правильные.

geographicRef.SetWellKnownGeogCS("WGS84");
geographicRef.SetUTM(37, true);

N>UTM и так в метрах. И что за масштабный коэффициент?
PARAMETER["scale_factor",0.9996] gdalinfo.exe выдает эту инфу. Подумал может надо вручную размеры в метрах на это умножать.

N>А параметры для SetGeoTransform — они откуда взялись?
3D модель привязана по EXIF GPS из снимков, развернута по Север-Юг и сверху сделано ортофото т.е. фактически для каждого пикселя известна его координата и высота. Но я пока указываю в SetGeoTransform UTM верхней точки и размер ортофото на местности.
Не в курсе, вроде можно по XYZ сетке привязывать, чтобы geoTiff был с картой высот?

Здравствуйте, Serpuh, Вы писали:
S>и размеры от экранной линейки в QGIS вроде правильные

Я бы для проверки взял несколько точек с GPS координатами и посмотрел где их QGis нарисует — если они на растре там где надо, то все нормально "привязалось".

S>А вот с зоной UTM непонятно, можно ли средствами GDAL перевести latitude longtitude в UTM зону типа такого? http://www.apsalin.com/utm-zone-finder.aspx

Насколько я знаю, готового в GDAL — нет. Только если самому написать. В википедии грид с зонами есть — https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
У UTM-а обычно зона прописана и "на лету" ее никто не меняет.

S>PARAMETER["scale_factor",0.9996] gdalinfo.exe выдает эту инфу. Подумал может надо вручную размеры в метрах на это умножать.

Этот параметр ко все координатной системе относится. Если нужно больше деталей — нужно раскуривать детали преобразования в UTM.
Для единиц там другой параметр — UNIT["metre",1,...]. Но все преобразования из системы в систему должен сделать GDAL/proj4.

S>3D модель привязана по EXIF GPS из снимков, развернута по Север-Юг и сверху сделано ортофото т.е. фактически для каждого пикселя известна его координата и высота. Но я пока указываю в SetGeoTransform UTM верхней точки и размер ортофото на местности.
S>Не в курсе, вроде можно по XYZ сетке привязывать, чтобы geoTiff был с картой высот?

Насколько я знаю, GDAL с XYZ не работает для растра. Привязка — двумерная, там и параметров-то нет для третьей координаты

Здравствуйте, nrwl, Вы писали:
N>Насколько я знаю, GDAL с XYZ не работает для растра. Привязка — двумерная, там и параметров-то нет для третьей координаты

Такой код нашел, можно через GCP высоту задать. Но если задать к примеру 1млн. GCP не грохнется ли алгоритм.

    GDAL_GCP *pasGCPs;
    ...
//1st point
    //Unique identifier, often numeric
    pasGCPs->pszId = "0";
    //Informational message or ""
    pasGCPs->pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs->dfGCPPixel = 0;
    //Line (y) location of GCP on raster
    pasGCPs->dfGCPLine = 0;
    //X position of GCP in georeferenced space
    pasGCPs->dfGCPX = 590000.000;
    //Y position of GCP in georeferenced space
    pasGCPs->dfGCPY = 4928000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs->dfGCPZ = 0;
    ...
    GDALSetGCPs(dst, GCP_count,pasGCPs,"WGS84");

Здравствуйте, Serpuh, Вы писали:

S>Здравствуйте, nrwl, Вы писали:
N>>Насколько я знаю, GDAL с XYZ не работает для растра. Привязка — двумерная, там и параметров-то нет для третьей координаты

S>Такой код нашел, можно через GCP высоту задать. Но если задать к примеру 1млн. GCP не грохнется ли алгоритм.

https://en.wikipedia.org/wiki/World_file — GDAL в результате вот такую матрицу сделает. Там третью координату прицепить-то некуда. Можно посмотреть что именно они с третьей координатой из GCP делают.

От:	Serpuh	фотомер.рф
Дата:	28.08.17 10:08
Оценка:

	От:	nrwl
	Дата:	29.08.17 15:53
	Оценка:

От:	Serpuh	фотомер.рф
Дата:	30.08.17 07:57
Оценка:

	От:	nrwl
	Дата:	30.08.17 23:26
	Оценка:

От:	Serpuh	фотомер.рф
Дата:	31.08.17 07:47
Оценка: