[Open source, Программирование, Геоинформационные сервисы, Визуализация данных, Научно-популярное] Вычислительная геология и визуализация: пример Python 3 Jupyter Notebook
Автор
Сообщение
news_bot ®
Стаж: 6 лет 9 месяцев
Сообщений: 27286
Сегодня вместо обсуждения геологических моделей мы посмотрим пример их программирования в среде Jupyter Notebook на языке Python 3 и с библиотеками Pandas, NumPy, SciPy, XArray, Dask Distributed, Numba, VTK, PyVista, Matplotlib. Это довольно простой ноутбук с поддержкой многопоточной работы и возможностью запуска локально и в кластере для обработки больших данных, отложенными вычислениями (ленивыми) и наглядной трехмерной визуализацией результатов. В самом деле, я постарался собрать разом целый набор сложных технических концепций и сделать их простыми. Для создания кластера на Amazon AWS смотрите скрипт AWS Init script for Jupyter Python GIS processing, предназначенный для единовременного создания набора инстансов и запуска планировщика ресурсов на главном инстансе.
Визуализация с помощью Visualization Toolkit(VTK) и PyVista это уже далеко не Matplotlib
Идея сделать такой пример возникла у меня давно, поскольку я регулярно занимаюсь разнообразными вычислительными задачами, в том числе для различных университетов и для геологоразведочной индустрии, и знаком очень близко с проблемами переносимости и поддерживаемости программ, а также проблемами работы с так называемыми большими данными (сотни гигабайт и терабайты) и визуализацией результатов. Так что само собой появилось желание сделать ноутбук-пример, в котором коротко и просто показать и красивую визуализацию и распараллеливание и ускорение кода Python и чтобы этот ноутбук можно было без изменений запустить как локально, так и на кластере. Все использованные библиотеки доступны уже много лет, но мало известны, или, как говорится, они остаются широко известными в узких кругах. Оставалось лишь найти подходящую задачку, на которой все это можно показать и это было, пожалуй, самым сложным — ведь мне хотелось, чтобы пример получился достаточно осмысленным и полезным. И вот такая задача нашлась — рассмотреть моделирование гравитационного поля на поверхности для заданной (синтетической в данном случае) модели плотности недр и некоторые последующие преобразования с вычислением фрактального индекса по компонентам пространственного спектра и кольцевого преобразования Радона, как его называют математики, или Хафа, согласно компьютерным наукам. Замечательно то, что с популярными библиотеками Python эти преобразования делаются буквально в несколько строчек кода, что особенно ценно для примера. Поскольку моделирование поля в каждой точке поверхности требует вычисления для всего трехмерного объема, мы будем обрабатывать гигантский объем данных. Для визуализации используем человеколюбивую обертку PyVista для библиотеки VTK — Visualization Toolkit, потому что писать код для последней это путь истинных джедаев… кто хочет лично в том убедиться, смотрите мой модуль к ParaView N-Cube ParaView plugin for 3D/4D GIS Data Visualization, написанный как раз на Python + VTK.
Теперь предлагаю проследовать по ссылке на страницу GitHub репозитория или сразу открыть ноутбук basic.ipynb Надеюсь, код достаточно просто читается, остановлюсь лишь на нескольких особенностях. Запускаемый в ноутбуке локальный кластер dask предназначен для работы на многоядерных компьютерах, а вот для работы в кластере потребуется настроить подключение к его планировщику. В упомянутом выше скрипте AWS Init script for Jupyter Python GIS processing есть соответствующие комментарии и ссылки. В коде мы используем векторизацию 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() весьма универсальная и одновременно обеспечивает векторизацию и поддержку параллельных ленивых вычислений dask для указанной коллбэк функции deltagrav. Хитрость заключается в том, что для выполнения вычислений на кубе для каждой точки поверхности нужно координаты поверхности передать в виде XArray массивов, а для использования dask они также должны быть dask массивами, что мы и обеспечиваем конструкциями da_x.isel(z=0).chunk(50) и da_y.isel(z=0).chunk(50), где 50 это размер блока по координатам x, y (подбирается в зависимости от размера массивов и количества доступных вычислительных потоков). Да, такая вот магия — достаточно лишь использовать вызов chunk() для XArray массива, чтобы автоматически превратить его в dask массив.
Обратим внимание, что dask-вычисления по умолчанию являются ленивыми (отложенными), то есть вызов функции forward_gravity() завершается почти мгновенно, но возвращаемый результат является лишь оберткой, которая инициирует вычисления только при непосредственном обращении к данным или вызовом load(). При интерактивной работе это очень удобно, так как мы можем написать сложный пайплайн с большими наборами данных и для проверки и визуализации выбрать лишь маленький его кусочек, а при необходимости и запустить вычисления на полном наборе данных. К примеру, мне часто приходится работать с NetCDF датасетами глобального рельефа планеты и прочими в сотни гигабайт на своем ноутбуке — визуализируя малую часть данных, а потом запускать уже готовый ноутбук в облаке для обработки всех данных. Таким образом, код для локальной работы и его продакшен версия ничем не отличаются. Главное, правильно настроить размеры dask блоков, иначе вся магия "сломается".
Код кольцевого преобразования простой и полностью основан на стандартной двумерной свертке с кольцевой маской. Вот с вычислением фрактального индекса есть небольшой нюанс. А именно, здесь используется гауссова полосовая фильтрация исходного растра для выделения спектральных компонентов (диапазонов пространственных частот), мощность которых вычисляется как стандартное отклонение значений полученных матриц. Для нашей синтетической модели можно бы и попроще, а вот для реальных данных этот метод оказывается оптимальным в силу своей устойчивости к очень зашумленным данным и отсутствия проблем с границами. Далее стандартным способом вычисляется наклон кривой на двойном логарифмическом графике (логарифм длины волны и логарифм мощности компоненты).
В заключение, приглашаю всех посетить GitHub репозитории с множеством геологических моделей и их визуализацией в Blender и ParaView, а также примерами различного анализа. Также смотрите готовые визуализации на YouTube канале.
===========
Источник:
habr.com
===========
Похожие новости:
- [JavaScript, Программирование, Совершенный код, ReactJS] Удвойте скорость написания кода на React с помощью этих простых трюков (перевод)
- [Научно-популярное] Теории заговора и научные гипотезы
- [Научно-популярное, Искусственный интеллект] Может ли ИИ заменить голоса актеров из «Симпсонов» (перевод)
- [Анализ и проектирование систем, Визуализация данных, Visual Basic for Applications, Подготовка технической документации] Хорошие BPM — инструменты, которых нет и нет. Моделирование процессов
- [Научно-популярное] Фотовольтаика
- [Высокая производительность, Программирование, Алгоритмы, WebAssembly] Разгоняем JS-парсер с помощью WebAssembly (часть 3: SIMD)
- [Научно-популярное] Умные окна заменят кондиционеры и лишат работы промышленных альпинистов
- [Open source, Разработка под Linux, Софт, IT-компании] Разработчик 7-Zip выпустил официальный билд для Linux спустя 22 года после выхода Windows-версии
- [C++, C, Программирование микроконтроллеров, Компьютерное железо] Размещение кучи FreeRTOS в разделе CCMRAM для STM32
- [Карьера в IT-индустрии, Научно-популярное, Мозг, Здоровье] Выгорание — норма? Часть 2. Подробно про фазы и «Что делать?»
Теги для поиска: #_open_source, #_programmirovanie (Программирование), #_geoinformatsionnye_servisy (Геоинформационные сервисы), #_vizualizatsija_dannyh (Визуализация данных), #_nauchnopopuljarnoe (Научно-популярное), #_pandas, #_numpy, #_scipy, #_xarray, #_dask, #_distributed, #_numba, #_vtk, #_pyvista, #_matplotlib, #_open_source, #_programmirovanie (
Программирование
), #_geoinformatsionnye_servisy (
Геоинформационные сервисы
), #_vizualizatsija_dannyh (
Визуализация данных
), #_nauchnopopuljarnoe (
Научно-популярное
)
Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы
Текущее время: 22-Ноя 14:52
Часовой пояс: UTC + 5
| Автор | Сообщение |
|---|---|
|
news_bot ®
Стаж: 6 лет 9 месяцев |
|
|
Сегодня вместо обсуждения геологических моделей мы посмотрим пример их программирования в среде Jupyter Notebook на языке Python 3 и с библиотеками Pandas, NumPy, SciPy, XArray, Dask Distributed, Numba, VTK, PyVista, Matplotlib. Это довольно простой ноутбук с поддержкой многопоточной работы и возможностью запуска локально и в кластере для обработки больших данных, отложенными вычислениями (ленивыми) и наглядной трехмерной визуализацией результатов. В самом деле, я постарался собрать разом целый набор сложных технических концепций и сделать их простыми. Для создания кластера на Amazon AWS смотрите скрипт AWS Init script for Jupyter Python GIS processing, предназначенный для единовременного создания набора инстансов и запуска планировщика ресурсов на главном инстансе.  Визуализация с помощью Visualization Toolkit(VTK) и PyVista это уже далеко не Matplotlib Идея сделать такой пример возникла у меня давно, поскольку я регулярно занимаюсь разнообразными вычислительными задачами, в том числе для различных университетов и для геологоразведочной индустрии, и знаком очень близко с проблемами переносимости и поддерживаемости программ, а также проблемами работы с так называемыми большими данными (сотни гигабайт и терабайты) и визуализацией результатов. Так что само собой появилось желание сделать ноутбук-пример, в котором коротко и просто показать и красивую визуализацию и распараллеливание и ускорение кода Python и чтобы этот ноутбук можно было без изменений запустить как локально, так и на кластере. Все использованные библиотеки доступны уже много лет, но мало известны, или, как говорится, они остаются широко известными в узких кругах. Оставалось лишь найти подходящую задачку, на которой все это можно показать и это было, пожалуй, самым сложным — ведь мне хотелось, чтобы пример получился достаточно осмысленным и полезным. И вот такая задача нашлась — рассмотреть моделирование гравитационного поля на поверхности для заданной (синтетической в данном случае) модели плотности недр и некоторые последующие преобразования с вычислением фрактального индекса по компонентам пространственного спектра и кольцевого преобразования Радона, как его называют математики, или Хафа, согласно компьютерным наукам. Замечательно то, что с популярными библиотеками Python эти преобразования делаются буквально в несколько строчек кода, что особенно ценно для примера. Поскольку моделирование поля в каждой точке поверхности требует вычисления для всего трехмерного объема, мы будем обрабатывать гигантский объем данных. Для визуализации используем человеколюбивую обертку PyVista для библиотеки VTK — Visualization Toolkit, потому что писать код для последней это путь истинных джедаев… кто хочет лично в том убедиться, смотрите мой модуль к ParaView N-Cube ParaView plugin for 3D/4D GIS Data Visualization, написанный как раз на Python + VTK. Теперь предлагаю проследовать по ссылке на страницу GitHub репозитория или сразу открыть ноутбук basic.ipynb Надеюсь, код достаточно просто читается, остановлюсь лишь на нескольких особенностях. Запускаемый в ноутбуке локальный кластер dask предназначен для работы на многоядерных компьютерах, а вот для работы в кластере потребуется настроить подключение к его планировщику. В упомянутом выше скрипте AWS Init script for Jupyter Python GIS processing есть соответствующие комментарии и ссылки. В коде мы используем векторизацию 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() весьма универсальная и одновременно обеспечивает векторизацию и поддержку параллельных ленивых вычислений dask для указанной коллбэк функции deltagrav. Хитрость заключается в том, что для выполнения вычислений на кубе для каждой точки поверхности нужно координаты поверхности передать в виде XArray массивов, а для использования dask они также должны быть dask массивами, что мы и обеспечиваем конструкциями da_x.isel(z=0).chunk(50) и da_y.isel(z=0).chunk(50), где 50 это размер блока по координатам x, y (подбирается в зависимости от размера массивов и количества доступных вычислительных потоков). Да, такая вот магия — достаточно лишь использовать вызов chunk() для XArray массива, чтобы автоматически превратить его в dask массив. Обратим внимание, что dask-вычисления по умолчанию являются ленивыми (отложенными), то есть вызов функции forward_gravity() завершается почти мгновенно, но возвращаемый результат является лишь оберткой, которая инициирует вычисления только при непосредственном обращении к данным или вызовом load(). При интерактивной работе это очень удобно, так как мы можем написать сложный пайплайн с большими наборами данных и для проверки и визуализации выбрать лишь маленький его кусочек, а при необходимости и запустить вычисления на полном наборе данных. К примеру, мне часто приходится работать с NetCDF датасетами глобального рельефа планеты и прочими в сотни гигабайт на своем ноутбуке — визуализируя малую часть данных, а потом запускать уже готовый ноутбук в облаке для обработки всех данных. Таким образом, код для локальной работы и его продакшен версия ничем не отличаются. Главное, правильно настроить размеры dask блоков, иначе вся магия "сломается". Код кольцевого преобразования простой и полностью основан на стандартной двумерной свертке с кольцевой маской. Вот с вычислением фрактального индекса есть небольшой нюанс. А именно, здесь используется гауссова полосовая фильтрация исходного растра для выделения спектральных компонентов (диапазонов пространственных частот), мощность которых вычисляется как стандартное отклонение значений полученных матриц. Для нашей синтетической модели можно бы и попроще, а вот для реальных данных этот метод оказывается оптимальным в силу своей устойчивости к очень зашумленным данным и отсутствия проблем с границами. Далее стандартным способом вычисляется наклон кривой на двойном логарифмическом графике (логарифм длины волны и логарифм мощности компоненты). В заключение, приглашаю всех посетить GitHub репозитории с множеством геологических моделей и их визуализацией в Blender и ParaView, а также примерами различного анализа. Также смотрите готовые визуализации на YouTube канале. =========== Источник: habr.com =========== Похожие новости:
Программирование ), #_geoinformatsionnye_servisy ( Геоинформационные сервисы ), #_vizualizatsija_dannyh ( Визуализация данных ), #_nauchnopopuljarnoe ( Научно-популярное ) |
|
Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы
Текущее время: 22-Ноя 14:52
Часовой пояс: UTC + 5