I recently had to bring in some IHS well data from Petra (or IHS Enerdeq) into ArcMap. While I endured 2 days of Introduction to Petra for Geoscientists this summer, I easy and expertise when it comes to using Petra is rather limited. I agree with many online posters who complain about how clumsy, cluttered, and counter-intuitive PETRA is. But that’s another story. I was told by a well seasoned PETRA user that PETRA does not export to SHP file format, upon which she exported a well data table with lat/long’s for me as XLS. Whether the SHP limitation is indeed the case, I have not further investigated. I simply created an ESRI feature class from the LL data. Another needed step was exporting the well bore paths for directional wells into a polyline feature class for use in ArcMap. Again, the only way to get the data out of PETRA appeared to be as a table. I chose CSV.
The table has the following headers:
UWI | NAME | MD | TVD | EWOFFSET | NSOFFSET | DIP | AZM | TVDSS | XPATH | YPATH | SURNAME | SURSTAT
and provides a row of data for each survey depth, with UWI(API) repeating for all rows belonging to the same well.
To read this into Python, I did this:
import csv
### Reference to Input CSV File
infile = open('c:/temp/dirEXPORT.csv', "rb")
reader = csv.reader(infile)
rownum = 0
well = ""
Wells = []
points = []
### Read CSV File, Line by Line
for row in reader:
if rownum == 0: # First line in file: read headers
colcount = 0
for column in row:
print column,
if column == "UWI":
uwi = colcount
if column == "XPATH":
xpath = colcount
if column == "YPATH":
ypath = colcount
if column == "MD":
z = colcount
colcount = colcount + 1
rownum=+1
else:
print "\n"
### Check if row (line) is for the same well as row above
### if so, read coordinates from field and assign to xyz as list
if row[uwi] == well:
xyz = [row[xpath],row[ypath],row[z]]
points.append(xyz)
print row[uwi], xyz
rownum=+1
### If row does not belong to the same well as the the last row,
### then append coordinate list and well name to Wells list, reset
### coordinates list (=[]), and assign the coordinates to the new well
else:
Wells.append([well,points])
print "New Well"
points = []
well = row[uwi]
xyz = [row[xpath],row[ypath],row[z]]
points.append(xyz)
print row[uwi], xyz
rownum=+1
infile.close
Any remaining print statements allow me to watch what’s going on when the program is run. It’s my crude way of debugging code.
Once all the data is loaded into a list, I start some arcpy magic:
import arcpy
from arcpy import env
import os
env.overwriteOutput = True
env.workspace = "C:/temp"
cnt = 0
mypath = env.workspace
outputFile = r"wormtracks.shp"
template = r"c:\temp\template.shp"
fieldLength = 25
sr = arcpy.Describe(template).spatialReference
arcpy.CreateFeatureclass_management(mypath,outputFile,"POLYLINE","","","",sr)
arcpy.AddField_management(outputFile, 'API', "TEXT", "", "", fieldLength) # add API as new table field
rows = arcpy.InsertCursor(outputFile)
array_container = arcpy.Array()
for well in Wells[1:]: # Ignore first item in Wells list. Poor programming where an empty list is created.
well_api = well[0] # First item in list is the API
wellList = well[1] # 2nd item in list is a list of XYZ points
for pt in wellList:
point_object = arcpy.Point()
point_object.X = float(pt[0])
point_object.Y = float(pt[1])
print pt[0]
print pt[1]
array_container.append(point_object) # I couldn't figure out the difference between array.append and array.add
print point_object
del point_object
row = rows.newRow()
row.shape = arcpy.Polyline(array_container) # Do not use arcpy.Multipoint which created a really messy polyline w/points ouf of sequence.
row.API = well_api
rows.insertRow(row)
array_container.removeAll()
del row #unlock row
del rows #unlock table
del arcpy
del outputFile
Working in Python, I find myself googling for solutions to programming snags all the time. More and more often, my questions turn out to be about use of the Python Standard Library, be it hashlib for file comparison purposes or shutil for standard file system operations, or whatever. So I finally took the plunge and placed the Standard Library on my desk – the Addison-Wesley Professional/Developers Library text. I’ve said before, I learn more from reading 2 pages in a book than googling 20 pages.
but the format of the class (30 minutes teaching, 10 minutes hands-on), plus the quality of the manual, made for a real good package that I would recommend to anyone interested in FME. Well done ! — That said, I think FME is fantastic. I think I mentioned to both of you that I have pulled my hair out trying to wrap my brain around ESRI’s ModelBuilder more than once. FME makes more sense to me, is a lot more intuitive to my thinking and the work flows I deal with …
