Получить строка/столбец из Ш/Д

[wowka1319]

Есть GDALDataset. Как получить координаты в пикселях (причем в вещественных числах) из широты/долготы? Как обратно сделать я знаю, но вот так - нет. Кто-то мне давал ссылку на функцию int GDALInvGeoTransform(double * gt_in, double * gt_out), но это явно не оно.

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

Прежде чем заводить тему стоит посмотреть, нет ли похожих.
viewtopic.php?f=30&t=16017&p=108576

[wowka1319]

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

1) Для преобразования из географических координат в пиксельные OGRSpatialReference не подходит. Он вообще не может ничего преобразовывать. Скорее OGRCoordinateTransofrmation - да и то, только для преобразования между различными системами координат.

которое ничем помочь конечно же не может.
Поэтому вопрос остался.

[wowka1319]

Ну как же черт возьми получить X/Y в пикселях по широте/долготе???? Меня уже несколько недель мучает эта проблема. Отчаянно прошу помощи.... И нету на этом форуме темы с решением!!! Дайте простой код получения нужной информации из GDALDataSetH!!!!!!!!!!!!!!

[wowka1319]

up. Вопрос встал остро! хелпми!!!

[gamm]

up. Вопрос встал остро! хелпми!!!

код за вас никто писать, конечно, не будет, но схема действий такая (если нет извратов типа повернутого изображения)

1) достаете строку +proj и переводите долготу/широту в координаты проекции снимка
2) достаете bbox (охватывающий прямоугольник) Xmin,Ymin,Xmsx,Ymax, вычисляете dx,dy (я не уверен, но вроде bbox считается по внешним краям пикселей), dx=(Xmax-Xmin)/Nx, Nx x Ny - размеры снимка, и вроде dx, dy там готовые есть
3) вычисляете номер пикселя Xpix=(X-Xmin)/dx, нумерация с 0
4) profit

P.S. я все это делаю в R, там обращение в gdal через обертку, но принцип тот же.

[wowka1319]

1) достаете строку +proj и переводите долготу/широту в координаты проекции снимка

как? обратный процесс мне пол силу, но вот так - нет.
А дальше вроде бы понятно, но все-таки опишу чтоб вы если что поправили:
2) вычисляем размер пикселя (dx и dy) в условных единицах проекции снимка по X и Y (они разве не равны?)
3) отталкиваясь от Xmin и Ymax (не очень ясно понимаю где их брать, но думаю разберусь) считаю (X-Xmin)/dx и (Y-Ymin)/dy
4)
верно? ну прошу разъяснить первый пункт - именно в нем загвоздка!!!

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

Из пикселей в координаты

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

GDALDataset *pDataset = (GDALDataset *)GDALOpen("test.tiff", GA_ReadOnly);
double adfGeoTransform[6] = {0, 0, 0, 0, 0, 0};
pDataset->GetGeoTransform(adfGeoTransform);

double dfRow = 1.5, dfCol = 2.3;
double dfX, dfY;

GDALApplyGeoTransform(adfInvGeoTransform, dfRow, dfCol, &dfX, &dfY);

Из координат в пиксели

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

GDALDataset *pDataset = (GDALDataset *)GDALOpen("test.tiff", GA_ReadOnly);
double adfGeoTransform[6] = {0, 0, 0, 0, 0, 0};
pDataset->GetGeoTransform(adfGeoTransform);
double adfInvGeoTransform[6] = {0, 0, 0, 0, 0, 0};
GDALInvGeoTransform(adfGeoTransform, adfInvGeoTransform);

double dfRow, dfCol;
double dfX = 37.55, dfY = 55.37;

GDALApplyGeoTransform(adfInvGeoTransform, dfX, dfY, &dfRow, &dfCol);

Если координаты снимка отличаются от координат на входе - их надо перепроецировать в координаты снимка при помощи OGRCoordinateTransformation

[wowka1319]

увы но ваш код работает неверно:
Во-первых:
Из пикселей в координаты будет так: (это по мануалу GDAL)

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

    double x,y;
    double geoTransform[6];
    dataset_->GetGeoTransform(geoTransform);
    GDALApplyGeoTransform(geoTransform, x, y, &x, &y); //использую те же переменные
    const char* projection = dataset_->GetProjectionRef();
    OGRSpatialReference* spatialReference = new OGRSpatialReference(projection);
    if (spatialReference->IsProjected())
    {
        OGRSpatialReference* spatialReferenceGeogCS = spatialReference->CloneGeogCS();
        OGRCoordinateTransformation* coordinateTransformation = OGRCreateCoordinateTransformation(spatialReference, spatialReferenceGeogCS);
        if (coordinateTransformation) coordinateTransformation->Transform(1, &x, &y); //опять использую те же переменные
    }

    // теперь в x - долгота, а в y - широта в градусах!!!
    // проверено - РАБОТАЕТ!!!

не понимаю зачем вы вообще приводили первый код, ведь этого не требовалось

Во-вторых: ваш пример Из координат в пиксели работает НЕВЕРНО.
может вы не поняли задачу? повторю - из координат долгота/широта в градусах получить координаты изображения X/Y в пикселях, т.е. нужен код работающий ровно наоборот коду, приведенному выше (в этом сообщении).
Вопрос остался. Пожалуйста помогите - проблема серьезно встала, не дает работать дальше!!!!!!!

[gamm]

увы но ваш код работает неверно

работает верно, вы просто не поняли, что там написано - это код для перевода пикселей в координаты проекции тифа и обратно линейным преобразованием.

Если проекция долгота/широта, то получатся искомые пиксели. Если нет, то сначала нужно преобразовать градусы в проекцию тифа. У вас все части есть, в чем вопрос - непонятно.

[wowka1319]

Вот именно, что мне нужен был код (сотый раз повторяю) из долготы/широты (НЕ из проекции тифа в метрах, а именно градусов) в пиксели (наоборот не надо - уже давноооо есть). Так никто рабочий код и не привел(((
Но взглянув на обратный код (который я привел постом выше), я по аналогии все-таки добился своего.
Итак, код для тех возможно сталкнется с этой проблемой (думаю Qt типы никого не собьют с толку):

Из пикселей X/Y в градусы Д/Ш (то, что уже приводил постом выше)

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

QPointF MMGeoImage::geoCoords(const QPointF &pixelCoords)
{
    if (!dataset_) return QPointF();

    double point[2];
    point[0] = pixelCoords.x();
    point[1] = pixelCoords.y();

    double geoTransform[6];
    dataset_->GetGeoTransform(geoTransform);
    GDALApplyGeoTransform(geoTransform, point[0], point[1], &point[0], &point[1]);

    const char* projection = dataset_->GetProjectionRef();
    OGRSpatialReference* spatialReference = new OGRSpatialReference(projection);
    if (spatialReference->IsProjected())
    {
        OGRSpatialReference* spatialReferenceGeogCS = spatialReference->CloneGeogCS();
        OGRCoordinateTransformation* coordinateTransformation = OGRCreateCoordinateTransformation(spatialReference, spatialReferenceGeogCS);
        if (coordinateTransformation)
        {
            coordinateTransformation->Transform(1, &point[0], &point[1]);
            return QPointF(point[0], point[1]); //OK
        }
        else return QPointF();
    }
    else return QPointF();
}

dataset_ - это GDALDataset*

И, наконец, искомый код - Из градус Д/Ш в пиксели X/Y:

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

QPointF MMGeoImage::pixelCoords(const QPointF &geoCoords)
{
    if (!dataset_) return QPointF();

    double point[2];
    point[0] = geoCoords.x();
    point[1] = geoCoords.y();
    const char* projection = dataset_->GetProjectionRef();
    OGRSpatialReference* spatialReference = new OGRSpatialReference(projection);
    if (spatialReference->IsProjected())
    {
        OGRSpatialReference* spatialReferenceGeogCS = spatialReference->CloneGeogCS();
        OGRCoordinateTransformation* coordinateTransformation = OGRCreateCoordinateTransformation(spatialReferenceGeogCS, spatialReference);
        if (coordinateTransformation)
        {
            coordinateTransformation->Transform(1, &point[0], &point[1]);

            double geoTransform[6];
            dataset_->GetGeoTransform(geoTransform);
            double invGeoTransform[6];
            GDALInvGeoTransform(geoTransform, invGeoTransform);

            GDALApplyGeoTransform(invGeoTransform, point[0], point[1], &point[0], &point[1]);
            return QPointF(point[0], point[1]); //OK
        }
        else return QPointF();
    }
    else return QPointF();
}

dataset_ - это GDALDataset*