# Laying out a seismic survey¶

This notebook goes with a blog post — Laying out a seismic survey — on Wednesday 17 December 2014.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
from geopandas import GeoDataFrame

%matplotlib inline


## Set up a basic geometry¶

Next we can define the survey parameters.

In [2]:
# initial conditions
xmi = 575000   # leftmost corner of grid (m)
ymi= 4710000   # bottommost corner of grid (m)
SL = 600       # shot line interval (m)
RL = 600       # receiver line interval (m)
si = 100       # source point interval (m)
ri = 100       # receiver interval (m)
x = 3000       # width of survey (m)
y = 1800       # height of survey (m)


Some easy geometry: calculate the number of source and receiver lines we can fit into the shape we have defined, given the spacing we want:

In [3]:
# Number of receiver lines and source lines
rlines = int(y/RL) + 1
slines = int(x/SL) + 1

# Put recevier lines East-West, and shot lines North South
rperline = int(x/ri) + 2
sperline = int(y/si) + 2

# offset the receivers relative to the sources
shiftx = -si/2.
shifty = -ri/2.


This is a double list comprehension. It's a bit sticky. First, to illustrate the point, see how this one-liner makes a simple list:

In [4]:
[x**2 for x in range(10)]

Out[4]:
[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

Now check this out:

In [5]:
[x**y for x in range(4) for y in range(4)]

Out[5]:
[1, 0, 0, 0, 1, 1, 1, 1, 1, 2, 4, 8, 1, 3, 9, 27]

Now hopefully it's not too hard to see what we're doing here...

In [6]:
# Find x,y coordinates of recs and shots
rcvrx = [xmi + rcvr*ri + shiftx for line in range(rlines) for rcvr in range(rperline)]
rcvry = [ymi + line*RL - shifty for line in range(rlines) for rcvr in range(rperline)]

srcx = [xmi + line*SL for line in range(slines) for src in range(sperline)]
srcy = [ymi + src*si for line in range(slines) for src in range(sperline)]

In [7]:
def plot_geoms(xcoords, ycoords, color = 'none', size = 50, alpha = 0.5):
"""
A helper function to make it a bit easier to plot multiple things.
"""
plot = plt.scatter(xcoords,
ycoords,
c=color,
s=size,
marker='o',
alpha=alpha,
edgecolor='none'
)
return plot

In [8]:
fig = plt.figure(figsize=(15,10))
r = plot_geoms(rcvrx, rcvry, 'b')
s = plot_geoms(srcx, srcy, 'r')
plt.axis('equal')
plt.grid()


## Making a more robust data container¶

In [9]:
# Zip into x,y coordinate pairs.
rcvrxy = zip(rcvrx, rcvry)
srcxy = zip(srcx, srcy)

# Create a list of Point objects.
rcvrs = [Point(x,y) for x,y in rcvrxy]
srcs = [Point(x,y) for x,y in srcxy]

# Put lists into a GeoDataFrame.
sources = GeoDataFrame({'geometry': srcs})


If you want to, save the GeoDataFrames as shapefiles:

In [10]:
sources.to_file('sources.shp')


## Computing the midpoints¶

We're going to need Shapely's LineString geometry — essentially we ask for the 50% (0.5) mark along a line interpolated between two Point objects. The interpolate method returns another Point.

In [11]:
from shapely.geometry import LineString


Another nested list comp gets the job done. Python is awesome!

In [12]:
midpoints = [LineString([r, s]).interpolate(0.5, normalized=True) for r in rcvrs for s in srcs]


We won't get into this too much right now, but let's just check we got what we wanted. We can throw the midpoints into a GeoPandas GeoSeries object, essentially a 1D GeoDataFrame — a bit like a 1D labelled array — then use its built-in plot method.

In [13]:
from geopandas import GeoSeries

In [14]:
m = GeoSeries(midpoints)
m.plot()

Out[14]:
<matplotlib.axes.AxesSubplot at 0x10ee96c10>

Voilà!