«Но я лентяй! Я ненавижу работать! Особенно я ненавижу тяжёлый труд во всех его проявлениях! Хитрые и изящные решения — вот мой конёк!»

― Элиезер Юдковский, «Гарри Поттер и методы рационального мышления»

«Торжественно клянусь, что замышляю шалость, и только шалость!»

Что, если я скажу вам, что программирование для современного ученого и инженера — это не просто киллер-фича в резюме, а спасательный трос в мире больших данных, рутинных процессов и нестандартных задач. Именно этой идеи придерживается наше сообщество, одним из ключевых принципов которого является принцип «DIY» (Do it yourself — сделай сам). Идея этого принципа проста: в решении своих повседневных задач мы полагаемся не на монструозные коммерческие костыли в красивой обертке готовые решения, а оптимизируем свои решения с помощью копролитов и палок Python и данных.

Каждый день, от проекта к проекту, мы сталкивается с лавиной данных, которые зачастую невозможно быстро и качественно обработать стандартными инструментами. Готовые программные пакеты хороши для типовых задач, но они бессильны, когда нужно решить что-то нестандартное — будь то автоматическая обработка сотен скважин или поиск скрытых закономерностей в данных. Программирование преодолевает эти ограничения: оно превращает рутину в несколько строк кода, а нестандартные задачи — в вызов, который вы в силах принять.

Цель данной заметки не научить вас программированию, нет. Для этого существует множество отличных источников. Наша цель — дать вам пинок под зад краткое руководство к действию. На примере реальных данных показать, что программирование на Python это максимально просто. Просто настолько, что после прочтения этой заметки даже ваша бабушка пересядет со спиц в Jupyter Notebook.

Если после этого текста у тебя не зачесались руки написать хотя бы строчку кода — проверь пульс. Ты точно ещё жив? Если пульс есть, но желания так и не появилось — сгоняй нам за пивком, а мы пока начинаем.

От таблиц к пространственным данным

Прежде чем мы погрузимся в увлекательный мир программирования, хотелось бы отметить, что при подготовке данной заметки, были использованы геолого-геофизические данные, распространяемые с курсом практических упражнений по сейсмической интерпретации под авторством А.С. Кирилова и К.Е. Закревского. Со своей стороны, SciPunk выражает респект авторам курса за их работу и приверженность принципам свободного распространения обучающих материалов. С деятельностью авторов курса можно ознакомиться на их сайте или сообществе Вконтакте. А теперь стартуем!

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

Начнем с простого: импортируем таблицу со стратиграфическими отбивками из текстового файла формата CSV. Для работы с таблицами в Python есть замечательная библиотека Pandas. Для того, чтобы ей воспользоваться, необходимо для начала ее импортировать:

1

import pandas as pd # Импортируем библиотеку Pandas

Далее, прочитаем содержимое файла TOPs.dat с использованием функции read_csv() из библиотеки Pandas. Для этого нам потребуется всего одна строчка кода. Все что нужно сделать — передать в качестве аргументов в функцию read_csv() несколько параметров:

1. Путь до нашего файла 'data\wells\TOPs.dat'
2. Символ, использующийся в качестве разделителя в фале csv. В нашем случае это символ пробела sep = ' '
3. Дополнительно используем параметр skipinitialspace = True для игнорирования лишних пробелов

1
2
3

#Чтение данных из файла CSV
well_tops = pd.read_csv('data\wells\TOPs.dat', sep = ' ', skipinitialspace = True)
well_tops

Вуаля! Теперь в переменной well_tops хранится объект DataFrame с нашей таблицей стратиграфических отбивок.

Not Bad! С этим уже можно работать, но мы пойдем дальше. Отличительной особенностью геолого-геофизической информации является ее пространственная привязка. Для тех, кто забыл что это такое, ниже приведена база небольшая информационная справка.

[!info] Геопространственные данные Данные, которые связывают некоторую геологическую информацию с конкретным географическим положением, называются геопространственными, и являются основой геоинформационных систем.

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

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

Для работы с картографической системой координат важно знать, какая проекция использовалась для отображения пространственных объектов на плоскость. Существует множество различных проекций, каждая из которых подходит под определенные участки Земной поверхности, но ни одна из них не может точно описать всю поверхность, не внося некоторых искажений.

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

— Дед, ты бредишь? В таблице же есть столбцы с координатами! — скажете вы мне. — Yep — отвечу вам я.

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

Для описания геометрии пространственных типов данных в геоинформатике используются геометрические примитивы: точка Point, отрезок LineString, многоугольник Polygon. Примитивы представляют собой набор координат, характеризующих положение этих данных в одной из систем координат. В нашем случае стратиграфические отбивки представляют собой точечные данные, для описания которых подходит геометрический примитив Point.

Теперь давайте модифицируем нашу таблицу так, чтобы она отвечала всем требованиям пространственного типа данных. С этим нам поможет библиотека GeoPandas. По сути это те же яйца, только в профиль тот же Pandas, только с поддержкой всех особенностей пространственных данных. Потребуется еще две строчки кода (хотя можно было бы уложиться и в одну), не считая импорта библиотеки.

1
2
3
4
5
6

import geopandas as gpd

#Конвертируем информацию в пространственный тип данных GeoDataFrame (POINT)
geometry = gpd.points_from_xy(well_tops.X, well_tops.Y)
gdf_tops = gpd.GeoDataFrame(well_tops, crs = 'EPSG:32640', geometry = geometry)
gdf_tops

Все что нужно было сделать — это создать специальный объект с описанием геометрии наших данных с помощью функции points_from_xy(), на вход которой мы подали два столбца с координатами well_tops.X и well_tops.Y из нашей исходной таблицы. Затем, на основе нашего DataFrame из переменной well_tops мы создали WunderGeoDataFrame, передав дополнительно переменную с геометрией geometry и информацию о системе координат и проекции crs = 'EPSG:32640' через ее уникальный идентификационный код EPSG. Смотрим что получили:

Ничего не напоминает по структуре? Правильно, точечный shp-файл. Теперь столбцы Well, Surface, X, Y, Z — это атрибутивная информация, а geometry — столбец, который хранит геометрию наших точек aka отбивок. А чтобы окончательно вас в этом убедить, давайте сохраним полученный GeoDataFrame в наш первый самопальный shp-файл, вот так:

1
2

# Сохраняем отбивки в формате shp-файла
gdf_tops.to_file(driver='ESRI Shapefile', filename=r'well_tops.shp', crs='EPSG:32640')

Как вы уже могли догадаться, читать shp-файлы с помощью библиотеки GeoPandas не сильно сложнее. Просто передаем путь и имя нашего файла в функцию read_file():

1
2

# Читаем наш новоиспеченный shp-файл
gdf_tops = gpd.read_file('well_tops.shp')

Теперь посмотрим на данные в таблице. Предположим, нам нужно получить список уникальных значений столбца Surface — легко:

1

gdf_tops['Surface'].unique()

А теперь мы хотим визуализировать отбивки артинского яруса P1ar по всем скважинам. Для начала отфильтруем данные по нашему запросу, а результат сохраним в переменную p1ar_tops:

1

p1ar_tops = gdf_tops[gdf_tops['Surface'] == 'P1ar'].copy()

Затем, отобразим карту с отбивками используя встроенный метод plot(). При этом маркер отбивки покрасим цветом, согласно значению ее абсолютной отметки column = 'Z':

1
2
3
4
5

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
p1ar_tops.plot(column = 'Z', legend = True, cmap = 'Spectral_r', ax = ax)
ax.grid(True, linestyle = 'dotted')

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

Двумерные данные. Поверхности (гриды)

Не одними лишь таблицами живет геолог(иня)\геофизик(чка). Не менее популярный формат представления геолого-геофизической информации — двумерные поверхности (гриды): структурные карты, карты толщин, карты изохрон, матрицы геофизических полей (гравика, магнитка и пр.) и так далее. Их тоже хотелось бы читать из файла. Ну так давайте попробуем. Возьмем один из распространенных открытых форматов Landmark Z-Map Grid. Здесь тоже есть готовые решения, например библиотека Zmapio. Идем по ссылке, читаем документацию, повторяем:

1
2

from zmapio import ZMAPGrid
time_grid = ZMAPGrid(r'data\surfaces\P1ar_time.dat')

Чтобы нашим отбивкам не было одиноко, подгрузим поверхность изохрон (времен пробега волны) отражающего горизонта P1ar, полученную по данным сейсморазведки МОГТ-3D. Импортировали — теперь посмотрим на нее, используя встроенный метод plot():

1
2
3
4
5

fig, ax = plt.subplots()
time_grid.plot(axes=ax, cmap='Spectral_r')
ax.axis('scaled')
ax.grid(True, linestyle='dotted')
plt.colorbar()

Отлично! Визуализация это конечно здорово, но нам ведь с этими данными еще науку делать работать придется по мере углубления в материал (хотя, откровенно говоря, правильная визуализация данных — уже 70% успеха). Посмотрим, в каком виде zmapio может вернуть нам нашу поверхность. Идем в документацию, смотрим — есть два стула варианта:

Вариант 1: Можно получить двумерный numpy массив используя метод z_values, вот так: time_grid.z_values Вариант 2: Работать с поверхностью, как с базой данных, сконвертировав ее в Pandas DataFrame, вот так: time_grid.to_dataframe()

В первом случае нам придется обращаться к данным через нумерованные индексы, которые каким-то образом придется постоянно соотносить с координатами. Во втором - работать с плоским массивом данных, но с возможностью поименованного (непосредственно по координатам) обращения к данным. Такая вот вилка в глаз Мортона (или аксиома Эскобара, кому как угодно).

Нам оно не надо. Благо братики постряпали удобную библиотеку Xarray, которая лишит нас этой головной боли. Xarray — это сын маминой подруги в мире многомерных массивов Python. Простыми словами, Xarray — многомерный Pandas, который позволяет работать с массивами, используя поименованные индексы. Не переживайте, если пока не понятно, сейчас будем разбираться.

Итак, для удобства последующей работы создадим двумерный xarray из numpy массива нашей поверхности изохрон, предварительно импортировав необходимую библиотеку:

1
2
3
4
5
6
7

import xarray as xr

time_grid_xr = xr.DataArray(time_grid.z_values.T,
                            name='time',
                            coords=[time_grid.y_values[:, 0],
                                    time_grid.x_values[0]],
                            dims=['y', 'x'])

Заглянем внутрь полученного массива xarray, вызвав нашу переменную time_grid_xr:

Что имеем? Внутри нашего xarray есть двумерный массив времен пробега волны и два массива с координатами X и Y, используя которые мы можем обращаться к нашим данным. Например, теперь мы легко можем получить значения времен в точках скважин (по их координатам) или снять значения времен вдоль линии сейсмического профиля. Но прежде чем мы это сделаем, давайте отобразим наши данные на одной карте:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

fig, ax = plt.subplots()

# Визуализируем цветовую расскраску карты по уровням с шагом 10 мс.
time_grid_xr.plot.contourf(ax=ax, cmap='Spectral_r',
                           levels = np.arange(-1700, -1900, -10))
# Визуализируем изолии времен с шагом 20 мс.
time_grid_xr.plot.contour(colors='black', ax=ax, linewidths=0.5,
                          levels = np.arange(-1700, -1900, -20),
                          linestyles='solid')
# Визуализируем отбивки по горизонту P1ar
p1ar_tops.plot(ax=ax, markersize = 25, label = 'P1ar tops', zorder=10)

ax.legend() # Добавлям легенду нашей карты
ax.grid(True, linestyle='dotted') # Добавлям координатную сетку

Здесь мы применили чуть более подробное оформление нашего графика:

Возвращаемся к нашим индексам. Если бы мы работали с поверхностью как с массивом numpy, мы должны были бы помнить, что “координаты” массива numpy (индексы) всегда начинаются со значения 0 и идут с шагом 1. Предположим: мы хотим получить значения времени пробега волны с карты изохрон в точке скважины с координатой X и Y. Для того, чтобы это сделать, мы должны были бы сначала определить порядковый номер (индекс) этой точки в системе координат numpy массива, и только по этому индексу получить значение с карты. Сложно, согласитесь? Вместо этого мы предусмотрительно сконвертировали наш numpy массив в формат xarray и теперь можем обращаться к данным напрямую по значению координат, а вся головная боль будет происходить под капотом библиотеки xarray. В конце концов мы с вами геологи\геофизики, а не пограмисты за 300кк в наносек — нам и своих проблем хватает. Итак, давайте таки получим значения времен в точках скважин с карты изохрон и запишем их в новую колонку times таблицы с отбивками по артинскому ярусу p1ar_tops:

1
2
3
4
5
6
7
8

# Подготовим два массива с координатами скважин в формате xarray
x_points = xr.DataArray(p1ar_tops.X.values, dims='points')
y_points = xr.DataArray(p1ar_tops.Y.values, dims='points')

# Формируем запрос для получения значений времен из массива xarray по координатам
time_points = time_grid_xr.sel(x = x_points, y = y_points, method='nearest')
# Записываем полученные значения времен в новый столбец таблицы p1ar_tops
p1ar_tops['times'] = time_points.values

Смотрим что получили:

Атлична-атлична! А теперь построим график зависимости “время-глубина” (диаграмма рассеяния aka кросс-плот) в точках скважин, например. Полученные уравнения регрессии братья сейсмики иногда используют для пересчета своих карт времен (изохрон) в карты глубин (структурные карты).

1
2
3
4
5
6

import seaborn as sns

sns.lmplot(x='times', y='Z', data=p1ar_tops)
plt.grid(True, linestyle='dotted')
plt.title('Time vs. Depth (P1ar)')
plt.show()

В этот раз для визуализации графика используем библиотеку Seabron. Это библиотека визуализации красивых статистических графиков, основанная на библиотеке matplotlib.

Красивое… Ну что ж, легкой поступью мы с вами подошли к теме анализа данных. Но не будем забегать вперед! Оставим эту тему на следующий урок по DIY. Обязательно научимся оценивать уравнения линейной регрессии, займемся прогнозами (в том числе с помощью машинного обучения) и оценкой неопределенностей. А пока продолжаем разбираться со входными данными.

Данные сейсморазведки. Сейсмические разрезы формата SEG-Y

Карту изохрон мы разобрали, теперь остановимся на данных сейсморазведки, по которым эта карта была построена. Начнем с временных сейсмических разрезов МОГТ-2D. Стандартом хранения сейсмических данных является бинарный формат SEG-Y, разработанный международным обществом геофизиков-разведчиков — Society of Exploration Geophysicists (SEG). При желании, со спецификациями данного формата можно ознакомиться на сайте сообщества.

Благо об этом уже кто-то позаботился за нас, ознакомился со спецификациями к данному типу файлов и любезно предоставил библиотеку для работы с ними. На самом деле таких библиотек на данный момент существует несколько. ObsPy — одна из таких библиотек. Она ориентирована на обработку сейсмологических данных и предоставляет возможность чтения данных в формате SEG-Y. В данном уроке мы воспользуемся возможностями именно этой библиотеки.

[!attention] Как известно, временной сейсмический разрез представляет собой набор так называемых трасс (временных рядов), которые получаются путем суммирования сейсмограмм общей глубинной точки (ОГТ). Поскольку библиотека ObsPy ориентирована, в первую очередь, на работу с данными сейсмологии, то с сейсмическими разрезами она работает именно в контексте временных рядов. Это может показаться не совсем удобно для работы с 2D разрезами и совершенно непригодно для данных 3D, но для понимания структуры формата SEG-Y данная библиотека подходит как нельзя лучше. Более совершенные библиотеки для работы с данными сейсморазведки будут разобраны в следующих Dungeon master-классах

От слов к делу: прочитаем данные SEG-Y файла (временного сейсмического разреза МОГТ-2D). Используем функцию _read_segy из библиотеки ObsPy:

1
2
3
4

from obspy.io.segy.segy import _read_segy

stream = _read_segy(r'data\seismic\22.segy', headonly=True) # Читаем
stream # Смотрим что получили

Функция возвращает нам объект, который содержит в себе 958 трасс, о чем нам любезно сообщает интерпретатор Python: 958 traces in the SEG Y structure. К конкретной трассе можно обратиться по ее индексу (начиная с нулевого). Так и поступим, заодно сразу ее визуализируем:

1
2
3
4
5
6
7

# Выдернем нулевую трассу из 958
one_trace = stream.traces[0] # Крендель в квадратных скобках и есть индекс

plt.figure(figsize=(16,1))
plt.plot(one_trace.data)
plt.grid(True, linestyle='dotted')
plt.show()

Так и есть - трасса. По оси Y здесь амплитуда сигнала, а вот с осью X есть нюанс - вместо привычного двойного времени пробега волны (TWT) здесь отражены дискреты сигнала. Чтобы это исправить, нам нужно получить больше информации о сигнале, как минимум шаг его дискретизации.

— А где его взять?!!11одинодин — спросите вы.

Дополнительная информация о сигнале хранится в заголовках SEG-Y файла, их есть несколько:

1. Текстовый заголовок файла: stream.textual_file_header. Не обязателен к заполнению, нет четкого стандарта его заполнения, а посему часто бывает пустым или содержит нечто бессмысленное
2. Бинарный заголовок файла: stream.binary_file_header имеет четкую структуру и именно там можно найти то, что нас интересует
3. Кроме того, у каждой трассы есть свой заголовок, доступ к которому можно получить через stream.traces[0].header. Информация о шаге дискретизации тут тоже присутствует

Теперь давайте пофиксим нашу трассу с учетом вновь полученной информации:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

import numpy as np

# Шаг дискретизации сигнала
dt = stream.binary_file_header.sample_interval_in_microseconds / 1e+03
# Количество отчетов в трассе
samples = stream.binary_file_header.number_of_samples_per_data_trace
twt = np.arange(0, samples * dt, dt) # Массив двойного времени волны (TWT)

one_trace = stream.traces[0]

plt.figure(figsize=(16,1))
plt.plot(twt, one_trace.data)
plt.grid(True, linestyle='dotted')
plt.show()

Теперь все верно. Длина записи сигнала в данном случае составляет 3.6 секунды. Осталось получить разрез целиком. Поскольку ObsPy работает с трассами, то пройдемся по ним в цикле, а результат сохраним в наш любимый xarray:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# Собираем разрез из трасс с помощью цикла
data = np.stack(t.data for t in stream.traces)

seismic_xr = xr.DataArray(data.T,
                          name='amplitude',
                          coords=[-twt, np.arange(1, len(stream.traces) + 1)],
                          dims=['TWT', 'Pickets'])
                          
fig, ax = plt.subplots(figsize=(9,6))
seismic_xr.plot(vmin=-40, vmax=40, cmap='RdBu_r', ax=ax)
ax.grid(True, linestyle='dotted')

Супер! Осталось разобраться с геометрией профиля. Координаты каждой трассы можно получить через заголовки самих трасс stream.traces[0].header. В нашем случае, координаты были записаны здесь:

X: stream.traces[0].header.source_coordinate_x Y: stream.traces[0].header.source_coordinate_y

Кстати, никогда не стесняйтесь экспериментировать с кодом: заглядывайте внутрь объектов, смотрите из чего они состоят, обращайтесь к официальной документации и не забывайте яндексить гуглить (лучше всего на английском языке. Так сложилось, что вероятность получить ответ на свой вопрос в англоязычном сегменте интернета гораздо выше).

С местоположением координат внутри файла разобрались. Но не будем же мы считывать координаты каждого профиля по отдельности? Как-никак мы с вами тут про автоматизацию говорили. Говорили? Ну было и было…

Я фантазирую так: как будто мы написали скрипт, который шарится внутри заданной папки в поисках SEG-Y файлов, попутно собирая из них необходимую информацию (включая геометрию профилей), а на выходе выдает нам GeoDataFrame со схемой изученности данными сейсморазведки.

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

1. Написать функцию, которая будет осуществлять поиск файлов с заданным расширением файлов внутри заданной папки, а на выходе выдавать нам DataFrame с находками
2. Написать функцию, которая будет последовательно проходить по списку найденных SEG-Y файлов из DataFrame, считывать из них необходимую информацию, а на выходе формировать GeoDataFrame с геометрией сейсмических профилей

Звучит неплохо. Давайте начнем с первого пункта:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

import os

def scan_files(path, type_files = ('.sgy', '.segy')):
    
    path_files = [] # Сюда будем сохранять пути файлов
    name_files = [] # Cюда имена файлов
    
    # Рекурсивный поиск файлов в заданной директории 
    for foldername, subfolder, filename in os.walk(path):
	    # Отфильтруем лишние файлы
        sgy = filter(lambda x: x.endswith(type_files), filename)
        
        for filename in sgy: # Пройдемся циклом по оставшимся файлам
	        # Формируем полный путь до текущего файла
            current_file = foldername + '\\' + filename 
            path_files.append(current_file) # Сохраняем путь в список
            name_files.append(filename) # Сохраняем имя файла в список
            
    # Формируем DataFrame из списков и возращаем из функции         
    return pd.DataFrame({'File': name_files, 'Path': path_files})

Итак, мы написали функцию scan_files, которая реализует функционал из первой подзадачи (см. выше). Функция ожидает на вход два аргумента:

1. Путь в формате строки, где необходимо искать файлы — path
2. Тип файлов в формате кортежа Pyhton, которые мы собираемся искать

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

1
2

sgy_files = scan_files(r'data\seismic', type_files = ('.sgy', '.segy'))
sgy_files

Отлично! Теперь мы можем использовать цикл для перебора содержимого полученного DataFrame, и считать необходимую информацию из найденных файлов SEG-Y:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

geometry = [] # Список для координат профилей
file = [] # Список для имен файлов
correct_sgy = [] # Список полного пути до файлов
incorrect_sgy = [] # Список кривых файлов, которые не открылись

# Пройдемся циклом по строкам DataFrame найденных файлов
for row in sgy_files.iterrows():
    try:
        stream = _read_segy(row[1]['Path'], headonly=True)
    except:
        print ('Ошибка чтения файла: ' + str(row[1]['Path']))
        incorrect_sgy.append(row[1]['Path'])
    else:
        x_coord = [t.header.source_coordinate_x for t in stream.traces]
        y_coord = [t.header.source_coordinate_y for t in stream.traces]
        geometry_list = list(zip(x_coord, y_coord))
        geometry.append(geometry_list)
        correct_sgy.append(row[1]['Path'])
        file.append(row[1]['File'])
        print ('Файл обработан: ' + str(row[1]['Path']))
        
geometry = [LineString(line) for line in geometry]
seismic_gdf = gpd.GeoDataFrame({'File': file, 'Path': correct_sgy},
                               crs='EPSG:32640', geometry=geometry)
seismic_gdf

FUCK YEAH! Кажется у нас получилось. Конечно, я немного схалтурил и не стал оформлять код в функцию, а еще не вынул дополнительную важную информацию из заголовков, но моя основная задача заключалась в том, чтобы научить вас ловить рыбу, а не накормить ей. Think about it… Зато я добавил исключения try — except, чтобы скрипт не спотыкался на кривых SEG-Yях. В качестве домашнего задания оставляю вам доработку напильником данного скрипта с учетом вышесказанного и ваших смелых идей. Результаты домашки будем рады видеть в нашем telegram-канале, в комментариях под новостью о статье.

Кстати, GeoDataFrame с положением сейсмических профилей seismic_gdf можете сохранить в shp-файл, как мы делали это с отбивками (см. выше) и загрузить его в ваш ArcGIS QGIS проект.

А теперь предлагаю нанести все наши данные на одну обзорную схему:

1
2
3
4
5
6
7
8
9

fig, ax = plt.subplots()

time_grid_xr.plot.contourf(ax = ax, cmap = 'Spectral_r',
                          levels = np.arange(-1700, -1900, -10))
p1ar_tops.plot(ax = ax, markersize = 25, label = 'wells', zorder=10)
seismic_gdf.plot(ax = ax, color = 'black', label = 'seismic lines')

ax.legend()
ax.grid(True, linestyle='dotted')

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# Подготовим два массива с координатами профиля в формате xarray
x_picket = xr.DataArray(seismic_gdf.iloc[4].geometry.coords.xy[0], dims='points')
y_picket = xr.DataArray(seismic_gdf.iloc[4].geometry.coords.xy[1], dims='points')

# Формируем запрос для получения значений времен из массива xarray по координатам
time_section = time_grid_xr.sel(x = x_picket, y = y_picket, method='nearest')

fig, ax = plt.subplots(figsize=(9,6))
seismic_xr.plot(vmin=-40, vmax=40, cmap='RdBu_r', ax=ax)
ax.plot(time_section-35, '--', label='P1ar_time')
ax.legend()
ax.grid(True, linestyle='dotted')

Заключение

На сегодня все. Подведем итоги. Сегодня мы:

1. Научились читать табличные данные из файла CSV-формата c помощью Python и библиотеки Pandas
2. Вспомнили базу из геоинформатики и поговорили о геопространственных данных
3. Научились превращать скучные таблицы в пространственные с помощью библиотеки GeoPandas
4. Создали свой первый shp-файл
5. Загрузили свой первый shp-файл
6. Научились работать с поверхностями формата Z-Map
7. Познакомились с библиотекой Xarray, разобрались почему этот формат хорош для работы с гридами
8. Научились извлекать информацию из грида используя координаты
9. Вплотную подошли к основам анализа данных — линейной регрессии. Договорились вернуться к ней позже
10. Разобрались как устроен формат SEG-Y. Научились работать с сейсмическими разрезами с помощью библиотеки ObsPy
11. Автоматизировали нашу первую рутинную задачу — написали шныря по папкам, который умеет строить схему изученности данными сейсморазведки МОГТ-2D. Да, он все еще не идеален, так что не забываем про домашнее задание.

Основную цель, которую мы преследовали при подготовке вводного урока по «DIY» — показать, что навыки программирования сегодня являются своего рода суперсилой для современного ученого и инженера, а также дать заинтересованным хорошую отправную точку для экспериментов со своими данными и подготовить читателя к нашей постоянной рубрике «очумелые ручки».

Безусловно, не весь спектр данных был освещен в данной заметке. В следующих уроках мы обязательно разберемся с данными сейсморазведки МОГТ-3D и поговорим про работу с данными геофизических исследований в скважинах (ГИС).

NB: Если вы до сих пор тратите часы на монотонные клики в интерфейсе, значит, вы уже отстаёте. Пора брать контроль над данными в свои руки — и Python ваш лучший инструмент для этого. Не забывайте: в эпоху big data и искусственного «интеллекта» выживут только те, кто умеет говорить на языке машин.

«Шалость удалась!»