Вырезка области из данных на криволинейной сетке при помощи 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 У меня имеются координаты нужной области (с 65 по 75° с.ш. и с 20°з.д. по 20° в.д.), я хочу нанести на карту прямоугольник с выбранными координатами. Здесь я опять выбрал наиболее лёгкое для меня решение, я построил линии координатной сетки только для нужных координат и получил нужный мне прямоугольник.

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')

m.drawparallels(np.arange(65,76,1),labels=[1,0,0,0], fontsize =14)
m.drawmeridians(np.arange(-20,25,5),labels=[0,0,0,1], fontsize =14)

plt.show()

mam_coordinats

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

Далее я наудачу вырезал кусок из данных и наложил его на нашу карту,

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')

m.drawparallels(np.arange(65,76,1),labels=[1,0,0,0], fontsize =14)
m.drawmeridians(np.arange(-20,25,5),labels=[0,0,0,1], fontsize =14)
m.pcolormesh(lon[300:650,325:700],lat[300:650,325:700],data[4,300:650,325:700],shading='flat',cmap=plt.cm.jet,latlon=True)

plt.show()

получил такую картину:

mam_coordinats_data

Осталось только подогнать значения для вырезки (когда смотришь на карту, делать это в разы легче). В итоге я получил такой код

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

m.drawcoastlines()
m.fillcontinents(color='gray',lake_color='white')
m.drawparallels(np.arange(65,76,1),labels=[1,0,0,0], fontsize =14)
m.drawmeridians(np.arange(-20,25,5),labels=[0,0,0,1], fontsize =14)

m.pcolormesh(lon[350:770,235:500],lat[350:770,235:500],data[4,350:770,235:500],shading='flat',cmap=plt.cm.jet,latlon=True)

plt.show()

и такую карту:
mam_coordinats_data_new Соответственно, искомые границы будут [350:770,235:500]

Leave a Reply