использование GDAL на C/C++

Ответить
HelloKitty
Новоприбывший
Сообщения: 7
Зарегистрирован: 29 ноя 2016, 22:49
Репутация: 0

использование GDAL на C/C++

Сообщение HelloKitty » 20 июл 2017, 12:13

Добрый день!
Подскажите пожалуйста как правильно использовать привязку по контрольным точкам?
По 4 точкам - понятно, через GDALSetGeoTransform растр привязался. Но что, если точек будет много?
1. Создаю лист контрольных точек с координатами привязки.
2. Делаю GDALSetGCPs.

В итоге QGis опорных точек не видит.
Может что-то еще нужно сделать? Есть ли какой-то еще способ?

Александр Мурый
Гуру
Сообщения: 5059
Зарегистрирован: 26 сен 2009, 16:26
Статьи: 3
Проекты: 5/1
Репутация: 717
Ваше звание: званий не имею
Откуда: Москва

Re: использование GDAL на C/C++

Сообщение Александр Мурый » 20 июл 2017, 12:21

Приведите код, что ли. Или рассчитываете на коллективную телепатию?
И желательно ещё образец получившегося привязанного растра.
Редактор материалов, модератор форума

HelloKitty
Новоприбывший
Сообщения: 7
Зарегистрирован: 29 ноя 2016, 22:49
Репутация: 0

Re: использование GDAL на C/C++

Сообщение HelloKitty » 20 июл 2017, 13:42

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

    GDAL_GCP *pasGCPs = NULL;
    int GCP_count = 5;
    pasGCPs = (GDAL_GCP*)CPLRealloc(pasGCPs, sizeof(GDAL_GCP) * GCP_count);
    GDAL_GCP *zeroGCPs = NULL;
    GDALInitGCPs(GCP_count, pasGCPs );
    zeroGCPs = 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;

    pasGCPs+= sizeof(GDAL_GCP);
    
    //2st point

    //Unique identifier, often numeric
    pasGCPs->pszId = "1";
    //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 = PicHeight;
    //X position of GCP in georeferenced space
    pasGCPs->dfGCPX = 590000.000;
    //Y position of GCP in georeferenced space
    pasGCPs->dfGCPY = 4914000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs->dfGCPZ = 0;

    pasGCPs+= sizeof(GDAL_GCP);

    //3st point

    //Unique identifier, often numeric
    pasGCPs->pszId = "2";
    //Informational message or ""
    pasGCPs->pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs->dfGCPPixel = PicWidth;
    //Line (y) location of GCP on raster
    pasGCPs->dfGCPLine = 0;
    //X position of GCP in georeferenced space
    pasGCPs->dfGCPX = 609000.000;
    //Y position of GCP in georeferenced space
    pasGCPs->dfGCPY = 4928000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs->dfGCPZ = 0;

    pasGCPs+= sizeof(GDAL_GCP);

    //4st point

    //Unique identifier, often numeric
    pasGCPs->pszId = "3";
    //Informational message or ""
    pasGCPs->pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs->dfGCPPixel = PicWidth;
    //Line (y) location of GCP on raster
    pasGCPs->dfGCPLine = PicHeight;
    //X position of GCP in georeferenced space
    pasGCPs->dfGCPX = 609000.000;
    //Y position of GCP in georeferenced space
    pasGCPs->dfGCPY = 4914000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs->dfGCPZ = 0;

    pasGCPs+= sizeof(GDAL_GCP);

    //5st point center

    //Unique identifier, often numeric
    pasGCPs->pszId = "4";
    //Informational message or ""
    pasGCPs->pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs->dfGCPPixel = PicWidth/2;
    //Line (y) location of GCP on raster
    pasGCPs->dfGCPLine = PicHeight/2;
    //X position of GCP in georeferenced space
    pasGCPs->dfGCPX = 599500.000;
    //Y position of GCP in georeferenced space
    pasGCPs->dfGCPY = 4921000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs->dfGCPZ = 0;

    pasGCPs = zeroGCPs;


    GDALSetGCPs(dst, GCP_count,pasGCPs,"NAD27");
и растр пришлось сжать
Вложения
GeoDakota.zip
(4.75 МБ) 15 скачиваний

gamm
Гуру
Сообщения: 2588
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 542
Ваше звание: программист
Откуда: Казань

Re: использование GDAL на C/C++

Сообщение gamm » 20 июл 2017, 16:49

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

pasGCPs+= sizeof(GDAL_GCP);
это типа юмор такой? sizeof() вроде в байтах выдается, а += прибавляет в штуках (pasGCPs), программа должна просто упасть ...

HelloKitty
Новоприбывший
Сообщения: 7
Зарегистрирован: 29 ноя 2016, 22:49
Репутация: 0

Re: использование GDAL на C/C++

Сообщение HelloKitty » 21 июл 2017, 13:34

Да, извиняюсь :oops: Изначально было вот так: (потом началось экспериментирование и добавление новых косяков :) )

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

    int GCP_count = 5;
    GDAL_GCP *pasGCPs = NULL;
    pasGCPs = (GDAL_GCP*)CPLRealloc(pasGCPs, sizeof(GDAL_GCP) * GCP_count);
    GDALInitGCPs(GCP_count, pasGCPs );

    //1st point

    //Unique identifier, often numeric
    pasGCPs[0].pszId = "0";
    //Informational message or ""
    pasGCPs[0].pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs[0].dfGCPPixel = 0;
    //Line (y) location of GCP on raster
    pasGCPs[0].dfGCPLine = 0;
    //X position of GCP in georeferenced space
    pasGCPs[0].dfGCPX = 590000.000;
    //Y position of GCP in georeferenced space
    pasGCPs[0].dfGCPY = 4928000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs[0].dfGCPZ = 0;

    //pasGCPs +=  7 ;

    //2st point

    //Unique identifier, often numeric
    pasGCPs[1].pszId = "1";
    //Informational message or ""
    pasGCPs[1].pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs[1].dfGCPPixel = 0;
    //Line (y) location of GCP on raster
    pasGCPs[1].dfGCPLine = PicHeight;
    //X position of GCP in georeferenced space
    pasGCPs[1].dfGCPX = 590000.000;
    //Y position of GCP in georeferenced space
    pasGCPs[1].dfGCPY = 4914000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs[1].dfGCPZ = 0;

    //3st point

    //Unique identifier, often numeric
    pasGCPs[2].pszId = "2";
    //Informational message or ""
    pasGCPs[2].pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs[2].dfGCPPixel = PicWidth;
    //Line (y) location of GCP on raster
    pasGCPs[2].dfGCPLine = 0;
    //X position of GCP in georeferenced space
    pasGCPs[2].dfGCPX = 609000.000;
    //Y position of GCP in georeferenced space
    pasGCPs[2].dfGCPY = 4928000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs[2].dfGCPZ = 0;


    //4st point

    //Unique identifier, often numeric
    pasGCPs[3].pszId = "3";
    //Informational message or ""
    pasGCPs[3].pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs[3].dfGCPPixel = PicWidth;
    //Line (y) location of GCP on raster
    pasGCPs[3].dfGCPLine = PicHeight;
    //X position of GCP in georeferenced space
    pasGCPs[3].dfGCPX = 609000.000;
    //Y position of GCP in georeferenced space
    pasGCPs[3].dfGCPY = 4914000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs[3].dfGCPZ = 0;

    //5st point center

    //Unique identifier, often numeric
    pasGCPs[4].pszId = "4";
    //Informational message or ""
    pasGCPs[4].pszInfo = NULL;
    //Pixel (x) location of GCP on raster
    pasGCPs[4].dfGCPPixel = PicWidth/2;
    //Line (y) location of GCP on raster
    pasGCPs[4].dfGCPLine = PicHeight/2;
    //X position of GCP in georeferenced space
    pasGCPs[4].dfGCPX = 599500.000;
    //Y position of GCP in georeferenced space
    pasGCPs[4].dfGCPY = 4921000.000;
    //Elevation of GCP, or zero if not known
    pasGCPs[4].dfGCPZ = 0;



    GDALSetGCPs(dst, GCP_count,pasGCPs,"NAD27");
Результат такой же :(

Александр Мурый
Гуру
Сообщения: 5059
Зарегистрирован: 26 сен 2009, 16:26
Статьи: 3
Проекты: 5/1
Репутация: 717
Ваше звание: званий не имею
Откуда: Москва

Re: использование GDAL на C/C++

Сообщение Александр Мурый » 21 июл 2017, 14:04

Вы уверены, что правильно задан параметр pszGCPProjection для функции GDALSetGCPs?

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

GDALSetGCPs(dst, GCP_count,pasGCPs,"NAD27");
Описание pszGCPProjection из доков:
pszGCPProjection the new OGC WKT coordinate system to assign for the GCP output coordinates. This parameter should be "" if no output coordinate system is known.
Указанный у вас NAD27 - датум, т.е. только часть описания системы координат. Гугл подсказывает, что там можно задавать строку в формате PROJ типа "+proj=latlong +datum=WGS84".
Редактор материалов, модератор форума

HelloKitty
Новоприбывший
Сообщения: 7
Зарегистрирован: 29 ноя 2016, 22:49
Репутация: 0

Re: использование GDAL на C/C++

Сообщение HelloKitty » 21 июл 2017, 14:40

Александр Мурый писал(а): Указанный у вас NAD27 - датум, т.е. только часть описания системы координат. Гугл подсказывает, что там можно задавать строку в формате PROJ типа "+proj=latlong +datum=WGS84".
Возможно, но, как я понимаю, этот влияет только на то, распознает ли автоматически qgis, в какую проекцию трансформировать. Параметр можно вообще оставить пустым, тогда при добавлении растрового слоя в qgis просто нужно указать, какая это проекция. Поправьте, если не так.
Со строкой типа "+proj=latlong +datum=NAD27" все так же.

Александр Мурый
Гуру
Сообщения: 5059
Зарегистрирован: 26 сен 2009, 16:26
Статьи: 3
Проекты: 5/1
Репутация: 717
Ваше звание: званий не имею
Откуда: Москва

Re: использование GDAL на C/C++

Сообщение Александр Мурый » 21 июл 2017, 14:58

В свойствах растра вот такое:

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

Size is 9500, 7000
Coordinate System is:
LOCAL_CS["unnamed",
    UNIT["unknown",1]]
Origin = (-1476.500000000000000,1476.500000000000000)
Pixel Size = (2953.000000000000000,-2953.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   -1476.500,    1476.500) 
Lower Left  (   -1476.500,-20669523.500) 
Upper Right (28052023.500,    1476.500) 
Lower Right (28052023.500,-20669523.500) 
Center      (14025273.500,-10334023.500)
По моим ощущениям, нормальный привязанный растр размером 9500х7000 не может иметь таких координат углов. Откуда вообще брались координаты для привязки и в какой именно системе координат?

А что, если привязывать по 4-м точкам, всё нормально проходит?
Редактор материалов, модератор форума

HelloKitty
Новоприбывший
Сообщения: 7
Зарегистрирован: 29 ноя 2016, 22:49
Репутация: 0

Re: использование GDAL на C/C++

Сообщение HelloKitty » 21 июл 2017, 15:30

Александр Мурый писал(а): По моим ощущениям, нормальный привязанный растр размером 9500х7000 не может иметь таких координат углов. Откуда вообще брались координаты для привязки и в какой именно системе координат?

А что, если привязывать по 4-м точкам, всё нормально проходит?
Система координат NAD27, это Северная Дакота
по 4 точкам:

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

    double BottomLeftX  = 590000.000;
    double BottomLeftY  = 4914000.000;
    double TopLeftX     = 590000.000;
    double TopLeftY     = 4928000.000;

    double BottomRightX = 609000.000;
    double BottomRightY = 4914000.000;
    double TopRightX    = 609000.000;
    double TopRightY    = 4928000.000;

    double MapWidth  = abs( TopRightX - BottomLeftX  );
    double MapHeight = abs( TopRightY - BottomLeftY );
    double WidthResolution = MapWidth/PicWidth;
    double HeightResolution =  (-1) * MapHeight/PicHeight;

    double adfGeoTransform[6] = { TopLeftX, WidthResolution, 0, TopLeftY, 0, HeightResolution };
    GDALSetGeoTransform( dst, adfGeoTransform);
получается вот такой растр.
Координаты были взяты из скачанного примера файла с привязкой
Вложения
GEO_Dakota_4_points.zip
(4.75 МБ) 15 скачиваний

HelloKitty
Новоприбывший
Сообщения: 7
Зарегистрирован: 29 ноя 2016, 22:49
Репутация: 0

Re: использование GDAL на C/C++

Сообщение HelloKitty » 21 июл 2017, 15:43

Если что, то вот он, этот примерчик
Вложения
spearfish_toposheet.zip
(4.3 МБ) 16 скачиваний

Serpuh
Новоприбывший
Сообщения: 11
Зарегистрирован: 27 авг 2017, 14:46
Репутация: 0
Откуда: Москва

Re: использование GDAL на C/C++

Сообщение Serpuh » 28 авг 2017, 09:46

HelloKitty писал(а):
20 июл 2017, 12:13
По 4 точкам - понятно, через GDALSetGeoTransform растр привязался. Но что, если точек будет много?
Подскажите, а в чем смысл множества контрольных точек? Только в их высоте?

Ответить

Вернуться в «GDAL/OGR»