MapServer + GeoTIFF, проблемы с охватом и переходом между СК

[jekhor]

День добрый. Пытаюсь загнать в MapServer лист топокарты. Имеется .gif с привязкой для Ozi. Конвертирую его в GeoTIFF с помощью http://ozi2geotiff.svn.sourceforge.net/ ... geotiff.py . Получаю TIFF:

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

$ gdalinfo 050k--n35-079-4.tiff 
Driver: GTiff/GeoTIFF
Files: 050k--n35-079-4.tiff
Size is 2788, 3107
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (3032446.759999999776483,7105097.959999999962747)
Pixel Size = (10.699968699999999,-10.699968699999999)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DOCUMENTNAME=/tmp/mapconverter_4bJ7DP.tif
  TIFFTAG_SOFTWARE=ImageMagick 6.3.7 12/10/08 Q16 http://www.imagemagick.org
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 3032446.760, 7105097.960) ( 27d14'27.36"E, 53d50'18.43"N)
Lower Left  ( 3032446.760, 7071853.157) ( 27d14'27.36"E, 53d39'41.21"N)
Upper Right ( 3062278.273, 7105097.960) ( 27d30'32.09"E, 53d50'18.43"N)
Lower Right ( 3062278.273, 7071853.157) ( 27d30'32.09"E, 53d39'41.21"N)
Center      ( 3047362.516, 7088475.559) ( 27d22'29.72"E, 53d45'0.16"N)

Затем создаю MAP-файл:

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

MAP
        IMAGETYPE PNG24
#       EXTENT 0.00 0 2788.0 3107.0
        EXTENT 27.0 53.0 28.0 54.0
        SIZE 800 600 
        SHAPEPATH "/home/ftp/maps/poehali/1_50000/GeoTIFF"
        IMAGECOLOR 255 255 255

        LAYER
                NAME "mil"
                DATA "050k--n35-079-4-rgb.tiff"
                TYPE RASTER
                STATUS ON
                PROCESSING  "BANDS=1,2,3"
                OFFSITE 0 0 0
        END     
END

Смотрю слой в браузере --- ничего не показывает (белая картинка). Опытным путём было выяснено, что если extent задать в точках исходного изображения (см. закомментированную строчку) --- то показывает. Что я делаю не так?

[Максим Дубинин]

EXTENT явно не верный, причем ни один ни другой
у вас же

PROJECTION["Mercator_1SP"]

и глядите какие точки

pper Left ( 3032446.760, 7105097.960) ( 27d14'27.36"E, 53d50'18.43"N)

а вы пытаетесь ему то пиксельный охват задать, то географическую СК

выясните правильный extent с помощью gdalinfo

[jekhor]

Тьфу, не сообразил, в какой проекции задаётся. Вот что бывает, если голову не включать.

Тогда второй вопрос:

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

jek@jeknote:~$ cs2cs +init=epsg:4284 +to +init=epsg:4326
27d0.0000E 55d10.0000N
27dE    55d10'N 0.000

/usr/share/proj/epsg:
# Pulkovo 1942
<4284> +proj=longlat +ellps=krass +no_defs  <>

# WGS 84
<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs  <>

Почему координаты в разных системах одинаковые? Неужели расхождение столь мало?

Собственно, почему спрашиваю (хотя вышесказанное может и не относиться к этому): конвертирую пачку карт с poehali.org (с привязкой в Ozi) в GeoTIFF, указав проекцию epsg:28405, кормлю это всё MapServer'у, затем запрашиваю через WMS в географической системе WGS84, сравниваю с GPS треком --- расхождение метров этак 150. В Ozi всё красиво, если открыть GeoTIFF в qGIS --- спроектированные координаты совпадают с написанными на карте.

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

 # Pulkovo 1942 / Gauss-Kruger zone 5
<28405> +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=5500000 +y_0=0 +ellps=krass +units=m +no_defs  <>

MAP-файл:

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

MAP
        PROJECTION
                "init=epsg:28405"
        END

        IMAGETYPE PNG24
        EXTENT 5301128.910000 5911358.720279 5697295.208612 6209713.230000
        UNITS meters
        SIZE 800 600
        SHAPEPATH "/home/ftp/maps/poehali/1_50000/GeoTIFF"
        IMAGECOLOR 255 255 255

        WEB
                METADATA
                        "wms_title"     "My WMS server"
                        "wms_onlineresource" "http://localhost/cgi-bin/mapserv?map=/home/jek/work/osm/mapserver/mil500.map&"
                        "wms_srs"       "EPSG:28405 EPSG:4326"
                END
        END

        LAYER
                PROJECTION
                        "init=epsg:28405"
                END

                NAME "milidx"
                STATUS ON
                TILEINDEX "mil_idx.shp"
                TILEITEM "LOCATION"
                TYPE RASTER
                METADATA
                        "wms_title" "milidx"
                        "wms_srs"       "EPSG:28405"
                END
        END
END

[jekhor]

Кажется, нашёл нужный набор глобу...тьфу, проекций.

projection = "+init=epsg:3335"
src_projection = "+init=epsg:4179"

Как-то так. Там много разных Pulkovo 1942, оказывается

[Максим Дубинин]

решились ли все остальные упомянутые проблемы?

[jekhor]

решились ли все остальные упомянутые проблемы?

Да. По поводу проекций и разницы при переводе координат, похоже, надо будет как следует почитать документацию по PROJ.4. Но цели своей я достиг --- получил WMS из генштабовских карт с правильной привязкой

[glebius]

Хмм, параллельно занимаемся одним и тем же. Я даже решил по этому поводу создать логин на форуме Может быть объединим усилия?

У меня есть квадраты N-36, N-37, O-36, O-37. Конвертация из ozi в геотифф уже вторые сутки длится, пока N-36 и частично N-37 готовы. Как получится сервер wms, то буду готов держать его включенным круглые сутки. Правда с местом на нём туго, всего под карты 18 Гб есть.

Сейчас пытаюсь освоить mapserver, проблемы те же, что у вас. Вообще в ГИС я полный новичок, и различные проекции для меня пока тёмный лес.

Но цели своей я достиг --- получил WMS из генштабовских карт с правильной привязкой

А у вас сейчас один tif хостится, или же создали shape (надеюсь правильно назвал) и имеете обширную карту составленную из множества листов? В конечном итоге хотелось бы добиться последнего.

[jekhor]

Как выяснилось, всё очень просто. Только если конвертируете ozi2geotiff.py, то его поправить надо на тему проекций. Вот моя версия: http://jekhor.info/files/maps/ozi2geotiff.py . Я делал так, чтобы конечная проекция совпадала с той, которая используется на картах Генштаба, чтобы координаты были теми же.

Про проекции читать где-то тут на ГИСлабе, я разобрался почти с нуля за пол-дня со всем.

После конвертации в GeoTIFF делаете индексный файл (gdaltindex mil_idx.shp *.tiff) и прописываете его в MAP-файл (я выше приводил пример, в нём нужно 28405 заменить на 3335, а так там всё рабочее). Всё. Я это делал, чтобы использовать в качестве подложки в josm (редактор OpenStreetmap), также проверял в qGIS.

Информация по MapServer взята в основном с этого сайта из переведённого Tutorial.

P.S.: Да, вам ещё понадобится обрезать поля карты, наверное. Я этим пока не занимался, не автоматизировал (благо в файлах привязки углы собственно карты есть).
P.P.S.: Насчёт публичного сервера аккуратно, не знаю, как у вас, но у нас в Беларуси буквально на днях Госкомимущество потребовало прикрыть карту на одном из сайтов, сделанную на базе OpenStreetmap (!). Цитата от хозяина сайта: "сказали что данная картографическая информация принадлежит белкартографии, "если хотите работать и спать спокойно - выкупаете права и ставите карту"".
P.P.P.S.: Да, суммарный охват проще всего узнать, натравив ogrinfo -al на индексный shape-файл.

[glebius]

Спасибо за советы! Сейчас буду разбираться. И я тоже это всё делаю для подложки в josm (тсс!)

Про публичный сервис я и не имел в виду. Просто общий сервер открытый исключительно для узкого круга пользователей. Для экономии места на винтах и времени потраченного на обработку и привязку. Да и cgiшник mapserv я бы не стал открывать всему миру даже если за ним не лежит никакой карты Не нравится мне указание конфигурационного файла клиентом

[jekhor]

Совсем упустил из виду. Это у меня 5-я зона, у вас 6-я и 7-я. Проекции будут следующими:


+proj=tmerc +lat_0=0 +lon_0=33 +k=1 +x_0=6500000 +y_0=0 +ellps=krass +towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 +units=m +no_defs

+proj=tmerc +lat_0=0 +lon_0=39 +k=1 +x_0=7500000 +y_0=0 +ellps=krass +towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 +units=m +no_defs

Возможно, параметры +towgs84 здесь надо поменять на другие (как я понял, смотреть сюда: viewtopic.php?p=2359#p2359), но у меня с этими хорошо попадает всё.

[glebius]

А в ozi .map файле проекция указана в пятой строке?

[jekhor]

В пятой строке, как я понял --- используемая географическая СК. Конечная проекция --- в девятой и там, где Projection Setup. В принципе, конечные координаты в .map никак не фигурируют, поэтому проекция может быть любой при конвертации. Логично использовать ту же, в которой карты, чтобы не было искажений растра, и чтобы было легко проверить по сетке правильность.

Тут важно понимать, что градусы и минуты в .map даются именно в той СК, которая в пятой строчке. Похоже, это как раз epsg:4179 (в Ozi --- Pulkovo 1942 (2)). Возможно, там тоже стоит посмотреть на параметры +towgs84, но при переходе в рамках одного датума это не критично, как я понял. Оригинальный ozi2geotiff.py не смотрит на используемые географическую СК и проекцию.

Советую сначала попробовать помучать 2-3 файла, а уже потом конвертировать все.

[glebius]

Среди всех .map что у меня есть, встречаются вот такие пятые строчки:

Pulkovo 1942 (1),, 0.0000, 0.0000,WGS 84
Pulkovo 1942 (1),WGS 84, 0.0000, 0.0000,WGS 84
Pulkovo 1942 (1),WGS 84, 0.0000, 0.0000,WGS 84
Pulkovo 1942 (2),, 0.0000, 0.0000,WGS 84
Pulkovo 1942 (2),WGS 84, 0.0000, 0.0000,WGS 84

Вообще я не понимаю, вроде как согласно http://gis-lab.info/qa/proj-sk-faq.html ГСК Pulkovo 1942 и WGS 84 взаимоисключающи, что же означают записи выше?

Совсем упустил из виду. Это у меня 5-я зона, у вас 6-я и 7-я. Проекции будут следующими:

А как выяснить у какого листа карт ГШ какая проекция?

[glebius]

Так вроде с Пулково 1942 и зонами Гаусса-Крюгера разобрался, теперь могу понять для какого листа какая проекция.

Но почему было принято решение всё переводить в epsg:4179?

[wildMan]

судя по описанию формата файла:

Datum settings
WGS 84,, 0.0000, 0.0000,WGS 84 - Note - for normal maps only the datum in the first field is used the rest of the parameters are for datum shifts and there is normally no need for these to be modified by the user.

строчка

Pulkovo 1942 (1),WGS 84, 0.0000, 0.0000,WGS 84

означает что карта привязана в Пулково но в программе будут отображаться координаты в WGS.
это несложно проверить. откройте карту в OZI и гляньте на координаты угла или центра.

возможно будет полезно - http://shurila.blogspot.com/2008/12/ozf ... -file.html