<jupyter_start><jupyter_text>2D Isostatic gravity inversion - Figures Este [IPython Notebook](http://ipython.org/videos.htmlthe-ipython-notebook) utiliza a biblioteca de código aberto [Fatiando a Terra](http://fatiando.org/)<jupyter_code>Figures - magma poor margin<jupyter_output><empty_output><jupyter_text>%matplotlib inlineimport numpy as npfrom scipy.misc import derivativeimport scipy as spyfrom scipy import interpolateimport matplotlibmatplotlib.use('TkAgg', force=True)import matplotlib.pyplot as pltfrom matplotlib.ticker import MaxNLocatorimport mathimport cPickle as pickleimport datetimeimport string as stfrom scipy.misc import imreadfrom __future__ import divisionfrom fatiando import gravmag, mesher, utils, gridderfrom fatiando.mesher import Prism, Polygonfrom fatiando.gravmag import prismfrom fatiando.utils import ang2vec, si2nt, contaminatefrom fatiando.gridder import regularfrom fatiando.vis import mplfrom numpy.testing import assert_almost_equalfrom numpy.testing import assert_array_almost_equalfrom pytest import raisesplt.rc('font', size=16) import functions as fc<jupyter_code>## Observation coordinates.<jupyter_output><empty_output><jupyter_text>Model`s limitsymin = 0.0ymax = 383000.0zmin = -1000.0zmax = 45000.0xmin = -100000.0xmax = 100000.0area = [ymin, ymax, zmax, zmin] ny = 150  number of observation datas and number of prisms along the profile coordinates defining the horizontal boundaries of the adjacent columns along the profiley = np.linspace(ymin, ymax, ny)  coordinates of the center of the columns forming the interpretation modeln = ny - 1dy = (ymax - ymin)/nycmin = ymin + 0.5*dyycmax = ymax - 0.5*dyyc = np.reshape(np.linspace(ycmin, ycmax, n),(n,1))x = np.zeros_like(yc)z = np.zeros_like(yc)-150.0  Edge extension (observation coordinates)sigma = 2.0edge = sigma*dy*n<jupyter_code>## Model parameters<jupyter_output><empty_output><jupyter_text>Model densities Indices and polygons relationship: cc = continental crust layer oc = ocean crust layer w = water layer s = sediment layer m = mantle layerdw = np.array([1030.0])ds0 = np.array([2350.0])ds1 = np.array([2600.0])dcc = np.array([2870.0])doc = np.array([2885.0])dm = np.array([3300.0])dc = dcc coordinate defining the horizontal boundaries of the continent-ocean boundaryCOT = 350000.0 list defining crust density variancedc = np.zeros_like(yc)aux = yc <= COTfor i in range(len(yc[aux])):    dc[i] = dccfor i in range(len(yc[aux]),n):    dc[i] = doc     defining sediments layers density matrixds = np.vstack((np.reshape(np.repeat(ds0,n),(1,n)),np.reshape(np.repeat(ds1,n),(1,n)))) S0 => isostatic compensation surface (Airy's model)S0 = np.array([44000.0]) original<jupyter_code>## Synthetic data<jupyter_output><empty_output><jupyter_text>gsyn = np.reshape(np.loadtxt('../data/magma-poor-margin-synthetic-gravity-data.txt'),(n,1))<jupyter_code>## Water bottom<jupyter_output><empty_output><jupyter_text>bathymetry = np.reshape(np.loadtxt('../data/etopo1-pelotas.txt'),(n,1))tw = 0.0 - bathymetry<jupyter_code>## True surfaces<jupyter_output><empty_output><jupyter_text>toi = np.reshape(np.loadtxt('../data/volcanic-margin-true-toi-surface.txt'),(n,1))true_basement = np.reshape(np.loadtxt('../data/volcanic-margin-true-basement-surface.txt'),(n,1))true_moho = np.reshape(np.loadtxt('../data/volcanic-margin-true-moho-surface.txt'),(n,1)) True reference moho surface (SR = S0+dS0)true_S0 = np.array([44000.0])true_dS0 = np.array([2200.0])orignal  layer sediments thicknessts0 = toi - tw True layer sediments thicknesstrue_ts1 = true_basement - toi True thickness sediments vectortrue_ts = np.vstack((np.reshape(ts0,(1,n)),np.reshape(true_ts1,(1,n)))) True layer mantle thicknesstrue_tm = S0 - true_moho true parameters vectorptrue = np.vstack((true_ts1, true_tm, true_dS0))<jupyter_code>## Initial guess surfaces<jupyter_output><empty_output><jupyter_text>initial guess basement surfaceini_basement = np.reshape(np.loadtxt('../data/volcanic-margin-initial-basement-surface.txt'),(n,1)) initial guess moho surfaceini_moho = np.reshape(np.loadtxt('../data/volcanic-margin-initial-moho-surface.txt'),(n,1)) initial guess reference moho surface (SR = S0+dS0)ini_dS0 = np.array([11500.0]) originalini_RM = S0 + ini_dS0<jupyter_code>## Known depths<jupyter_output><empty_output><jupyter_text>Known values: basement and moho surfacesbase_known = np.loadtxt('../data/volcanic-margin-basement-known-depths.txt')base_known = np.loadtxt('../data/volcanic-margin-basement-more-known-depths.txt')base_known_new = np.loadtxt('../data/volcanic-margin-basement-new-known-depths.txt')base_known = np.loadtxt('../data/volcanic-margin-basement-few-more-known-depths.txt')base_known_new = np.loadtxt('../data/volcanic-margin-basement-few-new-known-depths.txt')base_known_old = np.loadtxt('../data/volcanic-margin-basement-known-depths.txt')moho_known = np.loadtxt('../data/volcanic-margin-moho-known-depths.txt')<jupyter_code>## Initial guess data<jupyter_output><empty_output><jupyter_text>g0 = np.reshape(np.loadtxt('../data/magma-poor-margin-initial-guess-gravity-data.txt'),(n,1))<jupyter_code>## Inversion model<jupyter_output><empty_output><jupyter_text>g = np.reshape(np.loadtxt('../data/magma-poor-margin-predicted-gravity-data-alphas_-12(1)_-8(0)_-8(0)_-7(1).txt'),(n,1))p = np.reshape(np.loadtxt('../data/magma-poor-margin-parameter-vector-alphas_-12(1)_-8(0)_-8(0)_-7(1).txt'),(2*n+1,1))gama_list = np.loadtxt('../data/magma-poor-margin-gama-list-alphas_-12(1)_-8(0)_-8(0)_-7(1).txt')<jupyter_code>p = np.reshape(np.loadtxt('../data/volcanic-margin-parameter-vector-alphas_X_-8(1)_-8(1)_-7(2)_more-known-depths.txt'),(2*n+1,1))
g = np.reshape(np.loadtxt('../data/volcanic-margin-predicted-gravity-data-alphas_X_-8(1)_-8(1)_-7(2)_more-known-depths.txt'),(n,1))
#gama_list = np.loadtxt('../data/volcanic-margin-gama-list-alphas_X_-8(1)_-8(1)_-7(2)_more-known-depths.txt')<jupyter_output><empty_output><jupyter_text>g0 = g.copy()ini_basement = tw + ts0 + p[0:n]ini_moho = S0 - p[n:n+n]ini_dS0 = p[n+n] p = np.reshape(np.loadtxt('../data/magma-poor-margin-parameter-vector-alphas_-12(1)_-8(0)_-8(0)_-7(1)-sgm_001.txt'),(2*n+1,1))g = np.reshape(np.loadtxt('../data/magma-poor-margin-predicted-gravity-data-alphas_-12(1)_-8(0)_-8(0)_-7(1)-sgm_001.txt'),(n,1))gama_list = np.loadtxt('../data/magma-poor-margin-gama-list-alphas_-12(1)_-8(0)_-8(0)_-7(1)-sgm_001.txt')<jupyter_code>p = np.reshape(np.loadtxt('../data/volcanic-margin-parameter-vector-alphas_-10(2)_-8(1)_-8(1)_-7(2)-sgm_21_more-known-depths.txt'),(2*n+1,1))
g = np.reshape(np.loadtxt('../data/volcanic-margin-predicted-gravity-data-alphas_-10(2)_-8(1)_-8(1)_-7(2)-sgm_21_more-known-depths.txt'),(n,1))
#gama_list = np.loadtxt('../data/volcanic-margin-gama-list-alphas_-10(2)_-8(1)_-8(1)_-7(2)-sgm_21_more-known-depths.txt')<jupyter_output><empty_output><jupyter_text>Inverrsion resultsRM = S0 + p[n+n]basement = tw + ts0 + p[0:n]moho = S0 - p[n:n+n] print  p[n+n]<jupyter_code>## Lithostatic Stress<jupyter_output><empty_output><jupyter_text>sgm_true = 9.81*(10**(-6))*(dw*tw + ds0*ts0 + ds1*true_ts1 + dc*(S0-tw-ts0-true_ts1-true_tm)+dm*true_tm)sgm = 9.81*(10**(-6))*(dw*tw + ds0*ts0 + ds1*p[0:n] + dc*(S0-tw-ts0-p[0:n]-p[n:n+n])+dm*p[n:n+n])<jupyter_code>## Inversion model plot<jupyter_output><empty_output><jupyter_text>polygons_water = []for (yi, twi) in zip(yc, tw):    y1 = yi - 0.5*dy    y2 = yi + 0.5*dy    polygons_water.append(Polygon(np.array([[y1, y2, y2, y1],                                             [0.0, 0.0, twi, twi]]).T,                                               props={'density': dw - dcc}))polygons_sediments0 = []for (yi, twi, s0i) in zip(yc, np.reshape(tw,(n,)), np.reshape(toi,(n,))):    y1 = yi - 0.5*dy    y2 = yi + 0.5*dy    polygons_sediments0.append(Polygon(np.array([[y1, y2, y2, y1],                                                 [twi, twi, s0i, s0i]]).T,                                                    props={'density': ds0 - dcc}))    polygons_sediments1 = []for (yi, s0i, s1i) in zip(yc, np.reshape(toi,(n,)), np.reshape(basement,(n,))):    y1 = yi - 0.5*dy    y2 = yi + 0.5*dy    polygons_sediments1.append(Polygon(np.array([[y1, y2, y2, y1],                                                 [s0i, s0i, s1i, s1i]]).T,                                                   props={'density': ds1 - dcc}))    polygons_crust = []for (yi, si, Si, dci) in zip(yc, np.reshape(basement,(n,)), np.reshape(moho,(n,)), dc):    y1 = yi - 0.5*dy    y2 = yi + 0.5*dy    polygons_crust.append(Polygon(np.array([[y1, y2, y2, y1],                                             [si, si, Si, Si]]).T,                                                props={'density': dci - dcc})) polygons_mantle = []for (yi, Si) in zip(yc, np.reshape(moho,(n,))):    y1 = yi - 0.5*dy    y2 = yi + 0.5*dy    polygons_mantle.append(Polygon(np.array([[y1, y2, y2, y1],                                                 [Si, Si, (S0+p[n+n]), (S0+p[n+n])]]).T,                                                    props={'density': dm - dcc})) %matplotlib inlineplt.close('all')fig = plt.figure(figsize=(12,16))import matplotlib.gridspec as gridspecheights = [8, 8, 8, 1]gs = gridspec.GridSpec(4, 1, height_ratios=heights)ax1 = plt.subplot(gs[0])ax2 = plt.subplot(gs[1])ax3 = plt.subplot(gs[2])ax4 = plt.subplot(gs[3])ax1.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='--', linewidth=1)ax1.plot(0.001*yc, gsyn, 'or', mfc='none', markersize=8, label='simulated data')ax1.plot(0.001*yc, g0, '-b', linewidth=2, label='initial guess data')ax1.plot(0.001*yc, g, '-g', linewidth=2, label='predicted data')ax1.set_xlim(0.001*ymin, 0.001*ymax)ax1.set_ylabel('gravity disturbance (mGal)', fontsize=16)ax1.set_xticklabels(['%g'% (l) for l in ax1.get_xticks()], fontsize=14)ax1.set_yticklabels(['%g'% (l) for l in ax1.get_yticks()], fontsize=14)ax1.legend(loc='best', fontsize=14, facecolor='silver')ax2.plot(0.001*yc, sgm_true, 'or', mfc='none', markersize=8, label='simulated lithostatic stress')ax2.plot(0.001*yc, sgm, '-g', linewidth=2, label='evaluated lithostatic stress')ax2.set_xlim(0.001*ymin, 0.001*ymax)ax2.set_ylabel('lithostatic stress (MPa)', fontsize=16)ax2.set_xticklabels(['%g'% (l) for l in ax2.get_xticks()], fontsize=14)ax2.set_yticklabels(['%g'% (l) for l in ax2.get_yticks()], fontsize=14)ax2.legend(loc='best', fontsize=14, facecolor='silver')ax3.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=1)aux = yc <= COTfor (pwi) in (polygons_water):    tmpx = [x for x in pwi.x]    tmpx.append(pwi.x[0])    tmpy = [y for y in pwi.y]    tmpy.append(pwi.y[0])    ax3.plot(tmpx, tmpy, linestyle='None')    ax3.fill(tmpx, tmpy, color='lightskyblue')for (ps0i) in (polygons_sediments0):    tmpx = [x for x in ps0i.x]    tmpx.append(ps0i.x[0])    tmpy = [y for y in ps0i.y]    tmpy.append(ps0i.y[0])    ax3.plot(tmpx, tmpy, linestyle='None')    ax3.fill(tmpx, tmpy, color='tan')for (ps1i) in (polygons_sediments1):    tmpx = [x for x in ps1i.x]    tmpx.append(ps1i.x[0])    tmpy = [y for y in ps1i.y]    tmpy.append(ps1i.y[0])    ax3.plot(tmpx, tmpy, linestyle='None')    ax3.fill(tmpx, tmpy, color='rosybrown')for (pci) in (polygons_crust[:len(yc[aux])]):    tmpx = [x for x in pci.x]    tmpx.append(pci.x[0])    tmpy = [y for y in pci.y]    tmpy.append(pci.y[0])    ax3.plot(tmpx, tmpy, linestyle='None')    ax3.fill(tmpx, tmpy, color='orange')for (pcoi) in (polygons_crust[len(yc[aux]):n]):    tmpx = [x for x in pcoi.x]    tmpx.append(pcoi.x[0])    tmpy = [y for y in pcoi.y]    tmpy.append(pcoi.y[0])    ax3.plot(tmpx, tmpy, linestyle='None')    ax3.fill(tmpx, tmpy, color='olive')for (pmi) in (polygons_mantle):    tmpx = [x for x in pmi.x]    tmpx.append(pmi.x[0])    tmpy = [y for y in pmi.y]    tmpy.append(pmi.y[0])    ax3.plot(tmpx, tmpy, linestyle='None')    ax3.fill(tmpx, tmpy, color='pink')ax3.axhline(y=S0, xmin=ymin, xmax=ymax, color='w', linestyle='--', linewidth=3)ax3.plot(yc, tw, '-k', linewidth=3)ax3.plot(yc, toi, '-k', linewidth=3)ax3.plot(yc, true_basement, '-k', linewidth=3, label='true surfaces')ax3.plot(yc, true_moho, '-k', linewidth=3)ax3.plot(yc, ini_basement, '-.b', linewidth=3, label='initial guess surfaces')ax3.plot(yc, ini_moho, '-.b', linewidth=3)ax3.plot(yc, basement, '--w', linewidth=3, label='estimated surfaces')ax3.plot(yc, moho, '--w', linewidth=3)ax3.axhline(y=true_S0+true_dS0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=3)ax3.axhline(y=S0+ini_dS0, xmin=ymin, xmax=ymax, color='b', linestyle='-.', linewidth=3)ax3.axhline(y=S0+p[n+n], xmin=ymin, xmax=ymax, color='w', linestyle='--', linewidth=3)ax3.plot(base_known[:,0], base_known[:,1], 'v', color = 'yellow', markersize=15, label='known depths (basement)')ax3.plot(base_known_old[:,0], base_known_old[:,1], 'v', color = 'yellow', markersize=15, label='known depths (basement)')ax3.plot(base_known_new[:,0], base_known_new[:,1], 'v', color = 'magenta', markersize=15, label='more known depths (basement)')ax3.plot(moho_known[:,0], moho_known[:,1], 'D', color = 'lime', markersize=15, label='known depths (moho)')ax3.set_ylim((S0+p[n+n]), zmin)ax3.set_ylim((56000.0), zmin)ax3.set_xlim(ymin, ymax)ax3.set_xlabel('y (km)', fontsize=16)ax3.set_ylabel('z (km)', fontsize=16)ax3.set_xticklabels(['%g'% (0.001*l) for l in ax3.get_xticks()], fontsize=14)ax3.set_yticklabels(['%g'% (0.001*l) for l in ax3.get_yticks()], fontsize=14)ax3.legend(loc='lower right', fontsize=14, facecolor='silver')    X, Y = fig.get_dpi()*fig.get_size_inches()plt.title('Density contrast (kg/m$^{3}$)', fontsize=17)plt.title('Density (kg/m$^{3}$)', fontsize=17)ax4.axis('off')layers_list1 = ['water', 'sediment 1', 'sediment 2', 'continental', 'oceanic', 'mantle']layers_list2 = ['', '', '', 'crust', 'crust', '']colors_list = ['lightskyblue', 'tan', 'rosybrown', 'orange', 'olive', 'pink']density_list = ['-1840', '-520', '-270', '0', '15', '430'] originaldensity_list = ['1030', '2350', '2600', '2870', '2885', '3300']ncols = len(colors_list)nrows = 1h = Y / nrowsw = X / (ncols + 1)i=ncols-1for color, density, layers1, layers2 in zip(colors_list, density_list, layers_list1, layers_list2):    col = i // nrows    row = i % nrows    x = X - (col*w) - w        yi_line = Y    yf_line = Y - Y*0.15    yi_text1 = Y - Y*0.2    yi_text2 = Y - Y*0.28    yi_text3 = Y - Y*0.08        i-=1    poly = Polygon(np.array([[x, x+w*0.75, x+w*0.75, x], [yi_line, yi_line, yf_line, yf_line]]).T)    tmpx = [x for x in poly.x]    tmpx.append(poly.x[0])    tmpy = [y for y in poly.y]    tmpy.append(poly.y[0])    ax4.plot(tmpx, tmpy, linestyle='-', color='k', linewidth=1)    ax4.fill(tmpx, tmpy, color=color)    ax4.text(x+w*0.375, yi_text1, layers1, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')    ax4.text(x+w*0.375, yi_text2, layers2, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')    ax4.text(x+w*0.375, yi_text3, density, color = 'k', fontsize=(w*0.14), horizontalalignment='center', verticalalignment='center')plt.tight_layout()        mpl.savefig('../manuscript/figures/magma-poor-margin-grafics-estimated-model-alphas_3_2_0_1.png', dpi='figure', bbox_inches='tight')mpl.savefig('../manuscript/figures/magma-poor-margin-grafics-estimated-model-alphas_1_0_0_1-sgm_001.png', dpi='figure', bbox_inches='tight')plt.show()<jupyter_code># contrastes de densidade
print 'ddw = ', dw - dcc
print 'dds0 = ', ds0 - dcc
print 'dds1 = ', ds1 - dcc
print 'ddcc = ', dcc - dcc
print 'ddoc = ', doc - dcc
print 'ddm = ', dm - dcc<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1903930695629687, 'fraction_numerical': 0.047840703387638996, 'mean_word_length': 4.763412816691505, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '576522', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6650, 'n_tokens_neox': 6309, 'n_words': 1003, 'path': 'code/magma-poor-margin-model-(figures).ipynb', 'repo_name': 'pinga-lab/paper-isostatic-grav2D', 'size': '15468'}