Preliminaries - import dependencies



In [ ]:

    
%matplotlib inline

import numpy as np
import scipy as sp  
import matplotlib.pyplot as plt
from linecache import getline

Import a Petrel surface as IRAP grid file.

NB: IRAP grid files have a 4-line header and then Z values arranged in rows of 6. In a few instances, however, the IRAP grid exported did not have 6 values in the last row, in which case this method will break, unless NaNs are added (by hand or other means).

The surface was exported from Petrel as an IRAP grid; we will get some basic information from the file header The header looks like this: -nnn nInlines Xspacing Yspacing Xmin Xmax Ymin Ymax nXlines 0.000000 X Y 0 0 0 0 0 0 0



In [ ]:

    
# Import the grid, skipping the 4-row header
grid  = np.loadtxt("\\your_directory_path\grid_IRAP", skiprows=4)



In [ ]:

    
# Calculate the number of rows and columns in the raw file and total number of points
r, c = grid.shape
npoints = r*c

Get nInlines and nXlines, reshape grid and plot



In [ ]:

    
# get the string containing nInlines from line one of the header
grid1 = getline("grid_IRAP", 1)
grid2 = getline("grid_IRAP", 2)

# split values in the string and convert them to floating point numbers
g1a, g1b, g1c, g1d = (grid1.split(" "))
nInlines = float(g1b)

g2a, g2b, g2c, g2d = (grid2.split(" "))
nXlines = float(g2a)



In [ ]:

    
# It looks like we have 'npoints' grid points in total, organized in 'r' rows with c (this is always 6) values in each
# From the header we see that the grid had a size/shape of nInlines x nXlines.  
# Do we have all the points?

assert nInlines * nXlines == npoints
# it checks



In [ ]:

    
# Reshape the points into an r*c columns rectangular array 
newgrid = np.reshape(grid, (nInlines, nXlines))

# check the shape
newgrid.shape



In [ ]:

    
# N.B.
# use origin 'lower' or the image will be upside down (origin of an image is the top left)

# let's take a look
fig = plt.figure(figsize=(12, 7))

ax = fig.add_subplot(1, 1, 1)
ax.set_xticks([])
ax.set_yticks([])

plt.imshow(newgrid, cmap='gray_r',origin='lower')
plt.show()



In [ ]:

    
# Check for null values/invalid data
max = np.amax(newgrid) 
min = np.amin(newgrid)
print max, min

# it looks like the max value, 9999900.0, is the null



In [ ]:

    
# To eliminate the null value points, first we create a masked array 
# This way those values are not taken into account when calculating the mean value
finalgrid = np.ma.masked_where(newgrid == max, newgrid)
fgridmean = finalgrid.mean()

# Then we replace them with the mean value
finalgrid = finalgrid.filled(fgridmean)



In [ ]:

    
# Let's take a look again
fig = plt.figure(figsize=(12, 7))

ax = fig.add_subplot(1, 1, 1)
ax.set_xticks([])
ax.set_yticks([])

plt.imshow(finalgrid, cmap='gray_r', origin='lower')
plt.show()

Plotting grid with coordinates and contour overlay



In [ ]:

    
# get the string containing X,Y limits from line two of the header
grid_space = getline("grid_IRAP, 2)



In [ ]:

    
# split values in the string and conver them to floating point numbers
xmin, xmax, ymin, ymax = (grid_space.split(" "))
xmin, xmax, ymin, ymax = float(xmin), float(xmax), float(ymin), float(ymax)



In [ ]:

    
# create X, Y arrays
r1, c1 = finalgrid.shape

x = np.linspace(xmin,xmax,c1)
y = np.linspace(ymin,ymax,r1)



In [ ]:

    
# plot distance image versus X,Y coordinates and add contours

fig = plt.figure(figsize=(18, 12))

ax = fig.add_subplot(1, 1, 1)

plt.imshow(finalgrid, cmap='viridis',
           extent=[xmin, xmax, ymin, ymax], origin='lower')

levels = np.arange(5., 90., 10.)
CS = plt.contour(finalgrid, levels, colors='k',
                 origin='lower',
                 linewidths=1,
                 extent=[xmin, xmax, ymin, ymax])

plt.show()

# some useful reference:
# http://matplotlib.org/examples/pylab_examples/pcolor_demo.html
# http://matplotlib.org/examples/pylab_examples/contour_demo.html

Do something with the grid



In [ ]:

    
## filtered_grid = bla bla bla

Export results



In [ ]:

    
# 1 - Reverse reshape 
# 2 - Replace final values in the original grid, add headers, and export



In [ ]:

    
# 1 - Reverse reshape 
outgrid = np.reshape(filtered_grid, (r,c))

# Check
outgrid.shape



In [ ]:

    
# 2 - Copy the original grid, replace all non null values in it with distance values
output = grid
output[grid != max] = outgrid[grid != max]



In [ ]:

    
# Create output header

header =  "%s%s" % (getline("grid_IRAP", 1), "\n") 
header += "%s%s" % (getline("grid_IRAP", 2), "\n") 
header += "%s%s" % (getline("grid_IRAP", 3), "\n") 
header += "%s%s" % (getline("grid_IRAP", 4), "\n")



In [ ]:

    
# export the final map
# Open output file, write header to it, then add distance grid
           
with open("\\your_directory_path", "w") as f:
    f.write(header)
    np.savetxt(f,output,delimiter=' ', fmt='%f')