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