<jupyter_start><jupyter_text>2D Isostatic gravity inversion - Initial Guess Model 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>Initial Guess model of Pelotas Basin.<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 pltimport 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.0zmax = 40000.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([2855.0])dcc = np.array([2870.0])doc = np.array([2885.0])dm = np.array([3240.0])ds0 = np.array([2425.0])T07ds1 = np.array([2835.0])T07ds0 = np.array([2425.0])T06ds1 = np.array([2750.0])T06ds0 = np.array([2450.0])T05ds1 = np.array([2700.0])T05dc = 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) SR = S0+dS0 => reference Moho (Forward modeling)S0 = np.array([41000.0])dS0 = np.array([1000.0])<jupyter_code>## Observed and synthetic data<jupyter_output><empty_output><jupyter_text>gsyn = np.loadtxt('../data/pelotas-profile-synthetic-gravity-data.txt') gobs = np.reshape(np.loadtxt('../data/pelotas-profile-gz.txt'),(n,1))<jupyter_code>## For display<jupyter_output><empty_output><jupyter_text>img = imread('../../../../Desktop/L130.jpeg') shape=(719,1238)extent = [ymin, ymax, zmax, 0.0]def calculate_aspect(shape,extent):    dx = (abs(extent[1] - extent[0]))/float(shape[1])    dy = (abs(extent[3] - extent[2]))/float(shape[0])    return dx/dy<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>## Interpreted surfaces<jupyter_output><empty_output><jupyter_text>toi = np.reshape(np.loadtxt('../data/pelotas-profile-interpreted-toi-surface.txt'),(n,1))interpreted_basement = np.reshape(np.loadtxt('../data/pelotas-profile-interpreted-basement-surface.txt'),(n,1))interpreted_moho = np.reshape(np.loadtxt('../data/pelotas-profile-interpreted-moho-surface.txt'),(n,1)) reference moho surface (SR = S0+dS0)true_dS0 = np.array([2200.0])<jupyter_code># Initial guess surfaces
### Basement surface<jupyter_output><empty_output><jupyter_text>mpl.close('all')plt.figure(figsize=(12,10))mpl.subplot(2,1,1)mpl.title('Gravity disturbance data', fontsize=14)mpl.paths([[ymin, 0.]], [[ymax, 0.]], style='--k', linewidth=1)mpl.plot(0.001*yc, gobs, label='obs')mpl.xlim(0.001*ymin, 0.001*ymax)mpl.ylabel('gravity disturbance (mGal)', fontsize=16)mpl.xticks(fontsize=12)mpl.yticks(fontsize=12)mpl.legend(loc='best')axes = mpl.subplot(1,1,1)mpl.ylim(zmax, zmin)mpl.xlim(ymin, ymax)mpl.xticks(fontsize=12)mpl.yticks(fontsize=12)mpl.xlabel('y (m)')mpl.ylabel('z (m)')mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='-k', linewidth=1)mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=1.0)mpl.plot(yc, tw, '-b', linewidth=1)mpl.plot(yc, toi, '-b', linewidth=1)mpl.plot(yc, interpreted_basement, '--b', linewidth=1)mpl.m2km()basement_picks = mpl.draw_polygon(area, axes, marker='.', color='r') basement_picks<jupyter_code>basement_picks = np.array([[   1422.64201445,    4191.22257053],
       [  10332.51817002,    3247.05179878],
       [  18393.83469173,    2817.88326616],
       [  23485.19249491,    3161.21809225],
       [  26030.87139651,    3847.88774444],
       [  29000.83011503,    5049.55963577],
       [  35789.30718594,    8139.57307061],
       [  45972.02279231,   13032.09434244],
       [  53609.05949709,   16465.44260337],
       [  60397.536568  ,   18954.62009255],
       [  66761.73382198,   20156.29198388],
       [  75671.60997755,   20499.62680997],
       [  82460.08704847,   19984.62457083],
       [  88400.00448552,   18868.78638603],
       [  91794.24302097,   18182.11673384],
       [  93067.08247177,   18010.44932079],
       [  97309.88064109,   18439.61785341],
       [ 104522.63752893,   19727.12345126],
       [ 107492.59624746,   20070.45827735],
       [ 110886.83478291,   19898.79086431],
       [ 114281.07331837,   19040.45379907],
       [ 120645.27057235,   19641.28974474],
       [ 124463.78892474,   19812.95715778],
       [ 132100.82562951,   18439.61785341],
       [ 138040.74306656,   18868.78638603],
       [ 140586.42196816,   18868.78638603],
       [ 142707.82105282,   18353.78414689],
       [ 147374.89903907,   15607.10553814],
       [ 151193.41739146,   14662.93476638],
       [ 156284.77519464,   14577.10105986],
       [ 164770.37153328,   15177.93700552],
       [ 166467.49080101,   15263.77071205],
       [ 172831.68805499,   14662.93476638],
       [ 187257.20183068,   15607.10553814],
       [ 193197.11926773,   15263.77071205],
       [ 209744.03212808,   17066.27854904],
       [ 215259.66974819,   17237.94596208],
       [ 219502.46791751,   16894.61113599],
       [ 222472.42663604,   16293.77519033],
       [ 229260.90370695,   16465.44260337],
       [ 237322.22022866,   15435.43812509],
       [ 245383.53675037,   16207.9414838 ],
       [ 251323.45418742,   15607.10553814],
       [ 265748.96796311,   15950.44036423],
       [ 275507.40375254,   14062.09882072],
       [ 283144.44045732,   15006.26959248],
       [ 289508.6377113 ,   15092.103299  ],
       [ 295872.83496528,   16722.94372294],
       [ 300115.6331346 ,   17237.94596208],
       [ 304358.43130392,   17237.94596208],
       [ 308601.22947324,   16808.77742947],
       [ 312419.74782563,   15692.93924466],
       [ 319632.50471348,   12173.75727721],
       [ 324299.58269973,   11572.92133154],
       [ 328118.10105212,   11401.2539185 ],
       [ 354423.4497019 ,   11229.58650545],
       [ 366303.284576  ,    9512.91237498],
       [ 379455.95890089,    9512.91237498]])
# change the coordinates of the extremum points in order to
# avoid problems for constructing the interpolator
basement_picks[0,0] = ymin
basement_picks[-1,0] = ymax
basement = fc.surface_interpolate_function(basement_picks,yc)
for i in range(len(basement)):
    if basement[i] < toi[i]:
        basement[i] = toi[i]
# layer sediments thickness
ts0 = toi - tw
# layer igneous thickness
ts1 = basement - toi<jupyter_output><empty_output><jupyter_text>np.savetxt('../data/pelotas-profile-initial-basement-surface.txt', basement, fmt='%.18f')<jupyter_code># thickness sediments vector
ts = np.vstack((np.reshape(ts0,(1,n)),np.reshape(ts1,(1,n))))<jupyter_output><empty_output><jupyter_text>Moho surface<jupyter_code>mpl.close('all')

#mpl.subplot(2,1,1)
#mpl.title('Gravity disturbance data', fontsize=14)
#mpl.paths([[ymin, 0.]], [[ymax, 0.]], style='--k', linewidth=1)
#mpl.plot(0.001*yc, gobs, label='obs')
#mpl.xlim(0.001*ymin, 0.001*ymax)
#mpl.ylabel('gravity disturbance (mGal)', fontsize=16)
#mpl.xticks(fontsize=12)
#mpl.yticks(fontsize=12)
#mpl.legend(loc='best')

axes = mpl.subplot(1,1,1)
mpl.ylim(zmax, zmin)
mpl.xlim(ymin, ymax)
mpl.xticks(fontsize=12)
mpl.yticks(fontsize=12)
mpl.xlabel('y (m)')
mpl.ylabel('z (m)')
mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='-k', linewidth=1)
mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=1.0)
mpl.plot(yc, tw, '-b', linewidth=1)
mpl.plot(yc, toi, '-r', linewidth=1)
mpl.plot(yc, basement, '-b', linewidth=1)
#mpl.plot(yc, interpreted_basement, '--r', linewidth=1)
mpl.m2km()

moho_picks = mpl.draw_polygon(area, axes, color='r')
moho_picks<jupyter_output><empty_output><jupyter_text>moho_picks = np.array([[   3458.72396201,   31550.90311987],       [  14945.40824487,   29485.29631288],       [  40045.94056668,   28720.25675474],       [  63444.74188363,   29638.30422451],       [  74080.56066406,   26884.1618152 ],       [  93225.03446883,   27343.18555008],       [ 102159.12224439,   25201.07478728],       [ 111944.07552239,   24130.01940588],       [ 117900.13403943,   24665.54709658],       [ 130663.11657594,   27419.6895059 ],       [ 141724.36810759,   28184.72906404],       [ 150658.45588315,   24053.51545007],       [ 158741.67815628,   19998.80579191],       [ 170228.36243914,   20763.84535005],       [ 176609.8537074 ,   19463.27810121],       [ 184267.64322931,   19922.30183609],       [ 195328.89476095,   21758.39677564],       [ 206390.1462926 ,   21146.36512912],       [ 217025.96507303,   21681.89281982],       [ 227661.78385346,   20610.83743842],       [ 234043.27512172,   20534.33348261],       [ 242126.49739484,   20304.82161517],       [ 249784.28691675,   19004.25436632],       [ 257016.64368745,   19080.75832214],       [ 269779.62622396,   20840.34930587],       [ 275310.25198978,   20457.8295268 ],       [ 280840.87775561,   20610.83743842],       [ 291476.69653604,   21911.40468727],       [ 297007.32230186,   20534.33348261],       [ 311472.03584324,   20610.83743842],       [ 321682.42187246,   18545.23063144],       [ 329765.64414558,   17321.16733841],       [ 334445.40440897,   17091.65547097],       [ 340401.46292601,   18239.21480818],       [ 345506.65594062,   17780.19107329],       [ 354440.74371618,   18086.20689655],       [ 360822.23498444,   18621.73458725],       [ 368480.02450635,   16479.62382445],       [ 381668.43979408,   15944.09613375]])  change the coordinates of the extremum points in order to avoid problems for constructing the interpolatormoho_picks[0,0] = yminmoho_picks[-1,0] = ymax moho = fc.surface_interpolate_function(moho_picks,yc)for i in range(len(moho)):    if moho[i] < basement[i]:        moho[i] = basement[i]  layer mantle thicknesstm = S0 - moho layer crust thicknesstoc = moho - tw - ts0 - ts1<jupyter_code>np.savetxt('../data/pelotas-profile-initial-moho-surface.txt', moho, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>Known depths<jupyter_code>mpl.close('all')

axes = mpl.subplot(1,1,1)
mpl.ylim(zmax, zmin)
mpl.xlim(ymin, ymax)
mpl.xticks(fontsize=12)
mpl.yticks(fontsize=12)
mpl.xlabel('y (m)')
mpl.ylabel('z (m)')
mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='-k', linewidth=1)
#mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=1.0)
mpl.plot(base_known[:,0], base_known[:,1], 'vc', markersize=9.5, label='known depths (basement)')
mpl.plot(yc, tw, '-b', linewidth=1)
mpl.plot(yc, toi, '-b', linewidth=1)
mpl.plot(yc, interpreted_basement, '-b', linewidth=1)
mpl.plot(yc, basement, '-r', linewidth=1)
mpl.m2km()

base_picks = mpl.draw_polygon(area, axes, color='r')
base_picks
for i in range(len(base_picks)):
    j = int((base_picks[i][0] - ymin)/dy)
    print '[', np.reshape(yc,(n,))[j], ',' , np.reshape(basement,(n,))[j], '],'
print '[', np.reshape(yc,(n,))[3], ',' , np.reshape(basement,(n,))[3], '],'
print '[', np.reshape(yc,(n,))[9], ',' , np.reshape(basement,(n,))[9], '],'
#print '[', np.reshape(yc,(n,))[20], ',' , np.reshape(basement,(n,))[20], '],'
print '[', np.reshape(yc,(n,))[42], ',' , np.reshape(basement,(n,))[42], '],'
#print '[', np.reshape(yc,(n,))[59], ',' , np.reshape(basement,(n,))[59], '],'
print '[', np.reshape(yc,(n,))[89], ',' , np.reshape(basement,(n,))[89], '],'
print '[', np.reshape(yc,(n,))[112], ',' , np.reshape(basement,(n,))[112], '],'
print '[', np.reshape(yc,(n,))[128], ',' , np.reshape(basement,(n,))[128], ']'<jupyter_output><empty_output><jupyter_text>Basement surface known valuesbase_known = np.array([[ 8996.6442953 , 3369.12204792 ],[ 24419.4630872 , 3413.2275905 ],[ 109244.966443 , 19981.8301933 ],[ 230057.04698 , 16363.7183832 ],[ 289177.852349 , 15087.6420083 ],[ 330305.369128 , 11386.9799145 ]])np.savetxt('../data/pelotas-profile-basement-known-depths.txt', base_known, fmt='%.18f')  Basement surface known valuesbase_known_more = np.array([[ 8996.6442953 , 3369.12204792 ],[ 24419.4630872 , 3413.2275905 ],[ 52694.6308725 , 16054.3470372 ],[ 109244.966443 , 19981.8301933 ],[ 152942.95302 , 14633.4398588 ],[ 230057.04698 , 16363.7183832 ],[ 289177.852349 , 15087.6420083 ],[ 330305.369128 , 11386.9799145 ]])np.savetxt('../data/pelotas-profile-basement-more-known-depths.txt', base_known_more, fmt='%.18f')base_known_new = np.array([[ 52694.6308725 , 16054.3470372 ],[ 152942.95302 , 14633.4398588 ]])np.savetxt('../data/pelotas-profile-basement-new-known-depths.txt', base_known_new, fmt='%.18f')<jupyter_code>mpl.close('all')

axes = mpl.subplot(1,1,1)
mpl.ylim(zmax, zmin)
mpl.xlim(ymin, ymax)
mpl.xticks(fontsize=12)
mpl.yticks(fontsize=12)
mpl.xlabel('y (m)')
mpl.ylabel('z (m)')
mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='-k', linewidth=1)
mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=1.0)
mpl.plot(yc, tw, '-b', linewidth=1)
mpl.plot(yc, toi, '-b', linewidth=1)
mpl.plot(yc, interpreted_moho, '-b', linewidth=1)
mpl.plot(yc, moho, '-r', linewidth=1)
mpl.m2km()

moho_picks = mpl.draw_polygon(area, axes, color='r')
moho_picks
for i in range(len(moho_picks)):
    j = int((moho_picks[i][0] - ymin)/dy)
    print '[', np.reshape(yc,(n,))[j], ',' , np.reshape(moho,(n,))[j], '],'
print '[', np.reshape(yc,(n,))[60], ',' , np.reshape(moho,(n,))[60], '],'
print '[', np.reshape(yc,(n,))[85], ',' , np.reshape(moho,(n,))[85], '],'
print '[', np.reshape(yc,(n,))[133], ',' , np.reshape(moho,(n,))[133], ']'<jupyter_output><empty_output><jupyter_text>Moho surface known valuesmoho_known = np.array([[ 155513.422819 , 21618.1647126 ],[ 219775.167785 , 21405.0407579 ],[ 343157.718121 , 17991.3913566 ]])np.savetxt('../data/pelotas-profile-moho-known-depths.txt', moho_known, fmt='%.18f')<jupyter_code>## Initial guess data<jupyter_output><empty_output><jupyter_text>initial guess parameters vectorp0 = []p0 = np.vstack((ts1, tm, dS0))  prisms calculation by <NAME>aprism_w = fc.prism_w_function(xmax,xmin,dy,edge,dw,dcc,tw,yc)prism_s = fc.prism_s_function(xmax,xmin,dy,edge,ds,dcc,tw,p0,yc,ts0,two_layers=True)prism_c = fc.prism_c_function(xmax,xmin,dy,edge,S0,dcc,dc,tw,p0,yc,ts0,two_layers=True)prism_m = fc.prism_m_function(xmax,xmin,dy,edge,S0,dcc,dm,p0,yc)  z component of gravity calculation by <NAME>agzw = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_w)gzs0 = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_s[1])gzs1 = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_s[2])gzc = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_c)gzm = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_m)Observed data calculation:g0 = fc.g_function(x,yc,z,gzw,prism_s,prism_c,prism_m)g0 = gzw + gzs0 + gzs1 + gzc + gzm<jupyter_code>np.savetxt('../data/pelotas-profile-initial-guess-gravity-data.txt', g0, fmt='%.18f')
np.savetxt('../data/pelotas-profile-initial-guess-gravity-data-T07.txt', g0, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>Model plot<jupyter_code>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+dS0, S0+dS0]]).T,
                                                    props={'density': dm - dcc}))
%matplotlib inline
plt.close('all')
fig = plt.figure(figsize=(12,13))

import matplotlib.gridspec as gridspec
heights = [8, 8, 1]
gs = gridspec.GridSpec(3, 1, height_ratios=heights)
ax3 = plt.subplot(gs[0])
ax4 = plt.subplot(gs[1])
ax5 = plt.subplot(gs[2])

ax3.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='--', linewidth=1)
#ax3.plot(0.001*yc, gsyn, '-g', label='simulated data')
ax3.plot(0.001*yc, gobs, 'or', mfc='none', markersize=7, label='observed data')
ax3.plot(0.001*yc, g0, '-b', label='initial guess data')
ax3.set_xlim(0.001*ymin, 0.001*ymax)
ax3.set_ylabel('gravity disturbance (mGal)', fontsize=16)
ax3.set_xticklabels(['%g'% (l) for l in ax3.get_xticks()], fontsize=14)
ax3.set_yticklabels(['%g'% (l) for l in ax3.get_yticks()], fontsize=14)
ax3.legend(loc='best', fontsize=14, facecolor='silver')

ax4.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=1)
aux = yc <= COT
for (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])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.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])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.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])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.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])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.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])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.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])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='pink')

#ax4.axhline(y=S0, xmin=ymin, xmax=ymax, color='w', linestyle='--', linewidth=2)
ax4.plot(yc, tw, '-k', linewidth=3)
ax4.plot(yc, toi, '-k', linewidth=3)
ax4.plot(yc, interpreted_basement, '-k', linewidth=3, label='alternative interpretation surfaces')
ax4.plot(yc, interpreted_moho, '-k', linewidth=3)
ax4.plot(yc, basement, '-.b', linewidth=3, label='initial guess surfaces')
ax4.plot(yc, moho, '-.b', linewidth=3)
ax4.axhline(y=S0+true_dS0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=3)
ax4.axhline(y=S0+dS0, xmin=ymin, xmax=ymax, color='b', linestyle='-.', linewidth=3)
ax4.plot(base_known[:,0], base_known[:,1], 'v', color = 'yellow', markersize=15, label='known depths (basement)')
#ax4.plot(base_known_old[:,0], base_known_old[:,1], 'v', color = 'yellow', markersize=15, label='known depths (basement)')
#ax4.plot(base_known_new[:,0], base_known_new[:,1], 'v', color = 'magenta', markersize=15, label='more known depths (basement)')
ax4.plot(moho_known[:,0], moho_known[:,1], 'D', color = 'lime', markersize=15, label='known depths (moho)')
ax4.set_ylim((S0+dS0), zmin)
ax4.set_ylim((50500.0), zmin)
ax4.set_xlim(ymin, ymax)
ax4.set_xlabel('y (km)', fontsize=16)
ax4.set_ylabel('z (km)', fontsize=16)
ax4.set_xticklabels(['%g'% (0.001*l) for l in ax4.get_xticks()], fontsize=14)
ax4.set_yticklabels(['%g'% (0.001*l) for l in ax4.get_yticks()], fontsize=14)
ax4.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=18)
ax5.axis('off')

layers_list1 = ['water', 'sediment', 'SDR', 'continental', 'oceanic', 'mantle']
layers_list2 = ['', '', '', 'crust', 'crust', '']
colors_list = ['lightskyblue', 'tan', 'rosybrown', 'orange', 'olive', 'pink']
density_list = ['-1840', '-520', '-15', '0', '15', '370'] #original
#density_list = ['-1840', '-445', '-35', '0', '15', '370']
#density_list = ['1030', '2350', '2855', '2870', '2885', '3240']
ncols = len(colors_list)
nrows = 1
h = Y / nrows
w = X / (ncols + 1)

i=ncols-1
for 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.27
    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])
    ax5.plot(tmpx, tmpy, linestyle='-', color='k', linewidth=1)
    ax5.fill(tmpx, tmpy, color=color)
    ax5.text(x+w*0.375, yi_text1, layers1, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax5.text(x+w*0.375, yi_text2, layers2, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax5.text(x+w*0.375, yi_text3, density, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='center')

plt.tight_layout()    
    
#mpl.savefig('../manuscript/figures/pelotas-profile-initial-guess-model-grafics.png', dpi='figure', bbox_inches='tight')
#mpl.savefig('../manuscript/figures/pelotas-profile-initial-guess-model-grafics-T5-dsed-SDR.png', dpi='figure', bbox_inches='tight')
plt.show()      
%matplotlib inline
plt.close('all')
fig = plt.figure(figsize=(12,7))

import matplotlib.gridspec as gridspec
heights = [8, 1]
gs = gridspec.GridSpec(2, 1, height_ratios=heights)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

ax1.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=1)
aux = yc <= COT
for (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])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.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])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.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])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.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])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.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])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.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])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='pink')

#ax1.axhline(y=S0, xmin=ymin, xmax=ymax, color='w', linestyle='--', linewidth=2)
ax1.plot(yc, tw, '-k', linewidth=3)
ax1.plot(yc, toi, '-k', linewidth=3)
ax1.plot(yc, interpreted_basement, '-k', linewidth=3, label='alternative interpretation surfaces')
ax1.plot(yc, interpreted_moho, '-k', linewidth=3)
ax1.plot(yc, basement, '-.b', linewidth=3, label='initial guess surfaces')
ax1.plot(yc, moho, '-.b', linewidth=3)
ax1.axhline(y=S0+true_dS0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=3)
ax1.axhline(y=S0+dS0, xmin=ymin, xmax=ymax, color='b', linestyle='-.', linewidth=3)
ax1.plot(base_known[:,0], base_known[:,1], 'v', color = 'yellow', markersize=15, label='known depths (basement)')
#ax1.plot(base_known_old[:,0], base_known_old[:,1], 'v', color = 'yellow', markersize=15, label='known depths (basement)')
#ax1.plot(base_known_new[:,0], base_known_new[:,1], 'v', color = 'magenta', markersize=15, label='more known depths (basement)')
ax1.plot(moho_known[:,0], moho_known[:,1], 'D', color = 'lime', markersize=15, label='known depths (moho)')
ax1.set_ylim((S0+dS0), zmin)
ax1.set_ylim((50500.0), zmin)
ax1.set_xlim(ymin, ymax)
ax1.set_xlabel('y (km)', fontsize=16)
ax1.set_ylabel('z (km)', fontsize=16)
ax1.set_xticklabels(['%g'% (0.001*l) for l in ax1.get_xticks()], fontsize=14)
ax1.set_yticklabels(['%g'% (0.001*l) for l in ax1.get_yticks()], fontsize=14)
ax1.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=18)
ax2.axis('off')

layers_list1 = ['water', 'sediment', 'SDR', 'continental', 'oceanic', 'mantle']
layers_list2 = ['', '', '', 'crust', 'crust', '']
colors_list = ['lightskyblue', 'tan', 'rosybrown', 'orange', 'olive', 'pink']
density_list = ['-1840', '-520', '-15', '0', '15', '370'] #original
#density_list = ['-1840', '-430', '-40', '0', '15', '370']
#density_list = ['1030', '2350', '2855', '2870', '2885', '3240']

ncols = len(colors_list)
nrows = 1
h = Y / nrows
w = X / (ncols + 1)

i=ncols-1
for 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.27
    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])
    ax2.plot(tmpx, tmpy, linestyle='-', color='k', linewidth=1)
    ax2.fill(tmpx, tmpy, color=color)
    ax2.text(x+w*0.375, yi_text1, layers1, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax2.text(x+w*0.375, yi_text2, layers2, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax2.text(x+w*0.375, yi_text3, density, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='center')

plt.tight_layout()    
    
#mpl.savefig('../manuscript/figures/pelotas-profile-initial-guess-model.png', dpi='figure', bbox_inches='tight')
plt.show()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup