Привязка карт с поехали орг

[Goldfinch]

По GDAL не подскажете? По описаниям вроде бы появилась возможность конвертировать карту с привязкой Ozi в geotiff. Мне нужно в geotiff большой объем топокарт перевести. Ставил QGIS 1.7.0, пробовал под ним делать, пишет не понимает формат.

[Moskus]

Ну если в qgis то http://gis-lab.info/blog/2010-03/ozi-qgis/ Если gdal самостоятельно, то при наличии одноименных *.jpg и *.map должна работать строка: <i>gdalwarp -t_srs "EPSG:<b>тут_вписать_код</b>" input.jpg output.tif</i> где код - EPSG-номер системы координат, в которой находятся координаты в .map-файле. Если СК там разные, то предварительно раскидать файлы по спискам (которые пригодны для обхода файлов batch-скриптом) через какой-нибудь искатель файлов по содержимому.

[AlexWhiter]

Moskus, а насколько корректным будет использование в t_srs чего-либо глобального, типа WGS84 и Latlon - чтобы не цепляться за систему координат входных файлов?
Если я правильно понимаю, то в t_srs задается та система, которая будет именно на выходе, а входную программа определяет сама. Ведь так?

[Goldfinch]

<b>Moskus</b>, спасибо за наводку. Получил tiff, но почему-то без привязки. В GM можно это сделать без проблем, но хочется освоить свободный софт)

[Moskus]

Тут вот какое дело (распишу чуть поподробнее, чтоб было).
В промышленных форматах данные о привязке растра состоит из двух частей:

1. Опорные точки.
Они могут храниться в разных формах:
а) как соответствие пар координат в пикселях и пар координат в используемой СК (что допустимо для нетрансформированных растров),
б) например, как одно такое соответствие и параметры разрешения в единицах СК на пиксель (это допустимо только для растров, чья форма приведена в соответствие с используемой СК).
В свою очередь, эти данные могут быть внедрены в сам растровый файл, либо храниться отдельно (в world-файле, в .map-файле и т.п.).

2. Система координат.
Практически та же ситуация с СК. Она может быть указана (через ее название или код) либо описана, как это делается в WKT-формате (то, что содержится внутри .prj-файлов).
Точно также, эти данные могут быть внедрены в сам растровый файл или храниться отдельно.

Пара вида "растр+.map" с этой точки зрения представляет собой нетрансформированный растр и внешний файл, который содержит в себе набор из пар координат точек и название системы координат.

Следовательно, чтобы получить из "растр+.map", скажем, geotiff, нам, в принципе, достаточно только указать вручную, в какой СК указаны координаты. Но поскольку в geotiff должен храниться трансформированный растр, необходимо его еще и трансформировать.

Когда выполняется приведенный мной выше пример команды, происходит следующее: из .map берутся точки привязки, их исходная СК (та что s_srs) при этом принимается равной указанной целевой (t_srs) и растр трансформируется в соответствии с этим, после чего в выходной трансформированный geotiff попадает описание привязки уже в искомом формате и описание СК. То есть, указывая только t_src, мы не осуществляем перепроецирование в новую СК, а осуществляем только "выправление" растра, приводя его форму в соответствие с исходной СК, в которой мы имеем точки в .map.

Если есть желание сразу перепроецировать данные во что-то еще, то нужно указывать в s_srs ту СК, которая указана в .map, а в t_srs - ту, в которую хотим перепроецировать.
Однако у этого есть некоторый минус - если это сканированная карта вместе с рамкой, рамку ведь нужно еще и обрезать. В случае советской номенклатуры границы проходят по "круглым" значениям в градусах. А если перегнать растр из Пулково-1942 в WGS84, они такими круглыми уже не будут. Хотя все это, конечно, зависит от способа обрезки.

UPD: Особо адский случай, конечно, будет представлять конвертирование привязки тех растров, для которых "умельцы" правили datums.dat в OZI, т.е. название СК в .map есть, но чему оно реально соответствует - никому не известно. Такие файлы в сети встречаются. Как и те, которые бездумно привязаны с ошибочным датумом, что вызывает традиционный сдвиг.
Так что я бы, конечно, на полностью автоматическую обработку не стал особо надеяться...
Ну и, несмотря на свое положительное отношение к GDAL, для возни с .map предпочитаю GlobalMapper. Тем более что сбор мозаики там прост и приятен благодаря возможности выбирать все листы в строке или столбце и прописывать clipping им всем сразу по двум границам. Ведь, скажем, те же 200 листов - это всего лишь 4 строки + 5 столбцов.

[Moskus]

А GDAL (FWTools) свежий установлен?
Я попробую на выходных проверить, где может быть затык, потому что описанная процедура - не из документации (там этого нет просто).

[AlexWhiter]

Огромное спасибо за развернутый ответ!

[Moskus]

Не за что, тем более что вот народ еще жалуется, что не работает. :)
Буду проверять и смотреть, почему.

[Moskus]

Как говорится, докладаю. Взял с poehali.org лист 200k--n37-02, сам лист там в GIF, привязка выглядит вот так (кое-что пропущено, а файл я для собственного удобства переименовал): <pre>OziExplorer Map Data File Version 2.2 1.gif 1.gif 1 ,Map Code, Pulkovo 1942 (2),WGS 84, 0.0000, 0.0000,WGS 84 Reserved 1 Reserved 2 Magnetic Variation,,,E Map Projection,Transverse Mercator,PolyCal,No,AutoCalOnly,No,BSBUseWPX,No ...тут было 9 точек привязки... Projection Setup, 0.000000000, 39.000000000, 1.000000000, 500000.00, 0.00,,,,, Map Feature = MF ; Map Comment = MC These follow if they exist Track File = TF These follow if they exist Moving Map Parameters = MM? These follow if they exist MM0,Yes MMPNUM,4 MMPXY,1,102,101 MMPXY,2,2530,91 MMPXY,3,2559,2994 MMPXY,4,89,3003 MMPLL,1, 37.000000, 56.000000 MMPLL,2, 38.000000, 56.000000 MMPLL,3, 38.000000, 55.333333 MMPLL,4, 37.000000, 55.333333 MM1B,24.390244 MOP,Map Open Position,0,0 IWH,Map Image Width/Height,2612,3034 </pre><p>То есть это следующая система координат: Pulkovo-1942/Transverse Mercator, central meridian 39, false easting 500000, scale factor: 1 Что равно проекции Гаусса-Крюгера седьмой шестиградусной зоны, но без указания зоны в false easting (в отечественных картах оно таки используется, соответственно, false easting должно быть 7500000). Далее, я выяснил, как по-человечески обзавестись свежей версией GDAL. Под Windows бинарники собираются в составе двух инсталляторов: FWTools и OSGeo4W. Выбрал OSGeo4W, хотя в FWTools вроде бы тоже последняя стабильная версия входит. Проверяем версию GDAL: <pre>>gdalinfo --version GDAL 1.8.0, released 2011/01/12</pre><p>Интересуемся, что знает GDAL про растр+.map: >gdalinfo 1.gif Driver: GIF/Graphics Interchange Format (.gif) Files: 1.gif Size is 2612, 3034 Coordinate System is `' GCP Projection = PROJCS["unnamed", GEOGCS["Pulkovo 1942", DATUM["Pulkovo_1942", SPHEROID["Krassowsky 1940",6378245,298.3, AUTHORITY["EPSG","7024"]], TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12], AUTHORITY["EPSG","6284"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4284"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",39], PARAMETER["scale_factor",1], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]] GCP[ 0]: Id=, Info= (102,101) -> (375221.866305868,6210477.54231118,0) .... GCP[ 8]: Id=, Info= (2559,2994) -> (436537.571164519,6134902.2022114,0) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 3034.0) Upper Right ( 2612.0, 0.0) Lower Right ( 2612.0, 3034.0) Center ( 1306.0, 1517.0) Band 1 Block=2612x1 Type=Byte, ColorInterp=Palette Metadata: GIF_BACKGROUND=0 Color Table (RGB with 256 entries) 0: 99,145,164,255 ... ═ 255: 255,200,225,255 По-русски: GIF с 256-цветной палитрой, девять точек привязки, описание датума и проекции вполне адекватные (для меня небольшая новость, до того считал, что только точки читаются), но не совсем то, что мы хотим. ОК, значит нам нужно сначала развернуть 8-битный файлик до 24 бит, учитывая, что он будет вращаться и деформироваться, а потом трансформировать в проекцию, которая будет уместной и правильной. Пишем вот такую команду: <pre>>gdal_translate -expand rgba -a_nodata none -of GTiff 1.gif rgba.tif</pre><p>Этой командой конвертируем GIF в GeoTIFF, внедряя в его метаданные информацию о точках привязки и о проекции, а за одно делая его 32-битным с каналом прозрачности. Метаданные получаются такого вида: Driver: GTiff/GeoTIFF Files: rgba.tif Size is 2612, 3034 Coordinate System is `' GCP Projection = 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",39], PARAMETER["scale_factor",1], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]]] GCP[ 0]: Id=1, Info= (102,101) -> (375221.866305868,6210477.54231118,0) ... GCP[ 8]: Id=9, Info= (2559,2994) -> (436537.571164519,6134902.2022114,0) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=PIXEL Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 3034.0) Upper Right ( 2612.0, 0.0) Lower Right ( 2612.0, 3034.0) Center ( 1306.0, 1517.0) Band 1 Block=2612x1 Type=Byte, ColorInterp=Red Mask Flags: PER_DATASET ALPHA Band 2 Block=2612x1 Type=Byte, ColorInterp=Green Mask Flags: PER_DATASET ALPHA Band 3 Block=2612x1 Type=Byte, ColorInterp=Blue Mask Flags: PER_DATASET ALPHA Band 4 Block=2612x1 Type=Byte, ColorInterp=Alpha Если открыть этот файл в GlobalMapper, получится хрен знает что (он будет перевернут), но он еще не трансформирован в правильную форму. Сделать это легко: <pre>>gdalwarp -t_srs EPSG:28407 rgba.tif -of warped.tif</pre><p>Код EPSG найден на <A HREF="http://spatialreference.org">http://spatialreference.org</A> по запросу "Pulkovo 1942 Gauss-Kruger" (тут - внимательно, потому что там всяких вариантов много, нужны те, которые вот отсюда начинаются: <A HREF="http://spatialreference.org/ref/epsg/28 ... . В итоге имеем корректный файл, в котором незначащие пиксели - прозрачные, указана верная проекция и т.п. Для пущей красивости можно пошаманить с параметрами gdalwarp, такими, как - r (near|bilinear|cubic|cubicspline|lanczos) - режим ресемплинга -multi - параллельная обработка на многоядерном процессоре Результирующий warped.tif нормально открывается в GlobalMapper и вообще сам по себе стандартен.

[Moskus]

Использовать географическую проекцию смысла нет из-за проблем с масштабом. А WGS84/Mercator - можно. EPSG:3395