Вырезка области из данных на криволинейной сетке при помощи 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()
У меня имеются координаты нужной области (с 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()
Дальше нужно наложить на карту координаты криволинейной сетки, которые находятся у меня в переменных 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()
получил такую картину:
Осталось только подогнать значения для вырезки (когда смотришь на карту, делать это в разы легче). В итоге я получил такой код
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()
и такую карту:
Соответственно, искомые границы будут [350:770,235:500]