Esse exemplo foi tirado do cookbook: http://fatiando.readthedocs.org/en/latest/cookbook/mesher_prismmesh_topo.html



In [1]:

    
from fatiando import gridder, utils, mesher
from fatiando.vis import myv
from IPython.display import Image



In [2]:

    
x1, x2 = -100, 100
y1, y2 = -200, 200
bounds = (x1, x2, y1, y2, -200, 0)

x, y = gridder.regular((x1, x2, y1, y2), (50, 50))
height = (100 +
          -50 * utils.gaussian2d(x, y, 100, 200, x0=-50, y0=-100, angle=-60) +
          100 * utils.gaussian2d(x, y, 50, 100, x0=80, y0=170))

mesh = mesher.PrismMesh(bounds, (20, 40, 20))
mesh.carvetopo(x, y, height)



In [3]:

    
scene = myv.figure()
myv.prisms(mesh, linewidth=1)
myv.axes(myv.outline(bounds), fmt='%.0f')
myv.wall_north(bounds)
scene.scene.camera.position = [733.0947480478585, 519.11700319414456, -400.5263974989316]
scene.scene.camera.focal_point = [-1.0191166350008132, 22.817088165087696, -68.178053379279476]
scene.scene.camera.view_angle = 24.0
scene.scene.camera.view_up = [-0.30212417216680404, -0.17998297427981638, -0.93612344995834529]
scene.scene.camera.clipping_range = [503.99070671659922, 1505.6079017359789]
scene.scene.camera.compute_view_plane_normal()
scene.scene.render()
myv.savefig('mesh.png')
myv.show()









    



/home/leo/bin/anaconda/lib/python2.7/site-packages/mayavi/preferences/preference_manager.py:24: UserWarning: Module fatiando was already imported from /home/leo/src/fatiando/fatiando/__init__.pyc, but /home/leo/bin/anaconda/lib/python2.7/site-packages is being added to sys.path
  import pkg_resources
WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.



In [4]:

    
Image(filename='mesh.png')









    Out[4]:

O carvetopo cria uma variável mask no mesh. Esse mask é uma lista de prismas que devem ser desconsiderados. É assim que as funções que plotam e calculam sabem que devem ignorar aquela parte do mesh.



In [5]:

    
mesh.mask[-10:]  # Últimos 10 prismas que devem ser ignorados (estão acima da topografia)









    Out[5]:





[11443, 11444, 11445, 11446, 11447, 11448, 11465, 11466, 11467, 11468]

Quando algum código pede um elemento do mesh que está no mask, ele retorna None:



In [6]:

    
print(mesh[11443], mesh[11444])









    



(None, None)

Para inverter a topografia, você precisa inverter o mask. No Python existe o tipo set, que representa conjuntos. Eles são úteis para fazer coisas como união, diferença, intersecção, etc.



In [7]:

    
a = set(range(10))
a









    Out[7]:





{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}



In [8]:

    
a.difference([5, 6])









    Out[8]:





{0, 1, 2, 3, 4, 7, 8, 9}



In [9]:

    
antimask = list(set(xrange(mesh.size)).difference(mesh.mask))



In [10]:

    
antimask[:10]









    Out[10]:





[717, 718, 737, 738, 757, 758, 777, 778, 1477, 1478]



In [11]:

    
mesh.mask = antimask



In [12]:

    
scene = myv.figure()
myv.prisms(mesh, linewidth=1)
myv.axes(myv.outline(bounds), fmt='%.0f')
myv.wall_north(bounds)
scene.scene.camera.position = [733.0947480478585, 519.11700319414456, -400.5263974989316]
scene.scene.camera.focal_point = [-1.0191166350008132, 22.817088165087696, -68.178053379279476]
scene.scene.camera.view_angle = 24.0
scene.scene.camera.view_up = [-0.30212417216680404, -0.17998297427981638, -0.93612344995834529]
scene.scene.camera.clipping_range = [503.99070671659922, 1505.6079017359789]
scene.scene.camera.compute_view_plane_normal()
scene.scene.render()
myv.savefig('mesh-anti.png')
myv.show()



In [13]:

    
Image(filename='mesh-anti.png')









    Out[13]: