In [1]:
import re
import numpy as np
from io import StringIO
In [2]:
import geopandas as gp
from shapely.geometry import Point, LineString
In [3]:
import matplotlib.pyplot as plt
%matplotlib inline
In [4]:
points_file = "../data/Petrel_points_extract.txt"
Let's look at the start of the file:
In [5]:
with open(points_file) as f:
c = f.read()
print(c[:256])
We can make a function to parse the header into some variables, and read the bulk of the file – the data — into a NumPy array.
In [6]:
def gdf_from_points(points_file):
"""
Read a Petrel points file and create a GeoPandas DataFrame.
"""
with open(points_file) as f:
comments, fields = [], []
in_header = False
while True:
line = f.readline().strip()
if line.startswith('#'):
comments.append(line.strip('# '))
elif line.startswith('VERSION'):
version = line.split()[-1]
elif line.startswith('BEGIN'):
in_header = True
elif line.startswith('END'):
in_header = False
break
elif in_header:
fields.append(line.strip())
else:
break
d = f.read()
s = StringIO(d)
data = np.loadtxt(s)
geoseries = gp.GeoSeries([Point(r[0], r[1]).buffer(5) for r in data])
d = {'geometry': geoseries, 'twt': -data[:,2]}
return gp.GeoDataFrame(d)
Read the subset of data:
In [7]:
seafloor = gdf_from_points(points_file)
In [9]:
seafloor.plot()
plt.show()
In [10]:
horizon_file = "../data/Petrel_IESX2D_horizon_extract.txt"
In [12]:
with open(horizon_file) as f:
c = f.read()
print(c[:512])
Not totally sure what all these fields are. There seems to be two rows of header info for each horizon. This is followed by space-delimited data fields with columns as follows:
I have asked about this on Stack Exchange and on LinkedIn.
In [13]:
import re
s = "PROFILE Fault to seafloor TYPE 1 5 By Petrel 2014.2 (64-bit) 2d_ci7m_gf.ifdf m ms"
re.search(r'PROFILE (.+?) +TYPE', s).group(1)
Out[13]:
I think we can use the first integer as a sort of flag to determine the start of a new 'segment'.
I want to make a shapefile with line segments, I think. Not totally sure about that...
In [19]:
import operator, re
def gdf_from_iesx(filename, threed=False):
"""
Read a Petrel IESX file and create a GeoPandas DataFrame.
"""
with open(filename) as f:
points, linestrings, names = [], [], []
mins, maxs = [], []
min_sfls, max_sfls = [], []
minx, miny = [], []
last_cdp = 0
skip = False
while True:
line = f.readline().strip()
if not line:
# End of file
break
elif line.startswith('EOD'):
# End of horizon
last_cdp = 0 # Force capture
elif line.startswith('SNAPPING'):
continue
elif line.startswith('PROFILE'):
name = re.search(r'PROFILE (.+?) +TYPE', line).group(1)
# Some 'label' horizons slipped though, skip 'em.
if name.startswith('---'):
skip = True
else:
skip = False
print()
print(name, end="")
else:
if skip == True:
continue
line = line.split()
x, y = float(line[0]), float(line[1])
twtt = float(line[4])
if threed:
this_cdp = int(line[5]) + int(line[9])
else:
this_cdp = int(line[7])
if abs(this_cdp - last_cdp) < 2:
# Then it's a regular line, so keep adding
points.append(Point(x, y, twtt))
last_cdp = this_cdp
else:
if len(points) < 2:
last_cdp = this_cdp
continue
print('.', end="")
# Capture what we have
linestrings.append(LineString(points))
names.append(name)
zs = [p.z for p in points]
# We want the value and index of the min and max in the list of z values
min_idx, min_val = min(enumerate(zs), key=operator.itemgetter(1))
max_idx, max_val = max(enumerate(zs), key=operator.itemgetter(1))
# Record the values and shallowest point in our master lists
mins.append(min_val)
maxs.append(max_val)
p = points[min_idx]
minx.append(p.x)
miny.append(p.y)
# Make GeoDataFrames from the points yielding the min and max
min_gdf = gp.GeoDataFrame({'geometry':gp.GeoSeries(p)})
max_gdf = gp.GeoDataFrame({'geometry':gp.GeoSeries(points[max_idx])})
# Reset segment and carry on
last_cdp = this_cdp
points = [Point(x, y, twtt)]
return gp.GeoDataFrame({'geometry': linestrings,
'name': names,
'min': mins,
'max': maxs,
'minx': minx,
'miny': miny
})
Test on a small bit:
In [20]:
hrz = gdf_from_iesx(horizon_file)
In [21]:
hrz.head()
Out[21]:
Add a CRS...
In [23]:
from fiona.crs import from_epsg
hrz.crs = from_epsg(26720)
Try a quick plot...
In [26]:
fig = plt.figure(figsize=(18,10))
ax = hrz.plot(column='min')
plt.show()
Write a file if you like.
In [27]:
# hrz.to_file('../data/horizons_2d_UTM20N_NAD27.shp')
A format for importing wells into Petrel (or exporting from Petrel?). (Leaving this for later.)
# Petrel well head
VERSION 1
BEGIN HEADER
Well_ID
Well_sequence_number
Well_name
Surface X [m] NAD 83
Surface Y [m] NAD 83
KB elevation above sea level
TD [m]
Operator
END HEADER
' P-139 ' ' 139 ' ' ECE-13-P2_ST ' 529557 5045390 44.68 1,259.00 ' East Coast Energy '
' P-138 ' ' 138 ' ' ECE-13-P1 ' 529996 5045966 39.89 700 ' East Coast Energy '
' P-137 ' ' 137 ' ' Forent South Branch No.1 K-70-D/11-E-03 ' 496239 5003586 55.9 784 ' Forent Energy Ltd '
' P-136 ' ' 136 ' ' Forent Alton No. 1 E-49-C 11-E-03 ' 479105 5003441 40.23 1001 ' Forent Energy Ltd. '
' P-135 ' ' 135 ' ' E-38-A/11-E-5 ' 420980 5056480 123 946 ' Eastrock Resources '
' P-134 ' ' 134 ' ' ECE-11-01 ' 529976 5045865 43 678 ' East Coast Energy '
' P-133 ' ' 133 ' ' E-38-A/11-E-5 ' 443001 5015963 51 1726 ' Elmworth Energy Corporation '
In [ ]: