Построение изолиний на карте мира при помощи Python Basemap

Как известно, Python — один из самых популярных и мощных языков программирования с практически неограниченным функционалом. В данном HOWTO мы рассмотрим как с его помощью построить изолинии на карте мира по точечным данным.

Введение

В ходе работы мне потребовалось простое и масштабируемое решение для постройки карт метеополей по точечным данным (данных метеостанций). Существует программа OpenGrADS, которая позволяет работать с множеством форматов данных, однако решение для отображения станционных данных мне показалось не слишком элегантным. Мой выбор пал на python из-за простоты его использования, а также мощной библиотеки Basemap, обладающей необходимыми методами отображения данных.

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

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

Подготовка к работе

Для работы потребуется установленные пакеты numpy, scipy, matplotlib и basemap. Для установки первых трех достаточно последовательно ввести в терминале:

pip install numpy
pip install scipy
pip install matplotlib

Пакет Basemap можно скачать с официальной страницы на SourceForge. Для установки данного пакета требуется последовательно выполнить следующий набор команд:

cd basemap-*/geos-*/
export GEOS_DIR=/usr/local/geos
./configure --prefix=$GEOS_DIR
make
make install
cd basemap-*
python setup.py install

Здесь /usr/local/geos — директория, в которую вы хотите выполнить установку геофизических данных.

Если в ходе установки не возникло ошибок, вы можете проверить работоспособность данной библиотеки, запустив тестовый скрипт из папки examples:

cd examples
python simplest.py

Если все работает корректно, на экране вы должны увидеть следующую картину:

Построение изолиний

Можно заметить, что на изображении выше построена карта высот земной поверхности, представляющая собой изолинии геопотенциала, однако в данном примере набор широт и долгот — одномерный массив, а массив с данными — двумерный, а передо мной стояла задача построить карту по данными в точках (из одномерного массива).

Итак, допустим у нас есть файл с данными в следующем формате:

lat1 lon1 val1
lat2 lon2 val2
...
latN lonN valN

Разбирать механизм считывания и обработки данных я не буду, оставляя это в качестве самостоятельного упражнения для читателя. Скажу лишь, что сделать это можно, например, при помощи numpy.genfromtxt.

Когда данные считаны, пришло время приступить к их обработке. Идея метода заключается в том, что мы сначала создаем регулярную координатную сетку, на которую накладываем исходные значения, а затем интерполируем их в каждый узел созданной стеки.

numcols, numrows = 500, 500
xi = np.linspace(np.min(lons), np.max(lons), numcols)
yi = np.linspace(np.min(lats), np.max(lats), numrows)
xi, yi = np.meshgrid(xi, yi)

Здесь numcols, numrows — это количество узлов сетки по широте и долготе. Чем выше разрешение, тем выше точность отображаемых данных, но тем ниже скорость работы программы.

Если данные были успешно считаны и сохранены в массивы lats, lons и data, то можно создать вспомогательную двумерную сетку значений, на которую будут наложены исходные данные, а значения в каждом узле будут являться результатом интерполяции.

zi = gd((lons, lats), data, (xi, yi), method='linear')

В качестве параметра method можно указать значение nearest, linear или cubic, что соответствует интерполяции методом ближайшего соседа, линейной или кубической. Сравнить результат работы методов можно тут.

Все что остается — это передать полученные массивы методу Basemap.contourf, указав флаг latlon=True, чтобы не производить пересчет широты и долготы в координаты отображаемых пикселей.

m.contourf(xi, yi, data, clevs, cmap=plt.cm.jet, extend='both', latlon=True)

В результате вы должны получить примерно следующую картину:

Заключение

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

[(lat, lon, val), (lat, lon, val), ..., (lat, lon, val)]

второй — максимальное значение, которое будет отображено, все значения выше него будет обрезаны до этого значения. Файл с функцией и примером её вызова доступен на странице проекта на GitHub.