GDAL: Как получить ГК GeoTiff

[wowka1319]

Необходимо получать ГК пикселей (что бы потом по координатам накладывать на снимок названия объектов).
Пока остановился на

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

GDALRegister_GTiff();
	double adfGeoTransform[6];
	GDALDataset  *poDataset;

    poDataset = (GDALDataset *) GDALOpen("TEST_01.tiff", GA_ReadOnly);
    if (poDataset == NULL)
    {
		std::cout << "WTF?\n";
    }

Другими словами дальше туториала не ушел.
что дальше?

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

Приложите вывод gdalinfo для TEST_01.tiff в тему

[wowka1319]

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

Driver: GTiff/GeoTIFF
Files: TEST_01.tiff
Size is 4483, 2494
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Pulkovo 1942",
        DATUM["Pulkovo_1942",
            SPHEROID["Krassowsky 1940",6378245,298.2999999999985,
                AUTHORITY["EPSG","7024"]],
            AUTHORITY["EPSG","6284"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4284"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",111],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (73671.999948490688000,5838215.177363829700000)
Pixel Size = (50.000000000000057,-50.040096230954290)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   73672.000, 5838215.177) (112d 5'20.30"E, 52d40'0.97"N)
Lower Left  (   73672.000, 5713415.177) (112d 3'43.17"E, 51d32'43.98"N)
Upper Right (  297822.000, 5838215.177) (115d23'44.13"E, 52d35'25.17"N)
Lower Right (  297822.000, 5713415.177) (115d17'13.46"E, 51d28'19.01"N)
Center      (  185747.000, 5775815.177) (113d42'34.71"E, 52d 4'48.09"N)
Band 1 Block=4483x1 Type=UInt16, ColorInterp=Red
Band 2 Block=4483x1 Type=UInt16, ColorInterp=Green
Band 3 Block=4483x1 Type=UInt16, ColorInterp=Blue

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

Потом так:
1. Получаем матрицу преобразования
2. Собственно пересчитываем пикселы в СК изображения.

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

double        adfGeoTransform[6];
if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None ) //1
{
  double dfPixX(0), dfPixY(0);  //пиксельные координаты, например 0,0
  double dfSRSX, dfSRSY;        //координаты в системе координат изображения
  GDALApplyGeoTransform( adfGeoTransform, dfPixX, dfPixY, &dfSRSX, &dfSRSY ); //2
}
 

Полезные ссылки: http://www.gdal.org/gdal_tutorial.html

[wowka1319]

Спасибо. Я получил прямоугольные координаты. Как теперь их перевести в ГК? gdalinfo я указывал выше.

[Максим Дубинин]

координаты в спроецированной СК, например ГК (Гаусс-Крюгер) это и есть прямоугольные координаты.

[wowka1319]

но мне нужны Географические Координаты (Широта и долгота)! Как их получить?

[trir]

EPSG=4284

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

Как-то так:

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

        const char* pszProjection = poDataset ->GetProjectionRef();
        if(CPLStrnlen(pszProjection, 10) > 0)
        {
            OGRSpatialReference* pRasterSpatRef = new OGRSpatialReference(pszProjection);
            if(pRasterSpatRef->IsProjected())
            {
                OGRSpatialReference* pGSSpatRef = pRasterSpatRef->CloneGeogCS();
                OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation( pRasterSpatRef, pGSSpatRef );
                if(poCT)
                {
                    poCT->Transform(1, &dfSRSX, &dfSRSY);
                }
            }
        }