Лабораторная работа №5 - Работа с наложением тепловых карт в Google Maps & Leaflet JS

Цели данной работы:

1. Освоить создание визуальной тепловой карты в виде файла с растровым изображением на основе сеточного представления данных в виде таблицы (GDAL утилиты - входят в комплект установки QGIS - OSGeo4W Shell).
2. Освоить наложение растрового изображения на карты Google или другие карты с использованием библиотеки Leaflet JS.

Постановка задачи

Имеются:

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

• прогноз техногенной аварии или вредных выбросов на опасном объекте хозяйственной деятельности,
• оценка потенциального риска на территории,
• оценка уровня сигнала базовых станций радиосвязи.

Результаты представлены в виде квадратной равномерной сетки - таблицы, содержащей столбцы и строки с данными расчетной величины (X-ый столбец, Y-ая строка).
Направление увеличения номера столбца соответствует направлению на восток, увеличения номера строки - на юг.
Локальная система координат сонаправлена по оси X на восток, по оси Y на север.
Известны также шаг сетки в метрах, координаты юго-западного угла сетки в метрах (в локальной системе координат расчетов), а также широта-долгота цетра локальной системы координат расчетов (в мировой системе координат WGS'84).

Требуется:

Отобразить визуально результаты, наложив цветовой градиентный эквивалент расчетной сетки поверх карты Google Maps или других карт с использованием библиотеки Leaflet JS.
 

Описание инструментария

Google Maps JavaScript API предлагает для решения подобной задачи инструмент, позволяющий наложить на карту-подложку любой местности рисунок в формате PNG (полупрозрачный), с указанием четких границ наложения со всех четырех сторон света - севера, юга, запада и востока. Для отображения карт Google использует проекцию Меркатора, а для внутренней системы координат объектов - WGS'84. Обе эти системы координат - прямоугольные, что позволяет без труда использовать файлы с растровыми изображениями, которые по своему устройству также прямоугольны.

Таким образом, задача сводится к преобразованию сетки (таблицы) с результатми в растровый файл PNG. Но попутно следует выполнить еще несколько действий:

•  перепроецировать локальную систему координат в WGS'84 или в Меркатор;
•  задать прозрачность;
•  задать цветовой градиентный диапазон, отображающий соответствующие значения расчетной величины.

Весь набор инструментов для преобразования растровых данных доступен в утилитах GDAL (OSGeo4W shell).
 

Локальная система координат

За локальную систему координат удобно принять Равновеликую Азимутальную Проекцию Ламберта. Она является параметрической, то есть зависящей в данном случае от некоторой точки на земной поверхности, в близи которой искажения будут минимальны. Ее определение: "+proj=laea +lat_0=ШИРОТА +lon_0=ДОЛГОТА +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs".

Что хотим получить в итоге

Поехали

Загрузите и распакуйте архив с рабочими файлами.

Будем визуализировать гипотетическое распространение выбросов с ТЭЦ, находящейся на набережной Ижевского Пруда. Координаты где находится ТЭЦ - lat0,lon0 - необходимо будет получить самостоятельно.

Преобразуем in.asc >>> in.tif

Сначала превратим наши табличные данные в файл формата ESRI ASCII GRID. Это нам даст возможность манипулировать исходными данными, используя утилиты GDAL.
Для этого нам необходимо добавить несколько строк перед самими данными:

ncols        372
nrows        376
xllcorner    -650.0
yllcorner    -700.0
cellsize     10.0
...данные...

Здесь - число колонок, строк, координаты юго-западного угла в локальной системе координат (метры), шаг сетки (метры).

Получившийся файл - in.asc (есть в архиве) - мы будем далее использовать. Взгляните на него в текстовом редакторе.

Откройте командную строку с утилитами GDAL (OSGeo4W Shell), перейдите командами cd имя_покаталога в ваш каталог с распаковаными рабочими файлами.

Для преобразования наших данных в формат GeoTIFF, а также для перепроецирования из системы координат Ламберта в WGS'84, ну и также для задания прозрачности
наберите следующую команду gdalwarp вместе с указанными опциями:

gdalwarp -s_srs "+proj=laea +lat_0=ШИРОТА-0 +lon_0=ДОЛГОТА-0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs "EPSG:4326" -of gtiff -srcnodata 0 -dstalpha in.asc in.tif

(Замените Широту-0 и Долготу-0 на точку, где находится ТЭЦ)

Значения опций:
-s_srs - исходная система координат (проекция Ламберта)
-t_srs - целевая система координат (WGS'84)
-of - формат целевого растрового файла
-srcnodata - значение, обозначающее "нет данных" - будет использоваться для прозрачности
-dstalpha - создать альфа-канал для полупрозрачности
in.asc - исходный файл с данными (расширение *.asc говорит о том, что это ASCII GRID)
in.tif - целевой растровый файл.

После выполнения указанной команды мы получили растровый файл in.tif, который имеет формат GeoTIFF и пересчитан в систему координат WGS'84. Все пространственные привязки хранятся также в этом файле, и вы смело можете его загрузить как слой в любой ГИС-системе.
Чтобы посмотреть географические данные (и, в том числе и границы, которые нам понадобяться для наложения), необходимо выполнить команду:

gdalinfo in.tif

Полученный файл также имеет альфа-канал.
Количество пикселей в полученном файле будет отличаться от исходных данных и вычисляется автоматически. Вы можете указать требуемое число пикселей:

-ts Ширина Высота

Преобразуем in.tif >>> in.png

Полученный файл in.tif по-прежнему содержит исходные значения расчетной величины. Чтобы отобразить эти значения визуально, необходимо создать файл соответствия значения расчетной величины и отображаемого цвета. Давайте создадим текстовый файл color.txt, содержащий примерно такую информацию:

nv  255 255 255 0
30  255 200 200 50
100 255 120 120 150
150 255 80  80  170
200 255 50  50  200
255 255 10  10  230

здесь, первая колонка содержит значения расчетной величины, остальные четыре колонки - RGBA - красный, зеленый, синий и уровень непрозрачности соответственно в диапазоне 0-255.
Диапазон расчетной величины может быть любым (зависит от задачи). Необходимо задать несколько строк, а промежуточные значения будут синтерполированы.
nv - означает "нет данных".

Теперь, имея такой файл соответствия color.txt (можете взять из архива), мы можем преобразовать наш файл in.tif в in.png:

gdaldem color-relief in.tif color.txt in.png -of png -alpha

Полученный на данном этапе файл in.png уже можно посмотреть в любом графическом редакторе, так как он содержит цветовую информацию.

Наложение на карту Google Maps (или другие карты с использованием библиотеки Leaflet JS)

Осталось наложить файл in.png на карту-подложку Google Maps, указав географические границы наложения (north, south, west, east или для Leaflet JS [[south,west], [north, east]]) , которые мы можем узнать из команды gdalinfo.
Вы можете взять из архива с рабочими файлами заготовку izhevsk-google.html или izhevsk-leaflet.html. Не забудьте указать правильные границы наложения, имя файла, содержащего растровое наложение, и свой ключ от Google Maps API.

Закрепим

Давайте для закрепления материала еще раз повторим эти процедуры, но уже с другими данными и для другой территории.

Давайте наложим растровое изображение 2-км сетки на территорию вокруг столицы Исландии город Рейкъявика и посмотрим на отклонения сетки от линий меридианов и параллелей. Исландия взята как пример достаточно удаленной от экватора территории.

Исходный файл - grid.asc
Заготовка html с объектом-прямоугольником вдоль меридианов/параллелей - grid-google.html (grid-leaflet.html)
Файл соответствия цветов - color-grid.txt
Центр Рейкъявика - lat_0=64.125008, lon_0=-21.848468

Задайте дополнительно команде gdalwarp опцию -ts 2000 1000, чтобы сетка не искажалась.

Не забудьте при наложении правильно указать границы, имя PNG-файла и свой ключ от Google Maps API.

Удачи!

Полезные ссылки:

Использование GDAL для привязки растровых материалов

Формат ASCIIGRID с примерами

Примеры использования инструментов GDAL