Задание привязки растра через GDAL для использования в Mapinfo
-
- Активный участник
- Сообщения: 103
- Зарегистрирован: 09 окт 2009, 16:49
- Репутация: 10
Задание привязки растра через GDAL для использования в Mapinfo
Пытаюсь разобраться с 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 метров.
Подскажите, в чем ошибка?
Есть векторный файл в формате 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 метров.
Подскажите, в чем ошибка?
-
- Гуру
- Сообщения: 5292
- Зарегистрирован: 09 апр 2010, 19:30
- Репутация: 1015
- Ваше звание: просто мимо прохожу
- Откуда: Ё-бург
Re: Задание привязки растра через GDAL для использования в Mapinfo
попробуй +towgs84=24,-123,-94 -0.02,0.25,0.13,1.1
что за скрипты?
что за скрипты?
- Игорь Белов
- Гуру
- Сообщения: 2230
- Зарегистрирован: 04 янв 2011, 22:00
- Репутация: 1503
- Откуда: Казань
Re: Задание привязки растра через GDAL для использования в Mapinfo
Забудьте про формат PROJ.4, в данном случае он больше не работает.
Про параметр TOWGS84 тем более забудьте, трансформация теперь не является атрибутом проекции.
GDAL вслед за PROJ перешёл на формат описания систем координат WKT.
Для пересчета использем строку вида:
Файл msk45z2.wkt содержит описание системы координат «МСК-45 зона 2»:
P. S. Поправил BBOX для правильного отображения в QGIS.
Про параметр TOWGS84 тем более забудьте, трансформация теперь не является атрибутом проекции.
GDAL вслед за PROJ перешёл на формат описания систем координат WKT.
Для пересчета использем строку вида:
Код: Выделить всё
ogr2ogr -f "MapInfo File" -a_srs msk45z2.wkt out_file.tab in_file.tab
Код: Выделить всё
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]]]
- Вложения
-
- msk45z2.zip
- (604 байт) 297 скачиваний
The purpose of computing is insight, not numbers
-
- Активный участник
- Сообщения: 103
- Зарегистрирован: 09 окт 2009, 16:49
- Репутация: 10
Re: Задание привязки растра через GDAL для использования в Mapinfo
Оооо, спасибо огромное! Все заработало!
Только почему-то мапинфо отказывается показывать установленную систему координат, пишет просто "Широта/Долгота".
В результате пересчета файла получилось:
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 вроде нет такой настройки.
Только почему-то мапинфо отказывается показывать установленную систему координат, пишет просто "Широта/Долгота".
В результате пересчета файла получилось:
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 вроде нет такой настройки.
-
- Активный участник
- Сообщения: 103
- Зарегистрирован: 09 окт 2009, 16:49
- Репутация: 10
- Игорь Белов
- Гуру
- Сообщения: 2230
- Зарегистрирован: 04 янв 2011, 22:00
- Репутация: 1503
- Откуда: Казань
Re: Задание привязки растра через GDAL для использования в Mapinfo
Различаются:
- количество знаков после запятой в числе 64.033333333333;
- пределы Bounds.
The purpose of computing is insight, not numbers
-
- Активный участник
- Сообщения: 103
- Зарегистрирован: 09 окт 2009, 16:49
- Репутация: 10
Re: Задание привязки растра через GDAL для использования в Mapinfo
да, что-то знаки после запятой я плохо посчитал. Теперь мапинфо все распознает. Спасибо!
- dab
- Гуру
- Сообщения: 671
- Зарегистрирован: 16 дек 2011, 20:02
- Репутация: 170
- Ваше звание: Гуру
- Откуда: Москва
- Контактная информация:
Re: Задание привязки растра через GDAL для использования в Mapinfo
Игорь, а уточните, пожалуйста по параметру
BBOX[54.1,62.5,56.8,65.5]]]
При
PARAMETER["Longitude of natural origin",64.0333333333333,
А почему ... 62.5 ... 65.5 ... а не ... 62.5333333333333 ... 65.5333333333333 ... ?
BBOX[54.1,62.5,56.8,65.5]]]
При
PARAMETER["Longitude of natural origin",64.0333333333333,
А почему ... 62.5 ... 65.5 ... а не ... 62.5333333333333 ... 65.5333333333333 ... ?
- Игорь Белов
- Гуру
- Сообщения: 2230
- Зарегистрирован: 04 янв 2011, 22:00
- Репутация: 1503
- Откуда: Казань
Re: Задание привязки растра через GDAL для использования в Mapinfo
dab, BBOX нужен при выборе системы координат для отображения области покрытия на карте. Пределы даются в WGS 84, не в Пулково. Сотая доля градуса соответствует километру на местности, достаточно для описания области действия системы координат малого населённого пункта.
Пример: EPSG:2935 "Pulkovo 1942 / CS63 zone A1"
Пример: 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
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
Re: Задание привязки растра через GDAL для использования в Mapinfo
но в WKT ее можно как-то указать? или можно использовать только прошитые в EPSG глобусы?
- Игорь Белов
- Гуру
- Сообщения: 2230
- Зарегистрирован: 04 янв 2011, 22:00
- Репутация: 1503
- Откуда: Казань
Re: Задание привязки растра через GDAL для использования в Mapinfo
- Если в рамках темы мечтать, что GDAL будет писать в таблицы MapInfo пользовательский датум, начинающийся с девяток, то нет, GDAL этим не занимается.
- Понятно, что явно прописывать параметры трансформации, ассоциированной с датумом в EPSG, нет смысла. Софт позволяет выбирать трансформацию для датума из БД.
- Ответ на вопрос о возможности произвольной трансформации ДА.
Дано. В MapInfo имеются таблицы в проекции "GK zone 12 (Pulkovo 1942)", например.
Найти. Требуется создать СК такую, чтобы:
- объекты из таблиц MapInfo отображались в QGIS там, где их отображает MapInfo;
- при экспорте слоёв в таблицы 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
-
- Гуру
- Сообщения: 5292
- Зарегистрирован: 09 апр 2010, 19:30
- Репутация: 1015
- Ваше звание: просто мимо прохожу
- Откуда: Ё-бург
Re: Задание привязки растра через GDAL для использования в Mapinfo
по мне - так проще написать пару строчек на C# и решить задачу в лоб
-
- Гуру
- Сообщения: 4056
- Зарегистрирован: 15 окт 2010, 08:33
- Репутация: 1054
- Ваше звание: программист
- Откуда: Казань
-
- Гуру
- Сообщения: 5292
- Зарегистрирован: 09 апр 2010, 19:30
- Репутация: 1015
- Ваше звание: просто мимо прохожу
- Откуда: Ё-бург
Re: Задание привязки растра через GDAL для использования в Mapinfo
возможно тут можно что то сделать через MITAB_BOUNDS_FILE
- Игорь Белов
- Гуру
- Сообщения: 2230
- Зарегистрирован: 04 янв 2011, 22:00
- Репутация: 1503
- Откуда: Казань
Re: Задание привязки растра через GDAL для использования в Mapinfo
Что-то не пойму я постановку задачи.
Допустим, не устраивает точность привязки растров из Интернета, и векторные файлы содержат опознаки, к которым растры будут подтягиваться. Ну так геопривязка растров в этом случае — большая ручная работа с каждой картинкой индивидуально, какая уж тут автоматизация.
Значит, надо полагать, привязка растров устраивает? Но для чего тогда нужны какие-либо векторные файлы?
Склеим в SASPlanet в формат GeoTIFF растр topo1.tif. Трансформируем в МСК-45 зона 2:
Код: Выделить всё
gdalwarp -t_srs msk45z2.wkt topo1.tif topo2.tif
Лень каждый раз отвечать на вопрос о системе координат? OK, добавим параметр:
Код: Выделить всё
gdalwarp -t_srs msk45z2.wkt -co TFW=YES topo1.tif topo2.tif
Код: Выделить всё
2.6842605828
0.0000000000
0.0000000000
-2.6842605828
2299790.2694883929
486549.6184683086
Код: Выделить всё
!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 )
Можно написать скрипт, чтобы разом обрабатывал все файлы TFW в папке. Теперь точно всё! Или я неправильно понимаю задачу?
The purpose of computing is insight, not numbers
Кто сейчас на конференции
Сейчас этот форум просматривают: нет зарегистрированных пользователей и 18 гостей