Геокодирование спутниковых снимков: Тайлы

Это финальная статья о том, как геокодировать спутниковые снимки. Здесь я постараюсь описать, а что же делать дальше с полученным GeoTIFF файлом. Если интересно о том, как получить GeoTIFF, то можно почитать предыдущие статьи:

  1. Введение
  2. Опорные точки
  3. GeoTIFF
  4. Тайлы

QGIS

GeoTIFF, полученный на предыдущем шаге выглядит правдоподобно, но хочется убедиться, что проекция выполнена правильно. Самый простой способ - это наложить файл на настоящую карту. Для этого достаточно взять бесплатную QGIS и добавить слой OpenStreetMap:

После этого добавить слой из GeoTIFF файла:

Такое отображение не слишком полезно, поэтому нужно добавить прозрачности. В свойствах слоя нужно сделать общую прозрачность процентов на 60% и проставить “No Data Value = 0”.

Выглядит значительно лучше, но видно, что изображение совсем не совпадает с картой. В моём примере разница примерно 34км. Это достаточно много.

Ошибки геокодирования

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

  1. Ошибки модели SGP4
  2. Неправильно выбранная проекция лучей камеры на Землю
  3. Неправильно посчитанные кватернионы
  4. Неправильная ориентация самого спутника на орбите
  5. Неизвестный фактор

Каждая из этих причин может вносить существенную ошибку в геокодировании. И для разных спутников решение может быть совершенно разным. В случае Метеор-М2, я заметил, что в одной из работ указывается стандартное значение крена 2.5 градуса. Это значит, что камера спутника не точно смотрит в надир, а под небольшим углом. Я попытался смоделировать это в Rugged:

LOSBuilder losBuilder = new LOSBuilder(rawDirs);
losBuilder.addTransform(new FixedRotation("fixed-rotation", Vector3D.PLUS_I, FastMath.toRadians(2.5)));
TimeDependentLOS lineOfSight = losBuilder.build();

Пересоздав GeoTIFF, я получил следующее:

Разница составила около 10км. Это значительно лучше, чем было в прошлый раз. Видимо, я на правильном пути. Немножко “поиграв” с углами крена и тангажа, я высчитал их для снимка как 2.6 и 0.6:

LOSBuilder losBuilder = new LOSBuilder(rawDirs);
losBuilder.addTransform(new FixedRotation("fixed-rotation", Vector3D.PLUS_I, FastMath.toRadians(2.6)));
losBuilder.addTransform(new FixedRotation("fixed-rotation", Vector3D.PLUS_J, FastMath.toRadians(0.6)));
TimeDependentLOS lineOfSight = losBuilder.build();

В результате даёт:

Тайлы

Изображение идеально совпадает с картой. Что же ещё можно сделать? Рассказать об этом другим! Самый простой способ - это разрезать изображение на тайлы и наложить на одну из известных карт.

В GDAL есть встроенная программа gdal2tiles. С её помощью можно достаточно просто сгенерировать тайлы:

python3 gdal2tiles.py --profile=mercator -z 3-6 output.tif tiles

Эта команда создаст в папке tiles необходимые тайлы с уровнями детализации 3,4,5,6. Более детальные тайлы не имеет смысла создавать, так как разрешение МСУ-МР километр на пиксель.

После того, как созданы тайлы, их можно загружать на карту (да, снизу карта и её можно двигать и зумить):

Послесловие

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

Но я не собираюсь стоять на месте. Есть множество вещей, которые можно улучшить:

  1. Более точное геокодирование на основе береговых линий. Дело в том, что спутник может незначительно менять направление камеры на Землю. Гарантировать значения 2.6 и 0.6 нельзя, поэтому нужно сделать “доводку” изображения.
  2. Нужно создать сервис, который бы собирал изображения с Метеор-М2 и мог бы накладывать их на карту в реальном времени. Это вполне возможно, и я планирую этим заняться в ближайшее время.