from IPython.display import Image

from fatiando import gridder
from fatiando.vis import mpl
area = [-20, 20, -50, 50]
x, y = gridder.scatter(area, n=100)
data = x**2 + y**2
mpl.figure()
mpl.axis('scaled')
mpl.contourf(y, x, data, shape=(50, 50), 
    levels=30, interp=True)
mpl.colorbar(orientation='horizontal')
mpl.plot(y, x, '.k')
mpl.xlabel('y (East-West)')
mpl.ylabel('x (North-South)')
#mpl.savefig('gridding_plotting_contourf.png')
mpl.show()

Image(filename='gridding_plotting_contourf.png')

mpl.figure()
bm = mpl.basemap(area, projection='robin')
bm.drawmapboundary()
bm.drawcoastlines()
mpl.contourf(y, x, data, shape=(50, 50), 
    levels=30, interp=True, basemap=bm)
mpl.colorbar(orientation='horizontal')
#mpl.savefig('gridding_plotting_basemap.png')
mpl.show()

Image(filename='gridding_plotting_basemap.png')

from fatiando import mesher
model = [
    mesher.Prism(5, 8, 3, 7, 1, 7, 
        props={'density':200}),    
    mesher.Prism(1, 2, 4, 5, 1, 2, 
        props={'density':1000})]

from fatiando.vis import myv
bounds = [0, 10, 0, 10, 0, 10]
myv.figure()
myv.prisms(model, 'density')
myv.axes(myv.outline(bounds))
myv.wall_bottom(bounds)
myv.wall_north(bounds)
#myv.savefig('meshes_3dplotting_2prisms.png')
myv.show()

Image(filename='meshes_3dplotting_2prisms.png')

mesh = mesher.PrismMesh(bounds, shape=(3, 3, 3))
mesh.addprop('density', range(mesh.size))
myv.figure()
myv.prisms(mesh, 'density')
myv.axes(myv.outline(bounds))
#myv.savefig('meshes_3dplotting_mesh.png')
myv.show()

Image(filename='meshes_3dplotting_mesh.png')

from fatiando import utils
x, y = gridder.regular(bounds[:4], (50, 50))
heights = -5 + 5*utils.gaussian2d(x, y, 10, 5, 
    x0=10, y0=10)
mesh = mesher.PrismMesh(bounds, (20, 20, 20))
mesh.addprop('density', range(mesh.size))
mesh.carvetopo(x, y, heights)
myv.figure()
myv.prisms(mesh, 'density')
myv.axes(myv.outline(bounds))
myv.wall_north(bounds)
#myv.savefig('meshes_3dplotting_meshtopo.png')
myv.show()

Image(filename='meshes_3dplotting_meshtopo.png')

model = [
    mesher.Tesseroid(-60, -55, -30, -27, 500000, 0, 
        props={'density':200}),
    mesher.Tesseroid(-66, -55, -20, -10, 300000, 0, 
        props={'density':-100})]
fig = myv.figure(zdown=False)
myv.tesseroids(model, 'density')
myv.continents(linewidth=2)
myv.earth(opacity=1)
myv.meridians(range(0, 360, 45), opacity=0.2)
myv.parallels(range(-90, 90, 45), opacity=0.2)
# Rotate the camera to get a good view
scene = fig.scene
scene.camera.position = [21199620.406122234, 
    -12390254.839673528, -14693312.866768979]
scene.camera.focal_point = [-535799.97230670298, 
    -774902.33205294283, 826712.82283183688]
scene.camera.view_angle = 19.199999999999996
scene.camera.view_up = [0.33256519487680014, 
    -0.47008782429014295, 0.81756824095039038]
scene.camera.clipping_range = [7009580.0037488714, 
    55829873.658824757]
scene.camera.compute_view_plane_normal()
scene.render()
#myv.savefig('meshes_3dplotting_tesseroid.png')
myv.show()

Image(filename='meshes_3dplotting_tesseroid.png')

from fatiando import gravmag
area = [-80, -30, -40, 10]
shape = (50, 50)
lons, lats, heights = gridder.regular(area, shape, 
    z=2500000)
gz = gravmag.tesseroid.gz(lons, lats, heights, model)
mpl.figure()
bm = mpl.basemap(area, 'ortho')
bm.drawcoastlines()
bm.drawmapboundary()
bm.bluemarble()
mpl.title('Gravity anomaly (mGal)')
mpl.contourf(lons, lats, gz, shape, 30, basemap=bm)    
mpl.colorbar()
#mpl.savefig('gravmag_tesseroid_data.png')
mpl.show()

Image(filename='gravmag_tesseroid_data.png')

# Draw a polygon and make a polygonal prism
bounds = [-1000, 1000, -1000, 1000, 0, 1000]
area = bounds[:4]
mpl.figure()
mpl.axis('scaled')
vertices = mpl.draw_polygon(area, mpl.gca(), 
    xy2ne=True)
model = [mesher.PolygonalPrism(vertices, z1=0, 
    z2=500, props={'density':500})]
# Calculate the gravity anomaly
shape = (50, 50)
x, y, z = gridder.regular(area, shape, z=-1)
gz = gravmag.polyprism.gz(x, y, z, model)
mpl.figure()
mpl.axis('scaled')
mpl.title("Gravity anomaly (mGal)")
mpl.contourf(y, x, gz, shape, levels=30)
mpl.colorbar()
mpl.polygon(model[0], '.-k', xy2ne=True)
mpl.set_area(area)
mpl.m2km()
#mpl.savefig('forward_modeling_polyprism_data.png')
mpl.show()
myv.figure()
myv.polyprisms(model, 'density')
myv.axes(myv.outline(bounds), 
    ranges=[i*0.001 for i in bounds])
myv.wall_north(bounds)
myv.wall_bottom(bounds)
#myv.savefig('forward_modeling_polyprism_model.png')
myv.show()


Image(filename='forward_modeling_polyprism_model.png')

Image(filename='forward_modeling_polyprism_data.png')

# This cell is not include in the proceedings. 
# It saves the output of the previous cell so that I could run the code bellow with the same model if the notebook is restarted.
import cPickle as pickle
with open('forward_modeling_model.pickle', 'w') as f:
    pickle.dump([model, x, y, z, gz], f)

# This cell is not included in the proceedings.
# It loads the model saved in the previous cell.
# If you want to use the model shown in the proceedings, put the pickle file provided in the notebook directory and run this cell.
import cPickle as pickle
with open('forward_modeling_model.pickle') as f:
    model, x, y, z, gz = pickle.load(f)    
    bounds = [-1000, 1000, -1000, 1000, 0, 1000]
    area = bounds[:4]
    from fatiando import gravmag, gridder, mesher
    from fatiando.vis import mpl, myv
    shape = (50, 50)

estimate = gravmag.imaging.sandwich(x, y, z, gz, 
    shape, zmin=0, zmax=1000, nlayers=20, 
    power=0.1)
body = mesher.vfilter(1.3*10**8, 1.8*10**8,
    'density', estimate)
myv.figure()
myv.prisms(body, 'density', edges=False)
p = myv.polyprisms(model, 'density', 
    style='wireframe', linewidth=4)
p.actor.mapper.scalar_visibility = False
p.actor.property.color = (0, 0, 0)
myv.axes(myv.outline(bounds), 
    ranges=[i*0.001 for i in bounds])
myv.wall_north(bounds)
myv.wall_bottom(bounds)
#myv.savefig('gravmag_imaging.png')
myv.show()

Image(filename='gravmag_imaging.png')

# Make a mesh and a seed
mesh = mesher.PrismMesh(bounds, (15, 30, 30))
seeds = gravmag.harvester.sow([[200, 300, 100, 
    {'density':500}]], mesh)
myv.figure()
myv.prisms([mesh[s.i] for s in seeds])
p = myv.polyprisms(model, 'density', 
    style='wireframe', linewidth=3)
p.actor.mapper.scalar_visibility = False
p.actor.property.color = (0, 0, 0)
myv.axes(myv.outline(bounds), 
    ranges=[i*0.001 for i in bounds])
myv.wall_north(bounds)
myv.wall_bottom(bounds)
#myv.savefig('gravmag_harvester_seed.png')
myv.show()
# Now perform the inversion
data = [gravmag.harvester.Gz(x, y, z, gz)]
estimate = gravmag.harvester.harvest(data, seeds, 
    mesh, compactness=0.1, threshold=0.0001)[0]
mesh.addprop('density', estimate['density'])
body = mesher.vremove(0, 'density', mesh)
myv.figure()
myv.prisms(body, 'density')
p = myv.polyprisms(model, 'density', 
    style='wireframe', linewidth=4)
p.actor.mapper.scalar_visibility = False
p.actor.property.color = (0, 0, 0)
myv.axes(myv.outline(bounds), 
    ranges=[i*0.001 for i in bounds])
myv.wall_north(bounds)
myv.wall_bottom(bounds)
#myv.savefig('gravmag_harvester.png')
myv.show()
        

Image(filename='gravmag_harvester_seed.png')

Image(filename='gravmag_harvester.png')

import urllib
from fatiando import mesher, utils, seismic
from fatiando.vis import mpl
area = (0, 500000, 0, 500000)
shape = (30, 30)
model = mesher.SquareMesh(area, shape)
link = '/'.join(["http://fatiando.readthedocs.org",
    "en/Version0.1/_static/logo.png"])
urllib.urlretrieve(link, 'model.png')
model.img2prop('model.png', 4000, 10000, 'vp')
quake_locations = utils.random_points(area, 40)
receiver_locations = utils.circular_points(area, 20,
    random=True)
quakes, receivers = utils.connect_points(
    quake_locations, receiver_locations)
traveltimes = seismic.ttime2d.straight(model, 'vp', 
    quakes, receivers)
noisy = utils.contaminate(traveltimes, 0.001, 
    percent=True)

Image(filename='logo.png')

mesh = mesher.SquareMesh(area, shape)
slowness, residuals = seismic.srtomo.run(noisy, 
    quakes, receivers, mesh, smooth=10**6)
velocity = seismic.srtomo.slowness2vel(slowness)
mesh.addprop('vp', velocity)
# Make the plots
mpl.figure(figsize=(9, 7))
mpl.subplots_adjust(top=0.95, bottom=0.05, 
    left=0.05, right=0.95)
mpl.subplot(2, 2, 1)
mpl.title('Velocity model (m/s)')
mpl.axis('scaled')
mpl.squaremesh(model, prop='vp', cmap=mpl.cm.seismic)
mpl.colorbar(pad=0.01)
mpl.points(quakes, '*y', label="Sources")
mpl.points(receivers, '^g', label="Receivers")
mpl.m2km()
mpl.subplot(2, 2, 2)
mpl.title('Ray paths')
mpl.axis('scaled')
mpl.squaremesh(model, prop='vp', cmap=mpl.cm.seismic)
mpl.colorbar(pad=0.01)
mpl.paths(quakes, receivers)
mpl.points(quakes, '*y', label="Sources")
mpl.points(receivers, '^g', label="Receivers")
mpl.m2km()
mpl.subplot(2, 2, 3)
mpl.title('Estimated velocity (m/s)')
mpl.axis('scaled')
mpl.squaremesh(mesh, prop='vp', cmap=mpl.cm.seismic, 
    vmin=4000, vmax=10000)
mpl.colorbar(pad=0.01)
mpl.m2km()
mpl.subplot(2, 2, 4)
mpl.title('Residuals (s)')
mpl.hist(residuals, bins=10)
#mpl.savefig('seismic_tomo.png')
mpl.show()

Image(filename='seismic_tomo.png')