Сегодня вместо обсуждения геологических моделей мы рассмотрим пример их программирования в среде Jupyter Notebook на Python 3 и с библиотеками Pandas, NumPy, SciPy, XArray, Dask Distributed, Numba, VTK, PyVista, Matplotlib. Это достаточно простой ноутбук с поддержкой многопоточной работы и возможностью запуска локально и в кластере для обработки больших данных, ленивых вычислений (ленивых) и наглядной 3D визуализации результатов.
Фактически, я попытался взять целый ряд сложных технических концепций и сделать их простыми.
Чтобы создать кластер на Amazon AWS, смотрите скрипт Сценарий AWS Init для обработки ГИС Jupyter Python , предназначенный для одновременного создания набора экземпляров и запуска планировщика ресурсов на основном экземпляре.
Визуализация с использованием Visualization Toolkit (VTK) и PyVista далека от Matplotlib. Идея сделать такой пример пришла ко мне давно, так как я регулярно работаю над множеством вычислительных задач, в том числе для различных университетов и геологоразведочной отрасли, и хорошо знаком с проблемами переносимости и поддержки программ, а также проблемы работы с так называемыми большими данными (сотни гигабайт и терабайт) и визуализации результатов.
Поэтому естественно возникло желание сделать пример ноутбука, который бы кратко и просто показывал красивую визуализацию и распараллеливание и ускорение кода Python, и чтобы этот ноутбук можно было запускать без изменений как локально, так и на кластере.
Все использованные библиотеки доступны уже много лет, но малоизвестны или, как говорится, остаются широко известными в узких кругах.
Оставалось только найти подходящую задачу, на которой все это можно было бы показать, и это было, пожалуй, самое сложное — ведь хотелось, чтобы пример был достаточно осмысленным и полезным.
И была найдена такая задача – рассмотреть моделирование гравитационного поля на поверхности для заданной (синтетической в данном случае) модели плотности недр и некоторых последующих преобразований с расчетом фрактального индекса из составляющих пространственный спектр и преобразование кольца Радона, как его называют математики, или Хафа, согласно компьютерным наукам.
Замечательно то, что с помощью популярных библиотек Python эти преобразования выполняются всего за несколько строк кода, что особенно ценно для примера .
Поскольку моделирование поля в каждой точке поверхности требует расчетов для всего трехмерного объема, мы будем обрабатывать гигантский объем данных.
Для визуализации используем гуманную обертку ПиВиста для библиотеки ВТК — набор инструментов визуализации , ведь написание кода для последнего - это путь настоящего джедая.
кто хочет убедиться в этом сам, смотрите мой модуль для ParaView Плагин N-Cube ParaView для визуализации 3D/4D ГИС-данных , написанный на Python + VTK. Теперь я предлагаю вам подписаться перейдите по ссылке на страницу репозитория GitHub или сразу открыть ноутбук Basic.ipynb Я надеюсь, что код довольно легко читается; Я остановлюсь лишь на некоторых особенностях.
Локальный даск-кластер, работающий на ноутбуке, предназначен для работы на многоядерных компьютерах, но для работы в кластере потребуется настроить подключение к его планировщику.
В сценарии, указанном выше Сценарий AWS Init для обработки ГИС Jupyter Python есть соответствующие комментарии и ссылки.
В коде мы используем векторизацию NumPy, то есть сразу передаем массивы, а не скаляры, и пользуемся тем, что объекты XArray предоставляют доступ к внутренним объектам NumPy (object.values).
Ускорить код NumPy непросто, но используя Numba для такого кода, вы можете получить некоторый прирост скорости (возможно, даже около 15%):
Чтобы перебрать все точки на поверхности куба и вычислить для каждой такой точки гравитационное воздействие от каждой точки куба (это, кстати, получается пятикратный интеграл), запишем:from numba import jit @jit(nopython=True, parallel=True) def delta_grav_vertical(delta_mass, x, y, z): G=6.67408*1e-11 return -np.sum((100.*1000)*G*delta_mass*z/np.power(x**2 + y**2 + z**2, 1.5))
def forward_gravity(da):
(da_y, da_x, da_z) = xr.broadcast(da.y, da.x, da.z)
deltagrav = lambda x0, y0: delta_grav_vertical(da.values.ravel(), (da_x.values.ravel()-x0), (da_y.values.ravel()-y0), (da_z.values.ravel()-0))
gravity = xr.apply_ufunc(deltagrav, da_x.isel(z=0).
chunk(50), da_y.isel(z=0).
chunk(50), vectorize=True, dask='parallelized')
.
Здесь xarray.broadcast с линеаризацией массива функцией ravel() позволяет получить тройки координат для каждой точки куба из трех одномерных координат x, y, z. Выражения da_x.isel(z=0) и da_y.isel(z=0) извлекают координаты x, y верхней поверхности куба, на которой рассчитывается гравитационное поле (точнее, его вертикальная составляющая, поскольку оно именно это измеряется в практических исследованиях и такие данные доступны для анализа).
Функция xarray.apply_ufunc() очень универсальна и одновременно обеспечивает векторизацию и поддержку параллельных отложенных вычислений для указанной функции обратного вызова deltagrad. Хитрость в том, что для выполнения вычислений на кубе для каждой точки поверхности нужно передавать координаты поверхности в виде массивов XArray, а для использования dask они также должны быть массивами dask, что мы предоставляем с помощью da_x.isel(z =0).
кусковые конструкции.
(50) и da_y.isel(z=0).
chunk(50), где 50 — размер блока по координатам x, y (выбирается в зависимости от размера массивов и количества доступных вычислительных потоков).
Да, в этом вся магия — вам просто нужно использовать вызов chunk() для массива XArray, чтобы автоматически превратить его в массив dask. Обратите внимание, что dask-расчеты по умолчанию являются ленивыми, то есть вызов функции front_gradity() завершается практически мгновенно, но возвращаемый результат — это всего лишь обертка, которая инициирует вычисления только при прямом доступе к данным или вызове load().
Это очень удобно для интерактивной работы, так как мы можем написать сложный конвейер с большими наборами данных и выбрать для тестирования и визуализации лишь небольшую его часть, а при необходимости запустить расчеты на полном наборе данных.
Например, мне часто приходится работать с наборами данных NetCDF глобальной топографии планеты и другими в сотни гигабайт на своем ноутбуке - визуализируя небольшую часть данных, а затем запускать готовый ноутбук в облаке для обработки.
все данные.
Таким образом, код для локальной работы и его производственная версия ничем не отличаются.
Главное правильно задать размеры даск-блоков, иначе все волшебство «сломается».
Код кольцевого преобразования прост и полностью основан на стандартной 2D-свертке с кольцевой маской.
Есть небольшой нюанс с расчетом фрактального индекса.
А именно, он использует полосовую фильтрацию Гаусса исходного растра для выделения спектральных составляющих (пространственных частотных диапазонов), мощность которых рассчитывается как стандартное отклонение значений результирующих матриц.
Для нашей синтетической модели все могло быть проще, но для реальных данных этот метод оказывается оптимальным из-за устойчивости к очень зашумленным данным и отсутствия проблем с границами.
Далее стандартным способом рассчитывается наклон кривой на двойном логарифмическом графике (логарифм длины волны и логарифм мощности компонента).
В заключение приглашаю всех посетить Репозитории GitHub со множеством геологических моделей и их визуализацией в Blender и ParaView, а также примерами различных анализов.
Также смотрите готовые визуализации на сайте YouTube канал .
Теги: #Популярная наука #программирование #открытый исходный код #Визуализация данных #pandas #Numpy #Геоинформационные сервисы #Dask #matplotlib #SciPy #distributed #vtk #Numba #pyvista #xarray
-
Вр – Выпуск №5
19 Oct, 24 -
Прогрессивная Шкала Налогообложения
19 Oct, 24
