# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# # A mesh class for a tesseroid relief 
#
# There is a `PrismRelief` mesh class in [Fatiando a Terra](http://www.fatiando.org/) but we still need a `TesseroidRelief` for this inversion. This is a mesh of tesseroids distributed along an area. They undulate below and above a reference level, describing the relief of an interface. Tesseroids have either the top of bottom fixed to a reference height. The other end undulates along a relief. The `TesseroidRelief` class is defined in the [`mohoinv.py`](mohoinv.py) module. This notebook will show some of the features of the mesh and how to use it.

# ## Package imports

# Insert plots into the notebook
# %matplotlib inline

from __future__ import division, unicode_literals
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image
import multiprocessing
import seaborn  # Makes the default style of the plots nicer

# Load the required modules from Fatiando a Terra and show the specific version of the library used.

from fatiando import gridder, utils
from fatiando.gravmag import tesseroid
import fatiando

print("Using Fatiando a Terra version: {}".format(fatiando.__version__))

from mohoinv import TesseroidRelief

# ## Create some synthetic relief

# Define a regular grid.

# shape is nlat, nlon = the number of points in the grid
shape = (41, 31)
# Make a regular grid inside an area = (s, n, w, e)
area = (20, 60, -40, 40)
lat, lon, h = gridder.regular(area, shape, z=250e3)

# The model area is slightly larger because the points generated above are in the center of each cell.

dlat, dlon = gridder.spacing(area, shape)
s, n, w, e = area
modelarea = (s - dlat/2, n + dlat/2, w - dlon/2, e + dlon/2)

# Make a checker board relief undulating along a specified height reference.

f = 0.2
reference = -35e3
relief = 10e3*np.sin(0.5*f*lon)*np.cos(f*lat) + reference

plt.figure(figsize=(7, 3))
plt.title('Synthetic relief')
plt.axis('scaled')
plt.pcolormesh(lon.reshape(shape), lat.reshape(shape), relief.reshape(shape),
               cmap="RdYlBu_r")
plt.colorbar(pad=0.01).set_label('meters')
plt.xlim(lon.min(), lon.max())
plt.ylim(lat.min(), lat.max())
plt.tight_layout()

# Set a density contrast for the relief. The density contrast is negative if the relief is below the reference and positive otherwise.

density = 600*np.ones_like(relief)
density[relief < reference] *= -1

plt.figure(figsize=(7, 3))
plt.title('Relief density contrast')
plt.axis('scaled')
plt.pcolormesh(lon.reshape(shape), lat.reshape(shape), density.reshape(shape),
               cmap="RdBu_r")
plt.colorbar(pad=0.01).set_label(u'kg/m³')
plt.xlim(lon.min(), lon.max())
plt.ylim(lat.min(), lat.max())
plt.tight_layout()

# Now we can create a mesh.

sample_mesh = TesseroidRelief(modelarea, shape, relief, reference, {'density': density})

# ## Calculate the gravitational effect of this mesh
#
# The mesh behaves like a list of `Tesseroid` objects. So we can pass it to any function in Fatiando a Terra that takes such list as input.
#
# Below, we'll show an example of using the forward modeling functions in Fatiando to calculate the gravitational effect of the above relief.

ncpu = multiprocessing.cpu_count()
print('Number of cores: {}'.format(ncpu))

# Forward model the data on the grid generated above using  all available cores of the processor.

data = tesseroid.gz(lon, lat, h, sample_mesh, njobs=ncpu)

# The warnings above are because some small tesseroids (below 10 cm dimensions) are ignored to avoid numerical issues.

plt.figure(figsize=(7, 3))
plt.title('Relief gravity anomaly')
plt.axis('scaled')
plt.tricontourf(lon, lat, data, 40, cmap="RdBu_r")
plt.colorbar(pad=0.01).set_label(u'mGal')
plt.xlim(lon.min(), lon.max())
plt.ylim(lat.min(), lat.max())
plt.tight_layout()

StarcoderdataJupyterScriptsDedupFiltered