آموزش مدل سازی زمین شناسی ساختاری 3D در پایتون با Gempy
ماژول Gempy به عنوان کتاب منبع باز پایتون برای تولید مدل های زمین شناسی ساختاری کامل 3D است. این کتابخانه کامل برای ایجاد مدل های زمین شناختی از interfaces، گسل ها و جهت گیری های لایه است، همچنین توالی لایه های زمین شناسی را برای نشان دادن نفوذ و خطاهای سنگ انجام می دهد.
الگوریتم برای مدل سازی زمین شناسی مبتنی بر واسنجی cokriging جهانی با حمایت از کتابخانه های ریاضی Python تحت عناوین Numpy، PyMC3 و Theano است.
Gempy یک مدل شبکه ایجاد می کند که می تواند به عنوان بخش های 2D با Matplotlib و یا به عنوان اشیاء هندسی 3D به عنوان اشیاء VTK اجازه نمایندگی از مدل های زمین شناسی در Paraview برای برش، فیلتر، شفافیت و یک ظاهر طراحی مجدد را می دهند.
این آموزش یک نمونه اساسی از راه اندازی ساختار زمین شناسی طبقه بندی شده با 5 لایه و یک خطا است. به منظور ایجاد آموزش به طور کامل در دسترس اکثر کاربران، ما یک آموزش مکمل در مورد نحوه نصب Gempy در ویندوز با یک توزیع مخزن Anaconda ایجاد کرده ایم.
این آموزش دو بخش دارد:
ویدیو قسمت اول: نصب Gempy در ویندوز - زبان اصلی یوتیوب
(عدم دسترسی از آی پی ایران)
ویدیو قسمت دوم: مثال مدل سازی زمین شناسی با Gempy - زبان اصلی یوتیوب
(عدم دسترسی از آی پی ایران)
کد پایتون
این اسکریپت برای آموزش است:
راه اندازی محیط پایتون
در این بخش، کتابخانه های مورد نیاز را برای آموزش وارد می کنیم. به اسکریپت Gempy به Numpy و Matplotlib نیاز دارد. ما یک گزینه Jupyter را برای نمایش تعاملی گرافیک Matplotlib پس از سلول های اسکریپت (٪ matplotlib inline) پیکربندی می کنیم.
توجه داشته باشید که هشدارها فقط پیام هایی هستند که کاربر هنگام اجرای اسکریپت باید در نظر بگیرد، آنها به معنی شکست یک کد نیستند. از آنجا که این آموزش بر روی ویندوز برخی از کتابخانه های تکمیلی قادر به نصب نیست، اما عملکرد کلی از کد مدل سازی زمین شناسی کامل است.
# Embedding matplotlib figures in the notebooks
%matplotlib inline
# Importing GemPy
import gempy as gp
# Importing auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
ایجاد شیء مدل زمین شناسی و تعریف استاتیک
آموزش ایجاد یک شبکه از 100 * 100 ردیف x 100 لایه در طول گسترش 2km x 2km x 2km. مقادیر بالاتر هم امکان پذیر است اما زمان محاسبه بالاتر خواهد بود. سیستم مختصات از هماهنگی های محلی استفاده می کند، آموزش پیشرفته Gempy با مختصات UTM را ارزیابی می کید.
جهت گیری و تماس زمین شناسی از فایل های CSV وارد شده و تبدیل به یک فریم اطلاعات Pandas می شود. سپس مجموعه زمین شناسی (گسل ها / سازه ها) تعریف می شود و همچنین توالی تشریح زمین شناسی تعریف می شود.
مهم آن است که اشاره شود گسل ها باید مستقل قرار گیرند، جایی که جدید ترین ورودی است.
# Importing the data from CSV-files and setting extent and resolution
geo_data = gp.create_data([0,2000,0,2000,0,2000],[100,100,100],
path_o = "../Txt/simple_fault_model_orientations.csv", # importing orientation (foliation) data
path_i = "../Txt/simple_fault_model_points.csv") # importing point-positional interface data
gp.get_data(geo_data).loc[:,['X','Y','Z','formation']].head()
# Assigning series to formations as well as their order (timewise)
gp.set_series(geo_data, {"Fault_Series":'Main_Fault',
"Strat_Series": ('Sandstone_2','Siltstone', 'Shale', 'Sandstone_1')},
order_series = ["Fault_Series", 'Strat_Series'],
order_formations=['Main_Fault',
'Sandstone_2','Siltstone', 'Shale', 'Sandstone_1',
], verbose=0)
طرح توالی زمین شناسی
Gempy دارای برخی از ویژگی های مفید برای نشان دادن سری های زمین شناسی تعریف شده و توالی های شکل گیری است.
gp.get_sequential_pile(geo_data)
<gempy.plotting.sequential_pile.StratigraphicPile at 0x107149e8>
بررسی داده های ورودی
مجموعه داده های مختلف برای ساخت مدل زمین شناسی می تواند وجود داشت باشد. در هر صورت از ویژگی های ".get_" Gempy در دسترس باشد.
# Review of the centroid coordinates from the model grid
gp.get_grid(geo_data).values
array([[ 10., 10., 10.],
[ 10., 10., 30.],
[ 10., 10., 50.],
...,
[ 1990., 1990., 1950.],
[ 1990., 1990., 1970.],
[ 1990., 1990., 1990.]], dtype=float32)
# Defined interfases from the input CSV data
gp.get_data(geo_data, 'interfaces').loc[:,['X','Y','Z','formation']].head()
# Defined layer orientations from the input CSV data
gp.get_data(geo_data, 'orientations').loc[:,['X','Y','Z','formation','azimuth']]
نمایش گرافیکی داده های ورودی
در این قسمت، بازنمودهای 2D و 3D برای ارائه بینش مشخص و جهت ها انجام شد.
gp.plot_data(geo_data, direction='y')
E:\Software\Anaconda3\lib\site-packages\gempy\gempy_front.py:927: FutureWarning: gempy plotting functionality will be moved in version 1.2, use gempy.plotting module instead
warnings.warn("gempy plotting functionality will be moved in version 1.2, use gempy.plotting module instead", FutureWarning)
gp.plotting.plot_data_3D(geo_data)
تعامل زمین شناسی
هنگامی که داده های ورودی آماده می شوند، می توانیم داده ها و پارامترهای تعبیه را با روش InterpolatonData از کتابخانه Gempy تعریف کنیم.
مدل زمین شناسی تحت روش compute_model محاسبه شده است. نتایج فرایند مدل، سنگ شناسی و گسل ها در ابعاد آرایه همانند geo_data است.
interp_data = gp.InterpolatorData(geo_data, u_grade=[1,1], output='geology', compile_theano=True, theano_optimizer='fast_compile')
Compiling theano function...
Compilation Done!
Level of Optimization: fast_compile
Device: cpu
Precision: float32
Number of faults: 1
interp_data.geo_data_res.formations.as_matrix
<bound method NDFrame.as_matrix of value formation_number
Main_Fault 1 1
Sandstone_2 2 2
Siltstone 3 3
Shale 4 4
Sandstone_1 5 5
basement 6 6>
interp_data.geo_data_res.get_formations()
[Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]
Categories (5, object): [Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]
lith_block, fault_block = gp.compute_model(interp_data)
اکتشاف مدل لیتوگرافی
بلوک سنگ شناسی دارای دو بخش است، اول اطلاعات مربوط به شکل گیری سنگ شناسی، در حالی که دوم نشانگر جهت است. در این بخش، توزیع سنگ شناسی و اطلاعات جدا شده با گسل به صورت هیستوگرام نمایش داده می شود.
lith_block[0]
array([ 6.3131361 , 6.24877167, 6.19397354, ..., 2.00398016,
2.00626612, 2.00983 ], dtype=float32)
plt.hist(lith_block[0],bins=100)
plt.show()
plt.hist(fault_block[0],bins=10)
plt.show()
نمایندگی مدل زمین شناسی
به عنوان یک آرایه دیگر، بلوک های سنگی حاصل می توانند در Matplotlib نشان داده شوند. با این حال، Gempy روش های ویژه ای برای نمایندگی مقطعی دارد. با استفاده از ویدجت Jupyter نمایشی تعاملی از بخش های زمین شناسی در امتداد جهت Y با استفاده از handelbar برای تغییر در ردیف انجام می شود.
gp.plotting.plot_section(geo_data, lith_block[0], cell_number=50, direction='y', plot_data=False)
import ipywidgets as widgets
def plotCrossSection(cell):
gp.plotting.plot_section(geo_data, lith_block[0], cell_number=cell, direction='y', plot_data=False)
widgets.interact(plotCrossSection, cell=widgets.IntSlider(min=0,max=99,step=1,value=50) )
gp.plotting.plot_scalar_field(geo_data, lith_block[1], cell_number=50, N=6,
direction='y', plot_data=False)
plt.colorbar()
plt.show()
ver_s, sim_s = gp.get_surfaces(interp_data,lith_block[1],
fault_block[1],
original_scale=True)
gp.plotting.plot_surfaces_3D(geo_data, ver_s, sim_s)
فایل های ورودی
شما می توانید فایل های ورودی این آموزش را در این لینک دانلود کنید.
شناسه تلگرام مدیر سایت: SubBasin@
نشانی ایمیل: behzadsarhadi@gmail.com
(سوالات تخصصی را در گروه تلگرام ارسال کنید)
_______________________________________________________
سلام من از بعد مرحلهی افزودن مقادیر ارتفاعی dem، با ارور اجرا نشدن زون 39 روبهرو میشم . و ارور به این شکل هست که من رو به سایت Spatial Reference ارجا میده .چطور میتونم این مشکل رو رفع کنم؟