Методы: WinPython, python, numpy, scipy, matplotlib, basemap
Задачи:
1) Установить на windows Питон;
2) разобраться в основных модулях Питона, которые необходимы для построения "научного конвейера";
3) "Научный конвейер". На примере реанализа NCEP/NCAR разобрать полный цикл обработки гидрометеорологических и океанологических данных: от чтения netcdf файла до рисования карты, пригодной для презентации/публикации.
Хочу предложить читателям моего блога взглянуть на язык программирования Питон/Python. Недавно я начал активно использовать его для решения своих научных задач (построение карт различных метеорологических параметров по данным различных реанализов в формате netcdf). Python позволяет легко читать netcdf формат, что для пользователей windows очень актуально. Ncdump, конечно, хорош, но требует всё равно постобработки. К тому же большие файлы (а таких в научных задачах большинство) он "жуёт" плохо. Установить же библиотеку netcdf на windows, как я понимаю, задача достойная эпитета "эпическая".
А вот Питон/Python работает с netcdf легко и непринуждённо! К тому же сейчас для Питона написаны и свободно распространяются модули, которые позволяют создавать графику, в том числе и географические карты, очень высокого класса. Таким образом, с помощью Питона можно создать конвейер полного цикла по производству научной продукции (чтение данных, обработка и преобразование данных, анализ данных, графическое представление результатов).
Конечно, у Питона есть и недостатки. Но о них мы поговорим как-нибудь потом. А теперь перейдём к решению поставленных задач.
I. Установка Питона. IDLE, WinPython
Для начала нам нужен интерпретатор Питона. Можно скачать его с официального сайта. У Питона есть две ветви развития: вторая (Python 2.*) и третья (Python 3.*). Если кратко, то они отличаются синтаксисом некоторых основных функций интерпретатора Питона. Я буду писать о Python 3.*, который является Питона будущего поколения. Если интересно чем они отличаются, то вам сюда (англ.).
Скачав "Python 3.3.2 Windows X86-64 MSI Installer" вы получите стандартный интерпретатор Питона. В установочный пакет входит среда разработки IDLE, которая отличается незатейливостью и простой. На данном этапе вы получите "моторчик", основной блок, ядро, с помощью которого вы сможете создать первую тестовую программу "Hello, World!".
Для естественно-научной работы необходимы дополнительные пакеты-модули. Их несколько: numpy, scipy, matplotlib, basemap. Это необходимый минимум. Чтобы поставить их на windows, необходимо найти собранные для windows установочные пакеты. Обратите внимание на то, что одновременно должны совпадать архитектура (32 или 64), ветка Питона (2.* или 3.*) и текущая версия (3.2.* и 2.6.*, например) дополнительного пакета и интерпретатора. Совет: по умолчанию качайте пакеты и интерпретатор для 32-х разрядных систем. Как показывает практика, все вкусные пакеты появляются сначала для x32, а затем, может быть, выходят и для x64. А часто и не выходят. Скачать эти модули можно по этим ссылкам: Scipy, Numpy, Matplotlib, Basemap. После загрузки их нужно установить. Здесь всё просто, ошибиться невозможно. Если всё же вышла ошибка проверьте на совпадение "реквизиты" интерпретатора и устанавливаемого пакета.
Таким образом, необходимо вручную искать, скачивать дополнительные модули и устанавливать их на ядро-интерпретатор Питона. Недавно я обнаружил на сайте SciPy альтернативные сборки интерпретатора Питона. Они также есть и на python.org (Alternative Implementations). Они включают в себя много всяческих примочек (среда разработки, предустановленные модули, дополнительные возможности). Я попробовал несколько и остановился на WinPython. Скачиваем, устанавливаем и запускаем среду разработки Spyder.
В окне Console (справа внизу) появится сообщение-приглашение:
Imported NumPy 1.7.1, SciPy 0.12.0, Matplotlib 1.3.0 + guidata 1.6.1, guiqwt 2.3.1
Type "scientific" for more details.
Пишем scientific:
>>> scientific
This is a standard Python interpreter with preloaded tools for scientific
computing and visualization:
>>> import numpy as np # NumPy (multidimensional arrays, linear algebra, ...)
>>> import scipy as sp # SciPy (signal and image processing library)
>>> import matplotlib as mpl # Matplotlib (2D/3D plotting library)
>>> import matplotlib.pyplot as plt # Matplotlib's pyplot: MATLAB-like syntax
>>> from pylab import * # Matplotlib's pylab interface
>>> ion() # Turned on Matplotlib's interactive mode
>>> import guidata # GUI generation for easy dataset editing and display
>>> import guiqwt # Efficient 2D data-plotting features
>>> import guiqwt.pyplot as plt_ # guiqwt's pyplot: MATLAB-like syntax
>>> plt_.ion() # Turned on guiqwt's interactive mode
Within Spyder, this interpreter also provides:
* special commands (e.g. %ls, %pwd, %clear)
* system commands, i.e. all commands starting with '!' are subprocessed
(e.g. !dir on Windows or !ls on Linux, and so on)
То есть мало того, что нам не нужно искать и скачивать три пакета из четырёх, так они ещё и автоматически импортированы для удобства написания кода! Здорово! Однако последний пакет, который отвечает за рисование географических данных, не входит в эту сборку. Значит нужно скачать его и установить самостоятельно. Идём в Basemap и скачиваем установочный пакет для windows. Теперь его нужно установить в WinPython.
Установка пакетов для WinPython
Если просто запустить исполняемый файл, то возникнет ошибка - он будет искать стандартный интерпретатор Питона, установку которого мы рассмотрели выше. Для корректной установки модуля Basemap в WinPython, нужно запустить утилиту wpcp, которая находится в папке *\WinPython-64bit-3.3.2.2\scripts. Запустится оконное приложение, где нажав кнопку "add package" нужно выбрать скачанный установочный пакет. Завершите установку, нажав "Install package". Так можно устанавливать собранные для windows пакеты-модули в WinPython.
II. Обзор основных пакетов-модулей python для научной работы
Установив WinPython-64bit-3.3.2.2, например, и модуль Basemap, мы получили следующее:
1) интерпретатор питона - python.exe (*\WinPython-64bit-3.3.2.2\python-3.3.2.amd64\)
3) редактор Spyder
Интерпретатор даёт нам базовые возможности, с их помощью мы программируем.
Пакеты-модули дают доступ к специализированным функциям. Разберём их подробнее.
1. Numpy. Хорошее введение с примерами синтаксиса дано здесь. А если кратко, то numpy вводит в python класс массивов, позволяет осуществлять матричные операции, преобразование Фурье, работу с случайными числами. Чаще всего от numpy требуется завести массив. Например так:
>>> import numpy as np
>>> a = np.arange(-90,90,30)
>>> a
array([-90, -60, -30, 0, 30, 60])
>>> type(a)
<class 'numpy.ndarray'>
>>> a.ndim
1
>>> len(a)
6
>>> a.shape
(6,)
В результате будет создан объект "a" типа "numpy.ndarray". Это массив из 6 элементов, от [-90,90) c шагом 30. Его форма - одномерный массив, число "осей", соответственно, 1. Чтобы изменить форму массива нужно применить метод reshape.
>>> a.reshape(2,3)
array([[-90, -60, -30],
[ 0, 30, 60]])
Хорошую шпаргалку по осям numpy можно найти здесь.
Для новичков:
Функция len(имяобъекта) независима от имяобъекта и стоит выше него. Она заключает его в скобки. А вот метод имяобъекта.shape "проявляется" только после указания объекта (в нашем случае это массив "a"). Shape - это метод класса numpy.ndarray, то есть способность массива. Поэтому метод указывается после точки. Пример: машина - это класс. Цвет - метод машины. Оттенок - метод цвета. Чтобы получить название оттенка цвета машины, которую мы сейчас представили, согласно логике ООП и питона написать следующее:
Функция len(имяобъекта) независима от имяобъекта и стоит выше него. Она заключает его в скобки. А вот метод имяобъекта.shape "проявляется" только после указания объекта (в нашем случае это массив "a"). Shape - это метод класса numpy.ndarray, то есть способность массива. Поэтому метод указывается после точки. Пример: машина - это класс. Цвет - метод машины. Оттенок - метод цвета. Чтобы получить название оттенка цвета машины, которую мы сейчас представили, согласно логике ООП и питона написать следующее:
то_что_я_представляю = машина.цвет.оттенок
Если нужно что-то поменять в массиве по условию, то это легко сделать с помощью where (numpy.where):
>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> ii = np.where((a >= 2) & (a < 8))
>>> ii
(array([2, 3, 4, 5, 6, 7], dtype=int64),)
>>> a[ii] = a[ii] - 273
>>> a
array([ 0, 1, -271, -270, -269, -268, -267, -266, 8, 9])
ii - массив индексов, которые удовлетворили логическому условию в where.
2. Scipy. Это библиотека (модуль с большим числом функций) математических алгоритмов и удобных функций, базирующихся на модуле Numpy. Я уже кратко показывал возможности scipy для определения квантилей. Именно через модуль scipy можно работать с файлами формата netcdf. Покажу на простом примере. Скачаем данные реанализа NCEP/NCAR. Пусть это будет файл с данные многолетнего среднемесячного приземного давления (pres.mon.ltm.nc). Загружаем файл в папку, где мы сохраним файл питоновского кода.
Для новичков:
# - однострочный комментарий
'''
Всё что между тройными кавычками - многострочный комментарий
'''
Итак, код с комментариями.
# ------------------ Scipy 1 ---------------------------------------
from scipy.io import netcdf as nc # подключаем выбранный модуль netcdf из модуля scipy.io
fname = 'pres.mon.ltm' # имя файла
f = nc.netcdf_file(fname+'.nc') # открываем файл и получим f - идентификатор открытого файла
fkeys = list(f.variables.keys()) # Преобразуем в список кортеж имён переменных, содержащихся в файле
print(fkeys) # печатаем список переменных (для ветки python 2. скобки не нужны)
['time', 'valid_yr_count', 'lat', 'pres', 'lon', 'climatology_bounds']
'''
Определяем массивы времени, долгот, широт, давления.
Причём широты и данные по широтам считываем задом наперёд, то есть в обратном направлении: от -90 до 90 (изначально в реанализе они уложены от 90 до -90).
Это нужно для беспроблемного рисования карт с помощью модуля basemap.
'''
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][:]
v = f.variables['pres'][:]
v1 = f.variables['pres']
'''
Массив v и объект v разные. V - это данные, а v1 позволяет узнать коэффициент масштаба и добавочный коэффициент. Преобразовать массив v в нормальные значения можно так
"''
v = v*v1.scale_factor + v1.add_offset
f.close() # Закрываем файл
# -----------------------------------------------------------------
Итоговый код:
# ------------------ Scipy 2 ---------------------------------------
from scipy.io import netcdf as nc
fname = 'pres.mon.ltm'
f = nc.netcdf_file(fname+'.nc')
fkeys = list(f.variables.keys())
print(fkeys)
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][:]
z = f.variables['climatology_bounds'][:]
v = f.variables['pres'][:]
v1 = f.variables['pres']
v = v*v1.scale_factor + v1.add_offset
f.close()
# -----------------------------------------------------------------
fname = 'pres.mon.ltm'
f = nc.netcdf_file(fname+'.nc')
fkeys = list(f.variables.keys())
print(fkeys)
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][:]
z = f.variables['climatology_bounds'][:]
v = f.variables['pres'][:]
v1 = f.variables['pres']
v = v*v1.scale_factor + v1.add_offset
f.close()
# -----------------------------------------------------------------
3. Matplotlib. Это графический модуль, который позволяет строить как простые двумерных графики, так и трёхмерные изображение и даже географические карты.
Введение в matplotlib на русском разобрано здесь и у Колдунова. Для примера нарисуем график "роза ветров" (wind rose). Будем рисовать его в полярной проекции.
#----------- Matplotlib 1 -------------------------------------
from math import pi
import numpy as np
import matplotlib.pyplot as plt
angles = np.arange(0,2*pi,pi/4.0) # Задаём массив направлений
r = np.random.random(len(angles)) # ... и случайный массив значений
fig = plt.figure() # Создаём новое полотно, на котором будем рисовать
ax = fig.add_subplot(111, projection='polar') # ax - полотно в полярной проекции, занимающее
# всё пространство полотно fig
# Рисуем график и его название
ax.plot(angles, r, color='r', linewidth=1.)
ax.set_title("Wind rose", va='bottom')
plt.show()
Введение в matplotlib на русском разобрано здесь и у Колдунова. Для примера нарисуем график "роза ветров" (wind rose). Будем рисовать его в полярной проекции.
#----------- Matplotlib 1 -------------------------------------
from math import pi
import numpy as np
import matplotlib.pyplot as plt
angles = np.arange(0,2*pi,pi/4.0) # Задаём массив направлений
r = np.random.random(len(angles)) # ... и случайный массив значений
ax = fig.add_subplot(111, projection='polar') # ax - полотно в полярной проекции, занимающее
# всё пространство полотно fig
# Рисуем график и его название
ax.plot(angles, r, color='r', linewidth=1.)
ax.set_title("Wind rose", va='bottom')
plt.show()
Получаем следующее изображение:
В глаза сразу бросается 90 градусов вверху графика. Также наша роза не является замкнутой. Необходимо привести это изображение в порядок. Для этого добавим следующие строчки в код:
ax.set_theta_direction(-1) # теперь увеличение угла theta будет по часовой (-1 - увеличение угла
# идёт против часовой стрелки)
ax.set_theta_offset(pi/2.0) # изменим положение нуля (ставим его на положение 'N', которое
# сейчас равно 90 градусам)
ax.plot((angles[-1],angles[0]),(r[-1],r[0]), color='r', linewidth=1.) # Добавляем ещё один plot-линию,
# которая соединяет последнюю и первую точки.
Сохраняем получившийся рисунок в файл с расширением 'png' (можно и другой) с помощью
plt.savefig('rose.png')
Вот так намного лучше! А теперь выведем сразу две розы:
#----------- Matplotlib 2-------------------------------------
from math import pi
import numpy as np
import matplotlib.pyplot as plt
angles = np.arange(0, 2*pi, pi/4.)
fig = plt.figure(figsize=(10,5)) # зададим явно размер полотна, чтобы графики не перекрывались
for i in range(1,3,1):
r = np.random.random(len(angles)) # для каждой розы значения будут разные
ax = fig.add_subplot(120+i, projection='polar')
ax.plot(angles, r, color='r', linewidth=1.)
ax.plot((angles[-1],angles[0]),(r[-1],r[0]), color='r', linewidth=1.)
ax.set_theta_direction(-1)
ax.set_theta_offset(pi/2.0)
ax.set_title("Wind roses "+str(i), loc = 'center')
plt.show()
plt.savefig('rose_subplot.png')
#--------------------------------------------------------------
Здесь мы использовали subplot - разбили наше полотно на таблицу. Вписывать график мы будем уже в ячейку таблицы. Схема таблицы следующая:
SUBPLOT
ax = fig.add_subplot(121, projection='polar')
121 - rcn: r = 1, c = 2, n = 1
r - число строк в таблице
c - число столбцов в таблице
n - номер ячейки в получившейся таблице
Тогда, чтобы нарисовать графики друг под другом нужно будет задать значения для subplot 211 и 212. В приведенном выше коде этот номер изменяется с помощью переменной цикла i.
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
fname = 'pres.mon.ltm.nc'
f = nc.netcdf_file(fname)
fkeys = list(f.variables.keys())
print(fkeys)
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][:]
v = f.variables['pres'][:]
v1 = f.variables['pres']
f.close()
v = v*v1.scale_factor + v1.add_offset
m = Basemap(projection='mill')
m.drawcoastlines(linewidth=1)
m.drawcountries(linewidth=0.5)
m.drawmapboundary()
parallels = np.arange(-90,90+1,30)
meridians = np.arange(0,360+1,60)
m.drawparallels(parallels,labels=[True,True,True,True])
m.drawmeridians(meridians,labels=[True,False,False,True])
# Данные по широтам и долготам, взятые из netcdf файла необходимо перевести из векторов
# в матрицы. Делается это с помощью функции meshgrid
x1,y1 = np.meshgrid(x,y)
lons,lats = m(x1,y1) # задаем на картографической основе преобразованные в матрицы
# координаты
'''
Другой синтаксис перевода координат x и y в координаты, понятные карт. основе
m.contourf(lons,lats,v[1,:,:],40) # Рисуем на картографической основе m данные за февраль
# contourf рисует изолинии с послойной окраской через 40 единиц по сеточным данным
# Массив v[1,:,:] - срез от всего массива v, представляющий данные за один месяц февраль
plt.title('Test map N') # Название карты
plt.savefig('Test_basemap_map.png')
plt.show()
# ------------------------------------------------------------------------------
Запускаем скрипт и получаем ... вот это!
Что-то явно было сделано с ошибками. По порядку:
1) Нарисовано только половина карты!
2) Непонятные белые полосы сверху и снизу карты!
Начнём с главного. Код рисует половину карты потому что данная проекция рассчитана на изменение долгот в пределах от -180 до 180, а в наших данных реанализа долготы изменяются в других пределах, а именно от 0 до 357.5:
>>> y,x
[ 90. 87.5 85. 82.5 80. 77.5 75. 72.5 70. 67.5 65. 62.5
60. 57.5 55. 52.5 50. 47.5 45. 42.5 40. 37.5 35. 32.5
30. 27.5 25. 22.5 20. 17.5 15. 12.5 10. 7.5 5. 2.5
0. -2.5 -5. -7.5 -10. -12.5 -15. -17.5 -20. -22.5 -25. -27.5
-30. -32.5 -35. -37.5 -40. -42.5 -45. -47.5 -50. -52.5 -55. -57.5
-60. -62.5 -65. -67.5 -70. -72.5 -75. -77.5 -80. -82.5 -85. -87.5
-90. ] [ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 180. 182.5 185. 187.5 190. 192.5 195. 197.5
200. 202.5 205. 207.5 210. 212.5 215. 217.5 220. 222.5
225. 227.5 230. 232.5 235. 237.5 240. 242.5 245. 247.5
250. 252.5 255. 257.5 260. 262.5 265. 267.5 270. 272.5
275. 277.5 280. 282.5 285. 287.5 290. 292.5 295. 297.5
300. 302.5 305. 307.5 310. 312.5 315. 317.5 320. 322.5
325. 327.5 330. 332.5 335. 337.5 340. 342.5 345. 347.5
350. 352.5 355. 357.5]
Чтобы корректно отобразить данные необходимо перевести координаты. Это делается с помощью функции shiftgrid:
import mpl_toolkits.basemap as mtb
v[1,:,:],x = mtb.shiftgrid(180,v[1,:,:],x,start=False)
lon0 = 180 - конец или начало нового вектора долгот (см. start).
strat = False - заданная долгота lon0 (здесь равна 180) является концом нового вектора. Если True, то, соответственно, первым значением.
Карта выглядит неплохо, но что за белые линии (white stripes) вверху и снизу? Это следствие того, что вектор широт у наших данных идёт от 90 до -90, то есть с севера на юг. Проекция же корректно отображает данные, когда значение широты идут наоборот: от -90 до 90, то есть возрастают.
Чтобы наиболее безболезненно убрать этот дефект прочитаем данные задом наперед:
y = f.variables['lat'][::-1]
v = f.variables['pres'][:,::-1,:]
Для переменной v широты являются второй осью, поэтому конструкцию '::-1' нужно ставить туда.
Если приглядеться, то справа, в районе 180 градуса восточной долготы, заметна белая полосочка. Этот эффект есть разрыв данных на границе -180 и 180, что на самом деле есть одна долгота. На цилиндрической проекции это не очень заметно, но вот на другой:
m = Basemap(projection='ortho',lon_0=180,lat_0=30,resolution='l')
Здесь, в новой проекции, есть отличие от цилиндрических проекций, где задаются границы прямоугольника, lon_0 и lat_0 ( они задают координату центральной точки обзора).
Чтобы убрать этот разрыв есть специальная функция addcyclic:
from mpl_toolkits.basemap import addcyclic
v,x = addcyclic(v,x)
НО! В текущей версии эта функция работает только, если массив v - матрица, то есть имеет две оси. Так если вектор x имеет форму x.shape = (144,), то v.shape должен быть (n, 144).
В нашем случае это не так, у массива v есть третья ось. v.shape = (n, 144, 73).
На форуме я нашёл, как расширить возможности функции addcyclic.
Необходимо заменить строчки 4977-4998 в файле, который содержит код функции addcyclic. Пусть к файлу выглядит примерно так:
File "C:\WinPython-64bit-3.3.2.2\python-3.3.2.amd64\lib\site-ackages\mpl_toolkits\basemap\__init__.py", line 4997, in addcyclic
Нужно заменить функцию addcyclic на новую:
# ------------- new addcyclic -------------------------------
def addcyclic(arrin,lonsin):
"""
``arrout, lonsout = addcyclic(arrin, lonsin)``
adds cyclic (wraparound) point in longitude to ``arrin`` and ``lonsin``,
assumes longitude is the right-most dimension of ``arrin``.
"""
nlons = arrin.shape[-1]
newshape = list(arrin.shape)
newshape[-1] += 1
if ma.isMA(arrin):
arrout = ma.zeros(newshape,arrin.dtype)
else:
arrout = np.zeros(newshape,arrin.dtype)
arrout[...,0:nlons] = arrin[:]
arrout[...,nlons] = arrin[...,0]
if ma.isMA(lonsin):
lonsout = ma.zeros(nlons+1,lonsin.dtype)
else:
lonsout = np.zeros(nlons+1,lonsin.dtype)
lonsout[0:nlons] = lonsin[:]
lonsout[nlons] = lonsin[-1] + lonsin[1]-lonsin[0]
return arrout,lonsout
# -----------------------------------------------------------
Теперь addcyclic будет работать с n-мерными массивами. Правда ограничения на форму массивов всё же остались. Вектор, который будет замыкаться (x здесь) должен быть ПОСЛЕДНЕЙ осью вектора v. То есть (12,73,144) можно, а (12,144,73) нельзя замыкать с вектором x(144,). В нашем случае привнесенных изменений достаточно. Рисуем:
Осталось добавить цветовую шкалу.
cs = m.contourf(lons,lats,v[1,:,:],50)
cbar = plt.colorbar(cs,orientation = 'horizontal')
cbar.set_label('SLP, hPa')
или так, задав пользовательские границы изолиний и шкалы:
levels = (500,700,800,900,925,950,960,970,980,990,1000,1010,1020,1030,1040,1050,1100)
cs = m.contourf(lons,lats,v[1,:,:],levels)
cbar = plt.colorbar(cs,ticks=levels,orientation = 'horizontal')
cbar.set_label('SLP, hPa')
И весь код целиком:
# ------------------- Basemap 3 --------------------------------------
import mpl_toolkits.basemap as mtb
from scipy.io import netcdf as nc
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap,addcyclic
fname = 'pres.mon.ltm.nc'
f = nc.netcdf_file(fname)
fkeys = list(f.variables.keys())
print(fkeys)
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][::-1]
v = f.variables['pres'][:,::-1,:]
v1 = f.variables['pres']
f.close()
v = v*v1.scale_factor + v1.add_offset
v,x = addcyclic(v,x)
plt.figure(figsize=(15,15))
m = Basemap(projection='mill')
v[1,:,:],x = mtb.shiftgrid(180,v[1,:,:],x,start=False)
m.drawcoastlines(linewidth=1)
m.drawcountries(linewidth=0.5)
m.drawmapboundary()
parallels = np.arange(-90,90+1,30)
meridians = np.arange(0,360+1,60)
m.drawparallels(parallels,labels=[True,True,True,True])
m.drawmeridians(meridians,labels=[True,False,False,True])
x1,y1 = np.meshgrid(x,y)
lons,lats = m(x1,y1)
cs = m.contourf(lons,lats,v[1,:,:],50)
cbar = plt.colorbar(cs,orientation = 'horizontal')
cbar.set_label('SLP, hPa')
plt.title('Test pressure [hPa] map')
plt.savefig('Test_basemap_map.png')
plt.show()
# ----------------------------------------------------------------------
В глаза сразу бросается 90 градусов вверху графика. Также наша роза не является замкнутой. Необходимо привести это изображение в порядок. Для этого добавим следующие строчки в код:
ax.set_theta_direction(-1) # теперь увеличение угла theta будет по часовой (-1 - увеличение угла
# идёт против часовой стрелки)
ax.set_theta_offset(pi/2.0) # изменим положение нуля (ставим его на положение 'N', которое
# сейчас равно 90 градусам)
ax.plot((angles[-1],angles[0]),(r[-1],r[0]), color='r', linewidth=1.) # Добавляем ещё один plot-линию,
# которая соединяет последнюю и первую точки.
Сохраняем получившийся рисунок в файл с расширением 'png' (можно и другой) с помощью
plt.savefig('rose.png')
Вот так намного лучше! А теперь выведем сразу две розы:
#----------- Matplotlib 2-------------------------------------
from math import pi
import numpy as np
import matplotlib.pyplot as plt
angles = np.arange(0, 2*pi, pi/4.)
fig = plt.figure(figsize=(10,5)) # зададим явно размер полотна, чтобы графики не перекрывались
for i in range(1,3,1):
r = np.random.random(len(angles)) # для каждой розы значения будут разные
ax = fig.add_subplot(120+i, projection='polar')
ax.plot(angles, r, color='r', linewidth=1.)
ax.plot((angles[-1],angles[0]),(r[-1],r[0]), color='r', linewidth=1.)
ax.set_theta_direction(-1)
ax.set_theta_offset(pi/2.0)
ax.set_title("Wind roses "+str(i), loc = 'center')
plt.show()
plt.savefig('rose_subplot.png')
#--------------------------------------------------------------
SUBPLOT
ax = fig.add_subplot(121, projection='polar')
121 - rcn: r = 1, c = 2, n = 1
r - число строк в таблице
c - число столбцов в таблице
n - номер ячейки в получившейся таблице
Тогда, чтобы нарисовать графики друг под другом нужно будет задать значения для subplot 211 и 212. В приведенном выше коде этот номер изменяется с помощью переменной цикла i.
4. Basemap. Это дополнение к matplotlib, которое позволяет рисовать географические карты. Основу работы basemap составляет картографическая проекция, на которой будет изображаться то или иное поле. Разберём основные возможности Basemap на примере данных давления из реанализа NCEP/NCAR, которые мы обрабатывали в обзоре Scipy.
С помощью модуля basemap мы импортируем из 'mpl_toolkits.basemap' класс Basemap.
Зададим в качестве основы проекцию Миллера. Всего на выбор доступно более 30 проекций. Нарисуем карту земного шара в цилиндрической проекции, которой является проекция Миллера.
# ------------------- Basemap 1 --------------------------------------
import numpy as np
from mpl_toolkits.basemap import Basemap
m = Basemap(projection='mill')
m.drawcoastlines(linewidth=1) # Рисуем береговую линию
m.drawcountries(linewidth=0.5) # Рисуем границы стран
m.drawmapboundary() # Рисуем рамку вокруг карты
# Рисуем географическую сетку (долготы и широты)
parallels = np.arange(-90,90+1,30)
meridians = np.arange(0,360+1,60)
m.drawparallels(parallels,labels=[True,True,True,True])
m.drawmeridians(meridians,labels=[True,False,False,True])
# -----------------------------------------------------------------------
Получится следующее:
Раскрасим карту и добавим реки на карту:
m.drawmapboundary(fill_color='aqua') # Делаем заливку океана голубым цветом
m.fillcontinents(color='orange',lake_color='green') # Красим континенты в оранжевый,
# а озера - в зелёный.
m.drawrivers() # Добавляем реки
За дополнительными возможностями по рисованию картографической основы - сюда. Разобравшись с основами, перейдём к отображению информации на карте.
# ------------------- Basemap 2 --------------------------------------
from scipy.io import netcdf as nc
import numpy as npС помощью модуля basemap мы импортируем из 'mpl_toolkits.basemap' класс Basemap.
from mpl_toolkits.basemap import Basemap
Зададим в качестве основы проекцию Миллера. Всего на выбор доступно более 30 проекций. Нарисуем карту земного шара в цилиндрической проекции, которой является проекция Миллера.
# ------------------- Basemap 1 --------------------------------------
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
m = Basemap(projection='mill')
m.drawcoastlines(linewidth=1) # Рисуем береговую линию
m.drawcountries(linewidth=0.5) # Рисуем границы стран
m.drawmapboundary() # Рисуем рамку вокруг карты
# Рисуем географическую сетку (долготы и широты)
parallels = np.arange(-90,90+1,30)
meridians = np.arange(0,360+1,60)
m.drawparallels(parallels,labels=[True,True,True,True])
m.drawmeridians(meridians,labels=[True,False,False,True])
# -----------------------------------------------------------------------
Получится следующее:
Раскрасим карту и добавим реки на карту:
m.drawmapboundary(fill_color='aqua') # Делаем заливку океана голубым цветом
m.fillcontinents(color='orange',lake_color='green') # Красим континенты в оранжевый,
# а озера - в зелёный.
m.drawrivers() # Добавляем реки
За дополнительными возможностями по рисованию картографической основы - сюда. Разобравшись с основами, перейдём к отображению информации на карте.
# ------------------- Basemap 2 --------------------------------------
from scipy.io import netcdf as nc
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
fname = 'pres.mon.ltm.nc'
f = nc.netcdf_file(fname)
fkeys = list(f.variables.keys())
print(fkeys)
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][:]
v = f.variables['pres'][:]
v1 = f.variables['pres']
f.close()
v = v*v1.scale_factor + v1.add_offset
m = Basemap(projection='mill')
m.drawcoastlines(linewidth=1)
m.drawcountries(linewidth=0.5)
m.drawmapboundary()
parallels = np.arange(-90,90+1,30)
meridians = np.arange(0,360+1,60)
m.drawparallels(parallels,labels=[True,True,True,True])
m.drawmeridians(meridians,labels=[True,False,False,True])
# Данные по широтам и долготам, взятые из netcdf файла необходимо перевести из векторов
# в матрицы. Делается это с помощью функции meshgrid
x1,y1 = np.meshgrid(x,y)
lons,lats = m(x1,y1) # задаем на картографической основе преобразованные в матрицы
# координаты
'''
Другой синтаксис перевода координат x и y в координаты, понятные карт. основе
lons,lats = m(*np.meshgrid(x,y))
'''
m.contourf(lons,lats,v[1,:,:],40) # Рисуем на картографической основе m данные за февраль
# contourf рисует изолинии с послойной окраской через 40 единиц по сеточным данным
# Массив v[1,:,:] - срез от всего массива v, представляющий данные за один месяц февраль
plt.title('Test map N') # Название карты
plt.savefig('Test_basemap_map.png')
plt.show()
# ------------------------------------------------------------------------------
Запускаем скрипт и получаем ... вот это!
1) Нарисовано только половина карты!
2) Непонятные белые полосы сверху и снизу карты!
Начнём с главного. Код рисует половину карты потому что данная проекция рассчитана на изменение долгот в пределах от -180 до 180, а в наших данных реанализа долготы изменяются в других пределах, а именно от 0 до 357.5:
>>> y,x
[ 90. 87.5 85. 82.5 80. 77.5 75. 72.5 70. 67.5 65. 62.5
60. 57.5 55. 52.5 50. 47.5 45. 42.5 40. 37.5 35. 32.5
30. 27.5 25. 22.5 20. 17.5 15. 12.5 10. 7.5 5. 2.5
0. -2.5 -5. -7.5 -10. -12.5 -15. -17.5 -20. -22.5 -25. -27.5
-30. -32.5 -35. -37.5 -40. -42.5 -45. -47.5 -50. -52.5 -55. -57.5
-60. -62.5 -65. -67.5 -70. -72.5 -75. -77.5 -80. -82.5 -85. -87.5
-90. ] [ 0. 2.5 5. 7.5 10. 12.5 15. 17.5 20. 22.5
25. 27.5 30. 32.5 35. 37.5 40. 42.5 45. 47.5
50. 52.5 55. 57.5 60. 62.5 65. 67.5 70. 72.5
75. 77.5 80. 82.5 85. 87.5 90. 92.5 95. 97.5
100. 102.5 105. 107.5 110. 112.5 115. 117.5 120. 122.5
125. 127.5 130. 132.5 135. 137.5 140. 142.5 145. 147.5
150. 152.5 155. 157.5 160. 162.5 165. 167.5 170. 172.5
175. 177.5 180. 182.5 185. 187.5 190. 192.5 195. 197.5
200. 202.5 205. 207.5 210. 212.5 215. 217.5 220. 222.5
225. 227.5 230. 232.5 235. 237.5 240. 242.5 245. 247.5
250. 252.5 255. 257.5 260. 262.5 265. 267.5 270. 272.5
275. 277.5 280. 282.5 285. 287.5 290. 292.5 295. 297.5
300. 302.5 305. 307.5 310. 312.5 315. 317.5 320. 322.5
325. 327.5 330. 332.5 335. 337.5 340. 342.5 345. 347.5
350. 352.5 355. 357.5]
Чтобы корректно отобразить данные необходимо перевести координаты. Это делается с помощью функции shiftgrid:
import mpl_toolkits.basemap as mtb
v[1,:,:],x = mtb.shiftgrid(180,v[1,:,:],x,start=False)
lon0 = 180 - конец или начало нового вектора долгот (см. start).
strat = False - заданная долгота lon0 (здесь равна 180) является концом нового вектора. Если True, то, соответственно, первым значением.
Карта выглядит неплохо, но что за белые линии (white stripes) вверху и снизу? Это следствие того, что вектор широт у наших данных идёт от 90 до -90, то есть с севера на юг. Проекция же корректно отображает данные, когда значение широты идут наоборот: от -90 до 90, то есть возрастают.
Чтобы наиболее безболезненно убрать этот дефект прочитаем данные задом наперед:
y = f.variables['lat'][::-1]
v = f.variables['pres'][:,::-1,:]
Для переменной v широты являются второй осью, поэтому конструкцию '::-1' нужно ставить туда.
m = Basemap(projection='ortho',lon_0=180,lat_0=30,resolution='l')
Здесь, в новой проекции, есть отличие от цилиндрических проекций, где задаются границы прямоугольника, lon_0 и lat_0 ( они задают координату центральной точки обзора).
Чтобы убрать этот разрыв есть специальная функция addcyclic:
from mpl_toolkits.basemap import addcyclic
v,x = addcyclic(v,x)
НО! В текущей версии эта функция работает только, если массив v - матрица, то есть имеет две оси. Так если вектор x имеет форму x.shape = (144,), то v.shape должен быть (n, 144).
В нашем случае это не так, у массива v есть третья ось. v.shape = (n, 144, 73).
На форуме я нашёл, как расширить возможности функции addcyclic.
Необходимо заменить строчки 4977-4998 в файле, который содержит код функции addcyclic. Пусть к файлу выглядит примерно так:
File "C:\WinPython-64bit-3.3.2.2\python-3.3.2.amd64\lib\site-ackages\mpl_toolkits\basemap\__init__.py"
Spyder при запуске кода с базовой версией addcyclic выдаст ошибку, и можно будет щелкнуть по гиперссылке, чтобы этот файл открылся прям в Spyder. Ошибка будет такая:File "C:\WinPython-64bit-3.3.2.2\python-3.3.2.amd64\lib\site-ackages\mpl_toolkits\basemap\__init__.py", line 4997, in addcyclic
Нужно заменить функцию addcyclic на новую:
# ------------- new addcyclic -------------------------------
def addcyclic(arrin,lonsin):
"""
``arrout, lonsout = addcyclic(arrin, lonsin)``
adds cyclic (wraparound) point in longitude to ``arrin`` and ``lonsin``,
assumes longitude is the right-most dimension of ``arrin``.
"""
nlons = arrin.shape[-1]
newshape = list(arrin.shape)
newshape[-1] += 1
if ma.isMA(arrin):
arrout = ma.zeros(newshape,arrin.dtype)
else:
arrout = np.zeros(newshape,arrin.dtype)
arrout[...,0:nlons] = arrin[:]
arrout[...,nlons] = arrin[...,0]
if ma.isMA(lonsin):
lonsout = ma.zeros(nlons+1,lonsin.dtype)
else:
lonsout = np.zeros(nlons+1,lonsin.dtype)
lonsout[0:nlons] = lonsin[:]
lonsout[nlons] = lonsin[-1] + lonsin[1]-lonsin[0]
return arrout,lonsout
# -----------------------------------------------------------
Теперь addcyclic будет работать с n-мерными массивами. Правда ограничения на форму массивов всё же остались. Вектор, который будет замыкаться (x здесь) должен быть ПОСЛЕДНЕЙ осью вектора v. То есть (12,73,144) можно, а (12,144,73) нельзя замыкать с вектором x(144,). В нашем случае привнесенных изменений достаточно. Рисуем:
cs = m.contourf(lons,lats,v[1,:,:],50)
cbar = plt.colorbar(cs,orientation = 'horizontal')
cbar.set_label('SLP, hPa')
levels = (500,700,800,900,925,950,960,970,980,990,1000,1010,1020,1030,1040,1050,1100)
cs = m.contourf(lons,lats,v[1,:,:],levels)
cbar = plt.colorbar(cs,ticks=levels,orientation = 'horizontal')
cbar.set_label('SLP, hPa')
# ------------------- Basemap 3 --------------------------------------
import mpl_toolkits.basemap as mtb
from scipy.io import netcdf as nc
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap,addcyclic
fname = 'pres.mon.ltm.nc'
f = nc.netcdf_file(fname)
fkeys = list(f.variables.keys())
print(fkeys)
t = f.variables['time'][:]
x = f.variables['lon'][:]
y = f.variables['lat'][::-1]
v = f.variables['pres'][:,::-1,:]
v1 = f.variables['pres']
f.close()
v = v*v1.scale_factor + v1.add_offset
v,x = addcyclic(v,x)
plt.figure(figsize=(15,15))
m = Basemap(projection='mill')
v[1,:,:],x = mtb.shiftgrid(180,v[1,:,:],x,start=False)
m.drawcoastlines(linewidth=1)
m.drawcountries(linewidth=0.5)
m.drawmapboundary()
parallels = np.arange(-90,90+1,30)
meridians = np.arange(0,360+1,60)
m.drawparallels(parallels,labels=[True,True,True,True])
m.drawmeridians(meridians,labels=[True,False,False,True])
x1,y1 = np.meshgrid(x,y)
lons,lats = m(x1,y1)
cs = m.contourf(lons,lats,v[1,:,:],50)
cbar = plt.colorbar(cs,orientation = 'horizontal')
cbar.set_label('SLP, hPa')
plt.title('Test pressure [hPa] map')
plt.savefig('Test_basemap_map.png')
plt.show()
# ----------------------------------------------------------------------
