Getting PETRA well bore paths into ArcMAP using Python

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
About these ads

2 Comments

Filed under Uncategorized

2 responses to “Getting PETRA well bore paths into ArcMAP using Python

  1. H. Alex Huskey

    Arne I just discovered your blog and appreciate your sharing! Maybe I misunderstood, but I export horizontal wellbore plots and sfc locs as shapefiles almost every day from Petra and use them in ArcGIS/INFO/MAP/VIEWER (or Arc for short) to make better looking maps than Petra can dream of. Just use the Petra’s Thematic Mapper to export the shapefiles of “deviated well path” and surface location to a destination folder. Then find them with ArcCatalog & drag or send them to your Arc project layer pane. You must project them once in Arc, but then customize line width, color, etc. At 600 point JPEG, the charts look “maptastic!”

    • Arne

      Sorry for never replying to this… like I said (or did I?), I’m not much of a Petra user, and while my approach may seem laborious, it makes sense to me. I remembered your response a few days when a co-worker was trying to do the same and used TM. I need to compare the results of both some time. Thanks for your comment.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s