Fixed cookbook recipes and got rid of useless ones

Makefile

@@ -75,6 +75,7 @@ clean:
 	# The stuff fetched by the cookbook recipes
 	rm -rvf logo.png cookbook/logo.png
 	rm -rvf crust2.tar.gz cookbook/crust2.tar.gz
+	rm -rvf bouguer_alps_egm08.grd cookbook/bouguer_alps_egm08.grd
 
 clean-docs:
 	cd doc; make clean

cookbook/datasets_crust2.0_tesseroid.py

@@ -4,8 +4,7 @@
 """
 import time
 from multiprocessing import Pool
-from fatiando import gridder, utils, datasets
-from fatiando.gravmag import tesseroid
+from fatiando import datasets
 from fatiando.mesher import Tesseroid
 from fatiando.vis import mpl, myv
 
@@ -19,39 +18,3 @@
 myv.tesseroids(model, 'density')
 myv.continents(linewidth=3)
 myv.show()
-
-# Make the computation grid
-area = (-180, 180, -80, 80)
-shape = (100, 100)
-lons, lats, heights = gridder.regular(area, shape, z=250000)
-
-# Divide the model into nproc slices and calculate them in parallel
-
-
-def calculate(chunk):
-    return tesseroid.gz(lons, lats, heights, chunk)
-
-
-def split(model, nproc):
-    chunksize = len(model) / nproc
-    for i in xrange(nproc - 1):
-        yield model[i * chunksize: (i + 1) * chunksize]
-    yield model[(nproc - 1) * chunksize:]
-start = time.time()
-nproc = 8
-pool = Pool(processes=nproc)
-gz = sum(pool.map(calculate, split(model, nproc)))
-pool.close()
-print "Time it took: %s" % (utils.sec2hms(time.time() - start))
-
-mpl.figure(figsize=(10, 4))
-mpl.title('Crust gravity signal at 250km height')
-bm = mpl.basemap(area, 'robin')
-mpl.contourf(lons, lats, gz, shape, 35, basemap=bm)
-cb = mpl.colorbar()
-cb.set_label('mGal')
-bm.drawcoastlines()
-bm.drawmapboundary()
-bm.drawparallels(range(-90, 90, 45), labels=[0, 1, 0, 0])
-bm.drawmeridians(range(-180, 180, 60), labels=[0, 0, 0, 1])
-mpl.show()

cookbook/gravmag_eqlayer_joint.py

@@ -26,8 +26,7 @@
           + scale * EQLGravity(x2, y2, z2, gzz, layer, field='gzz'))
 regul = Smoothness2D(layer.shape)
 # Use an L-curve analysis to find the best regularization parameter
-solver = LCurve(
-    misfit, regul, [10 ** i for i in range(-30, -20)], jobs=8).fit()
+solver = LCurve(misfit, regul, [10 ** i for i in range(-30, -20)]).fit()
 layer.addprop('density', solver.estimate_)
 
 # Now I can forward model gz using my layer to produce an integrated map in a

cookbook/gravmag_harvester_grav.py

@@ -1,67 +1,70 @@
 """
 GravMag: 3D gravity inversion by planting anomalous densities using
-``harvester`` (simple example)
+``harvester`` (more complex interactive example)
 """
 from fatiando import gridder, utils
-from fatiando.gravmag import prism, harvester
-from fatiando.mesher import Prism, PrismMesh, vremove
+from fatiando.gravmag import polyprism, harvester
+from fatiando.mesher import PolygonalPrism, PrismMesh, vremove
 from fatiando.vis import mpl, myv
 
 # Create a synthetic model
-model = [Prism(250, 750, 250, 750, 200, 700, {'density': 1000})]
+bounds = [-10000, 10000, -10000, 10000, 0, 10000]
+vertices = [[-4948.97959184, -6714.64019851],
+            [-2448.97959184, -3141.43920596],
+            [2448.97959184,   312.65508685],
+            [6938.7755102,  5394.54094293],
+            [4846.93877551,  6228.28784119],
+            [2653.06122449,  3409.4292804],
+            [-3520.40816327, -1434.24317618],
+            [-6632.65306122, -6079.4044665]]
+model = [PolygonalPrism(vertices, 1000, 4000, {'density': 1000})]
 # and generate synthetic data from it
-shape = (25, 25)
-bounds = [0, 1000, 0, 1000, 0, 1000]
+shape = (20, 20)
 area = bounds[0:4]
 xp, yp, zp = gridder.regular(area, shape, z=-1)
 noise = 0.1  # 0.1 mGal noise
-gz = utils.contaminate(prism.gz(xp, yp, zp, model), noise)
-# plot the data
+gz = utils.contaminate(polyprism.gz(xp, yp, zp, model), noise)
+
+# Create a mesh
+mesh = PrismMesh(bounds, (25, 50, 50))
+# Wrap the data so that harvester can read it
+data = [harvester.Gz(xp, yp, zp, gz)]
+# Plot the data and pick the location of the seeds
 mpl.figure()
-mpl.title("Synthetic gravity anomaly (mGal)")
+mpl.suptitle("Pick the seeds (polygon is the true source)")
 mpl.axis('scaled')
 levels = mpl.contourf(yp, xp, gz, shape, 12)
 mpl.colorbar()
+mpl.polygon(model[0], xy2ne=True)
 mpl.xlabel('Horizontal coordinate y (km)')
 mpl.ylabel('Horizontal coordinate x (km)')
-mpl.m2km()
+seedx, seedy = mpl.pick_points(area, mpl.gca(), xy2ne=True).T
+# Set the right density and depth
+locations = [[x, y, 1500, {'density': 1000}] for x, y in zip(seedx, seedy)]
 mpl.show()
-
-# Inversion setup
-# Create a mesh
-mesh = PrismMesh(bounds, (25, 25, 25))
-# Wrap the data so that harvester can use it
-data = [harvester.Gz(xp, yp, zp, gz)]
-# Make the seed
-seeds = harvester.sow([[500, 500, 450, {'density': 1000}]], mesh)
-# Run the inversioin
+# Make the seed and set the compactness regularizing parameter mu
+seeds = harvester.sow(locations, mesh)
+# Run the inversion
 estimate, predicted = harvester.harvest(data, seeds, mesh,
-                                        compactness=0.5, threshold=0.0005)
-
+                                        compactness=0.05, threshold=0.0005)
 # Put the estimated density values in the mesh
 mesh.addprop('density', estimate['density'])
-
-# Plot the adjustment
+# Plot the adjustment and the result
 mpl.figure()
-mpl.title("True: color | Inversion: contour")
+mpl.title("True: color | Predicted: contour")
 mpl.axis('scaled')
 levels = mpl.contourf(yp, xp, gz, shape, 12)
 mpl.colorbar()
 mpl.contour(yp, xp, predicted[0], shape, levels, color='k')
 mpl.xlabel('Horizontal coordinate y (km)')
 mpl.ylabel('Horizontal coordinate x (km)')
 mpl.m2km()
-residuals = gz - predicted[0]
-mpl.figure()
-mpl.title('Residuals: mean=%g stddev=%g' % (residuals.mean(), residuals.std()))
-mpl.hist(residuals, bins=10)
-mpl.xlabel('Residuals (mGal)')
-mpl.ylabel('# of')
 mpl.show()
 # Plot the result
 myv.figure()
-myv.prisms(model, 'density', style='wireframe')
+myv.polyprisms(model, 'density', opacity=0.6, linewidth=5)
 myv.prisms(vremove(0, 'density', mesh), 'density')
+myv.prisms(seeds, 'density')
 myv.axes(myv.outline(bounds), ranges=[i * 0.001 for i in bounds], fmt='%.1f',
          nlabels=6)
 myv.wall_bottom(bounds)