Вырезка области из данных на криволинейной сетке при помощи Python

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

Имеем: большой массив с модельными данными для большой области океана. Для исследования требуется лишь небольшой участок, скажем с 65 по 75° с.ш. и с 20°з.д. по 20° в.д.

Задача: вырезать нужный нам маленький район из большого.

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

Для нас не важно получить именно прямоугольник с нужными с координатами, поэтому интерполяция в новую координатную сетку не нужна. Т.е. мы можем вырезать район любой формы, но должны выполняться два условия:

1. Вырезанный участок должен полностью покрыть заданный прямоугольник

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

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

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

Импортруем Basemap:

from mpl_toolkits.basemap import Basemap
%pylab inline

строим карту:

m = Basemap(projection='merc',llcrnrlat=60,urcrnrlat=78,\
llcrnrlon=-40,urcrnrlon=40,resolution='l')
figsize(15,15)

m.drawcoastlines()
m.fillcontinents(color='gray',lake_color='white')

plt.show()

map Читать далее «Вырезка области из данных на криволинейной сетке при помощи Python»