Задание привязки растра через GDAL для использования в Mapinfo

Системы координат, проекции, преобразования, привязка
spawner
Активный участник
Сообщения: 103
Зарегистрирован: 09 окт 2009, 16:49
Репутация: 10

Задание привязки растра через GDAL для использования в Mapinfo

Сообщение spawner » 11 авг 2020, 15:43

Пытаюсь разобраться с GDAL, для использования в скриптах, возникла проблема:

Есть векторный файл в формате tab, система координат МСК-45, зона 2, в tab файле система координат не установлена (план-схема)

Я пытаюсь через GDAL назначить для tab файла правильную систему координат (параметры которой я нашел где-то тут на форуме в формате в PROJ.4). Для пересчета использую строку вида:

ogr2ogr -f "MapInfo File" -a_srs "+proj=tmerc +lat_0=0 +lon_0=64.03333333333 +k=1 +x_0=2300000 +y_0=-5709414.70 +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs" out_file.tab in_file.tab

В конечном файле получается система координат:
CoordSys Earth Projection 8, 104, "m", 64.03333333333, 0, 1, 2300000, -5709414.7000000002

а должна быть:
CoordSys Earth Projection 8, 1001, "m", 64.03333333333, 0, 1, 2300000, -5709414.6900000004

Т.е. вместо проекции 1001 получается 104, и графика отлетает на 100 метров.

Подскажите, в чем ошибка?

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение trir » 11 авг 2020, 15:52

попробуй +towgs84=24,-123,-94 -0.02,0.25,0.13,1.1

что за скрипты?

Аватара пользователя
Игорь Белов
Гуру
Сообщения: 2229
Зарегистрирован: 04 янв 2011, 22:00
Репутация: 1501
Откуда: Казань

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение Игорь Белов » 11 авг 2020, 17:34

Забудьте про формат PROJ.4, в данном случае он больше не работает.
Про параметр TOWGS84 тем более забудьте, трансформация теперь не является атрибутом проекции.
GDAL вслед за PROJ перешёл на формат описания систем координат WKT.

Для пересчета использем строку вида:

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

ogr2ogr -f "MapInfo File" -a_srs msk45z2.wkt out_file.tab in_file.tab
Файл msk45z2.wkt содержит описание системы координат «МСК-45 зона 2»:

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

PROJCRS["Pulkovo 1942 / MSK45 zone 2",
    BASEGEOGCRS["Pulkovo 1942",
        DATUM["Pulkovo 1942",
            ELLIPSOID["Krassowsky 1940",6378245,298.3,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4284]],
    CONVERSION["MSK45 zone 2",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",64.0333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",2300000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",-5709414.7,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Russia - MSK45 zone 2"],
        BBOX[54.1,62.5,56.8,65.5]]]
P. S. Поправил BBOX для правильного отображения в QGIS.
Вложения
msk45z2.zip
(604 байт) 297 скачиваний
The purpose of computing is insight, not numbers

spawner
Активный участник
Сообщения: 103
Зарегистрирован: 09 окт 2009, 16:49
Репутация: 10

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение spawner » 11 авг 2020, 19:26

Оооо, спасибо огромное! Все заработало!

Только почему-то мапинфо отказывается показывать установленную систему координат, пишет просто "Широта/Долгота".

В результате пересчета файла получилось:
CoordSys Earth Projection 8, 1001, "m", 64.033333333333, 0, 1, 2300000, -5709414.6900000004 Bounds (-27700000, -20709414.69) (32300000, 9290585.31)

А в файле с установленной в mapinfo системой координат:
CoordSys Earth Projection 8, 1001, "m", 64.03333333333, 0, 1, 2300000, -5709414.6900000004 Bounds (-5949281.53901, -15711552.1878) (10549281.539, 4292722.80776)

Т.е. вроде все совпадает, кроме пределов. Не в курсе, где они в мапинфо настраиваются? В mapfinwo.prj вроде нет такой настройки.

spawner
Активный участник
Сообщения: 103
Зарегистрирован: 09 окт 2009, 16:49
Репутация: 10

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение spawner » 11 авг 2020, 19:28

trir писал(а):
11 авг 2020, 15:52
что за скрипты?
Хочу сделать, чтоб из SasPlanet снимок сразу в нужную систему координат пересчитывался.
Делать это руками не все могут обучиться :/ да и кликать надоело по кнопкам.

Аватара пользователя
Игорь Белов
Гуру
Сообщения: 2229
Зарегистрирован: 04 янв 2011, 22:00
Репутация: 1501
Откуда: Казань

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение Игорь Белов » 11 авг 2020, 19:44

spawner писал(а):
11 авг 2020, 19:26
почему-то мапинфо отказывается показывать установленную систему координат, пишет просто "Широта/Долгота"
Различаются:
  • количество знаков после запятой в числе 64.033333333333;
  • пределы Bounds.
Для MapInfo это основание считать системы координат разными. Если система координат таблицы отсутствует в файле MAPINFOW.PRJ, MapInfo показывает "Широта/Долгота".
The purpose of computing is insight, not numbers

spawner
Активный участник
Сообщения: 103
Зарегистрирован: 09 окт 2009, 16:49
Репутация: 10

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение spawner » 11 авг 2020, 20:20

:oops: да, что-то знаки после запятой я плохо посчитал. Теперь мапинфо все распознает. Спасибо!

Аватара пользователя
dab
Гуру
Сообщения: 671
Зарегистрирован: 16 дек 2011, 20:02
Репутация: 170
Ваше звание: Гуру
Откуда: Москва
Контактная информация:

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение dab » 11 авг 2020, 21:10

Игорь, а уточните, пожалуйста по параметру
BBOX[54.1,62.5,56.8,65.5]]]
При
PARAMETER["Longitude of natural origin",64.0333333333333,

А почему ... 62.5 ... 65.5 ... а не ... 62.5333333333333 ... 65.5333333333333 ... ?

Аватара пользователя
Игорь Белов
Гуру
Сообщения: 2229
Зарегистрирован: 04 янв 2011, 22:00
Репутация: 1501
Откуда: Казань

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение Игорь Белов » 12 авг 2020, 06:28

dab, BBOX нужен при выборе системы координат для отображения области покрытия на карте. Пределы даются в WGS 84, не в Пулково. Сотая доля градуса соответствует километру на местности, достаточно для описания области действия системы координат малого населённого пункта.

Пример: EPSG:2935 "Pulkovo 1942 / CS63 zone A1"

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

PROJCRS["Pulkovo 1942 / CS63 zone A1",
    BASEGEOGCRS["Pulkovo 1942",
        DATUM["Pulkovo 1942",
            ELLIPSOID["Krassowsky 1940",6378245,298.3,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4284]],
    CONVERSION["CS63 zone A1",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0.116666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",41.5333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",1300000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Asia - FSU - CS63 zone A1"],
        BBOX[41.37,39.99,43.59,43.04]],
    ID["EPSG",2935]]
The purpose of computing is insight, not numbers

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

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение gamm » 12 авг 2020, 10:25

Игорь Белов писал(а):
11 авг 2020, 17:34
трансформация теперь не является атрибутом проекции.
но в WKT ее можно как-то указать? или можно использовать только прошитые в EPSG глобусы?

Аватара пользователя
Игорь Белов
Гуру
Сообщения: 2229
Зарегистрирован: 04 янв 2011, 22:00
Репутация: 1501
Откуда: Казань

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение Игорь Белов » 12 авг 2020, 13:51

gamm писал(а):
12 авг 2020, 10:25
Игорь Белов писал(а):
11 авг 2020, 17:34
трансформация теперь не является атрибутом проекции
но в WKT ее можно как-то указать? или можно использовать только прошитые в EPSG глобусы?
  • Если в рамках темы мечтать, что GDAL будет писать в таблицы MapInfo пользовательский датум, начинающийся с девяток, то нет, GDAL этим не занимается.
  • Понятно, что явно прописывать параметры трансформации, ассоциированной с датумом в EPSG, нет смысла. Софт позволяет выбирать трансформацию для датума из БД.
  • Ответ на вопрос о возможности произвольной трансформации ДА.
Рассмотрим пример.
Дано. В MapInfo имеются таблицы в проекции "GK zone 12 (Pulkovo 1942)", например.
Найти. Требуется создать СК такую, чтобы:
  1. объекты из таблиц MapInfo отображались в QGIS там, где их отображает MapInfo;
  2. при экспорте слоёв в таблицы MapInfo правильно писался датум #1001.
Решение.
Прежде всего отметим, что условие первого пункта можно было бы достигнуть выбором датума "Pulkovo 1942(83)", ведь именно он в восточногерманской версии используется в MapInfo под номером #1001. Но при этом не выполняется второй пункт, поскольку в EPSG "Pulkovo 1942(83)" и "Pulkovo 1942" разные датумы, и GDAL будет писать в таблицы MapInfo #104, что неправильно.

Решение состоит в том, чтобы создать связанную систему координат на основе "Pulkovo 1942" с конформным трёхмерным преобразованием, значения параметров которого позаимствовать из EPSG:1675 "Pulkovo 1942(83) to WGS 84 (1)". Назовём эту систему координат, скажем, "MI1001 GK zone 12":

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

BOUNDCRS[
    SOURCECRS[
        PROJCRS["MI1001 GK zone 12",
            BASEGEOGCRS["Pulkovo 1942",
                DATUM["Pulkovo 1942",
                    ELLIPSOID["Krassowsky 1942",6378245,298.3,
                        LENGTHUNIT["metre",1,
                            ID["EPSG",9001]]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8901]]],
            CONVERSION["unknown",
                METHOD["Transverse Mercator",
                    ID["EPSG",9807]],
                PARAMETER["Latitude of natural origin",0,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",69,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",1,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["False easting",12500000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]]],
            CS[Cartesian,2],
                AXIS["(E)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]],
                AXIS["(N)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from MI1001 to WGS84",
        METHOD["Coordinate Frame rotation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",24,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",-123,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",-94,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",-0.02,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0.25,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0.13,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1.0000011,
            ID["EPSG",8611]]]]
The purpose of computing is insight, not numbers

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение trir » 12 авг 2020, 14:01

по мне - так проще написать пару строчек на C# и решить задачу в лоб

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

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение gamm » 12 авг 2020, 15:42

trir писал(а):
12 авг 2020, 14:01
проще
я раньше тоже так думал. Но потом понял, что у штатных средств есть свои достоинства, и "пусть цветут сто цветов" 😎

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение trir » 12 авг 2020, 15:50

возможно тут можно что то сделать через MITAB_BOUNDS_FILE

Аватара пользователя
Игорь Белов
Гуру
Сообщения: 2229
Зарегистрирован: 04 янв 2011, 22:00
Репутация: 1501
Откуда: Казань

Re: Задание привязки растра через GDAL для использования в Mapinfo

Сообщение Игорь Белов » 12 авг 2020, 19:02

spawner писал(а):
11 авг 2020, 19:28
Хочу сделать, чтоб из SasPlanet снимок сразу в нужную систему координат пересчитывался.
Делать это руками не все могут обучиться :/ да и кликать надоело по кнопкам.
spawner писал(а):
11 авг 2020, 15:43
Есть векторный файл в формате tab…
Я пытаюсь через GDAL назначить для tab файла правильную систему координат
Что-то не пойму я постановку задачи.

Допустим, не устраивает точность привязки растров из Интернета, и векторные файлы содержат опознаки, к которым растры будут подтягиваться. Ну так геопривязка растров в этом случае — большая ручная работа с каждой картинкой индивидуально, какая уж тут автоматизация.

Значит, надо полагать, привязка растров устраивает? Но для чего тогда нужны какие-либо векторные файлы?

Склеим в SASPlanet в формат GeoTIFF растр topo1.tif. Трансформируем в МСК-45 зона 2:

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

gdalwarp -t_srs msk45z2.wkt topo1.tif topo2.tif
Откроем topo2.tif в MapInfo. В нормальной версии будет распознана внутренняя привязка. Появится диалог выбора системы координат, укажем "МСК45 зона 2". Всё!

Лень каждый раз отвечать на вопрос о системе координат? OK, добавим параметр:

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

gdalwarp -t_srs msk45z2.wkt -co TFW=YES topo1.tif topo2.tif
Вместе с растром topo2.tif создастся файл привязки topo2.tfw. Содержимое выглядит примерно так:

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

2.6842605828
0.0000000000
0.0000000000
-2.6842605828
2299790.2694883929
486549.6184683086
Безо всяких GDAL'ов напишем простенький скрипт, превращающий это в файл topo2.TAB:

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

!table
!version 300
!charset WindowsCyrillic

Definition Table
  File "topo2.tif"
  Type "RASTER"
  (2299788.92735810,486550.960598600) (0,0) Label "Pt 1",
  (2299791.61161868,486550.960598600) (1,0) Label "Pt 2",
  (2299788.92735810,486548.276338017) (0,1) Label "Pt 3"
  CoordSys Earth Projection 8, 1001, "m", 64.0333333333, 0.0833333333, 1, 2300000, -5700200
  Units "m"
Координаты в скобках вычисляем так:

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

( $5 - $1 * 0.5 , $6 - $4 * 0.5 )
( $5 + $1 * 0.5 , $6 - $4 * 0.5 )
( $5 - $1 * 0.5 , $6 + $4 * 0.5 )
где цифра с долларом означает число из соответствующей строки topo2.tfw.

Можно написать скрипт, чтобы разом обрабатывал все файлы TFW в папке. Теперь точно всё! Или я неправильно понимаю задачу?
The purpose of computing is insight, not numbers

Ответить

Вернуться в «Координаты и привязка»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 9 гостей