RU:OrbView-3

From OpenStreetMap Wiki
Jump to navigation Jump to search

12 января 2012 года USGS открыло в Public Domain снимки со спутника OrbView-3, принадлежащего GeoEye.

Параметры снимков

  • Спутник: OrbView-3
  • Разрешение: 1 метр для ч/б, 4 метра для цветных
  • Покрытие: вся планета, кусочно
  • Период съёмки: 2003—2007 года
  • Облачность: 0—100%
  • Лицензия: Public Domain

Как загрузить

Заходите на сайт Earth Explorer. Увеличиваете на карте нужный регион (или вводите название города в поле поиска), и затем с нажатым Shift рисуете рамку или контур по точкам. Её можно поправить, тягая углы. Затем жмёте «Data Sets», и там внизу ставите галочку напротив «OrbView-3». В Additinal Criteria можно отфильтровать снимки по облачности: у многих она близка к 100%. Уровень обработки 1Gst значительно лучше 1A в холмистой местности, если такие снимки будут. Наконец, жмите Results и получайте список снимков, либо ничего.

Рядом с каждым снимком будет ряд значков. Нога рисует контуры заснятой области на карте, а листок справа от ноги накладывает сам снимок. Четвёртый значок с зелёной стрелкой даёт скачать zip-архив, но только если зарегистрироваться (бесплатно) на этом сайте. Значок с коробкой позволяет составить заказ из нескольких снимков, затем скачать их разом через их программу.

Как обработать

Подключение снимка делится на три этапа: орторектификация (изменение снимка с учётом рельефа), уменьшение глубины цвета (с 16 до 8 бит), загрузка в редактор. Для первых двух шагов потребуется GDAL, для последнего — либо настроенный python, либо mapserver. Подборка ссылок есть на сайте gdal.

Инструкция основана на странице вики гис-лаба и инструкции Miroff.

Скачивание рельефа

У нас есть два открытых источника данных о рельефе: SRTM и ASTER GDEM.

1. SRTM

  • на странице http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp смотрим на картинку и выписываем куда-нибудь номера рядов и столбцов с интересующей нас местностью;
  • скачиваем с http://gis-lab.info/data/srtm-tif выбранные zip-файлы, складываем в один каталог;
  • распаковываем
  • Если файлов несколько - объединяем все файлы SRTM в формате GeoTIFF в единый виртуальный растр:
 gdalbuildvrt srtm.vrt srtm*tif

2. ASTER GDEM

  • Идём на сайт Aster GDEM (http://gdem.ersdac.jspacesystems.or.jp/index.jsp), регистрируемся, заходим в раздел Search.
  • Выбираем "Select tiles by shapefile" и загружаем файл .shp из архива со снимком
  • Если файлов несколько - их надо склеить в один
 gdal_merge.py -o DEM_merged.tif  ASTGTM2_N5[45]E02[67]/*_dem.tif

Вычисление высоты геоида

Этот шаг можно пропустить, но он повысит точность орторектификации. По крайней мере теоретически.

Если в области снимка найдется подходящий трек с тремя координатами, надо забить данные одной его точки по ссылке http://www.unavco.org/community_science/science-support/geoid/geoid.html и получить высоту геоида в точке. Вряд ли она сильно будет варьироваться на протяжении одного снимка, так что одной точки нам хватит. Укажите высоту в параметре -to 'RPC_HEGHT=<height>' на следующем этапе или исключите этот параметр, если вы не смогли/не захотели вычислять высоту.

Орторектификация

Кроме SRTM, очевидно, потребуются два файла: снимок в tiff и инструкции по привязке в _rpc.txt.

 gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -co ZLEVEL=9 -multi -dstnodata 0 -srcnodata 0 -overwrite -ot Int16 \
  -t_srs epsg:3857 -rpc -to 'RPC_DEM=srtm.vrt' -to 'RPC_HEGHT=<height>' -r bilinear <imagery>.tif <imagery>.rec.tif

Динамический диапазон

Орторектифицированные данные получаются в 16-битном формате, который плохо конвертируется в тайлы. Поэтому потребуется преобразовать их к стандартному 8-битному уровню

gdal_translate -co TILED=YES -co COMPRESS=DEFLATE -co ZLEVEL=9 -scale 0 1024 -ot Byte <imagery>.rec.tif <imagery>.trans.tif

0 и 1024 являются практически стандартным преобразованием. Если получившаяся картинка по яркости вам не понравится - попробуйте другие уровни (Miroff дает числа 30 и 400)

Если установлен QGIS (кстати он понимает и показывает данные в 16-битном формате), то яркость можно подобрать с помощью инструмента "растяжение гистограммы" на панели, иконка появится если отметить "Вид-Панели инструментов-Растр" и установлен модуль GDAL Tools (Raster Tools). Попробуйте зазумится на самые яркие или темные места, или на границу снимков - если пытаетесь выровнять яркость, и нажмите на иконку "Локальное растяжение гистограммы". После подбора яркости сохраните проект .qgs и откройте этот файл в текстовом редакторе - найдите значения для параметра scale в секции <contrastEnhancementMinMaxValues><minMaxEntry>. После преобразования, если снимков несколько - можно создать виртуальный растр, предварительно расположив снимки в нужном порядке, следя как они перекрываются. После чего нарезать на тайлы на следующем шаге, указав в качестве обрабатываемого снимка - виртуальный растр *.vrt

Тайлы

Самый простой способ загрузить снимок в JOSM — нарезать его на тайлы и положить в каталог любого веб-сервера, например, Apache. Команда такая (нужен python!):

gdal2tiles.py -w none -n -z 12-18 <imagery>.trans.tif tiles

Можно попробовать улучшить качество добавив ключ -r lanczos или -r bilinear.

После её выполнения каталог tiles будет содержать все тайлы для снимка. Разумеется, тайлы от разных снимков можно объединять. После всего подключить источник в JOSM просто — нужно добавить TMS с адресом:

tms[18]:http://localhost/{zoom}/{x}/{-y}.png

или

tms[18]:file:/C:/path/to/{zoom}/{x}/{-y}.png

Единственный минус — что gdal2tiles делает набор png, тогда как для снимка лучше jpg (иначе получается почти гигабайт на снимок).

Впрочем, png можно немножко "похудеть" специальными программами. Например TruePNG http://x128.ho.ua/pngutils.html

cmd /k "for /r %i IN (*.png) do truepng /o4 /g1 %i"

UPD: Есть патч (спасибо PShA) добавляющий поддержку JPEG и возможность регулировать степень сжатия, так как в GDAL по умолчанию для JPEG используется QUALITY=75

WMS

Другой способ — поднять WMS на базе MapServer. Примерный map-файла для одного снимка:

MAP
  NAME "OrbView-3"
#  CONFIG "CPL_DEBUG" "ON"
  STATUS ON
  UNITS METERS
  SHAPEPATH "/home/user/orbview3"

  PROJECTION
    "init=epsg:4326"
  END
  EXTENT -180 -90 180 90
  TRANSPARENT ON

WEB
    TEMPLATE "template.html"
    IMAGEPATH "/home/user/orbview3/"
    IMAGEURL "/orbview3"
    METADATA
        "ows_enable_request"   "*"
        "wms_title"     "OrbView-3 3V050704P0000801711A520014701682M_001635922"
        "wms_onlineresource" "http://my-server/cgi-bin/mapserv"
        "wms_srs" "EPSG:3857 EPSG:4326"
    END
END

LAYER
    NAME "orbview3rec"
    TYPE RASTER
    STATUS ON
    OPACITY ALPHA
    DATA "3v050704p0000801711a520014701682m_001635922.rec.tif"
    PROJECTION
        "init=epsg:3857"
    END
    METADATA
        "wms_title" "OrbView-3 Rectified"
        "wms_srs" "EPSG:3857 EPSG:4326"
    END
END

END

Для объединения множества снимков можно использовать gdal_merge (имеет смысл только для рядом располагающихся снимков), vrt (для не очень большого количества снимков) или tileindex.

Готовый WMS

Выложите скачанный архив со снимком куда-нибудь в доступное место (например, на dropbox). Оповестите Alexandr Zeinalov (shurik в IRC). Подождите.

Информация о уже обработанных снимках есть на сайте, вместе с картой и списком. Для подключения слоя OrbView-3 в JOSM добавьте следующую подложку WMS:

wms:http://osm.sbin.ru/ov3/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=orbview3rec,orbview3pv&SRS={proj(EPSG:4326)}&WIDTH={width}&height={height}&BBOX={bbox}&format=image/jpeg

Исходный код доступен по адресу https://github.com/shurshur/ov3wms. Используется mapserver с tileindex.

Ссылки