GEOG 485:
GIS Programming and Automation

4.4 Writing geometries

PrintPrint

As you parse out geographic information from "raw" sources such as text files, you may want to convert it to a format that is native to your GIS. This section of the lesson discusses how to write vector geometries to ArcGIS feature classes. We'll read through the same GPS-produced text file from the previous section, but this time we'll add the extra step of writing each coordinate to a polyline shapefile.

You've already had some experience writing point geometries when we learned about insert cursors. To review, you use arcpy.Point() to create a Point object and then use this object in the tuple given to insertRow() for the geometry field referred to as "SHAPE@" (see 4.1).

# Create point
inPoint = arcpy.Point(-121.34, 47.1)
...

# Create new row
cursor.insertRow((inPoint,))

For polylines and polygons, you create multiple Point objects that you add to an Array object. Then you make a Polyline or Polygon object using the array. With polygons it's a good practice to make the end vertex the same as the start vertex if possible.

The code below creates an empty array and adds three points using the Array.add() method. Then the array is used to create a Polyline object.

The first parameter you pass in when creating a polyline is the array containing the points for the polyline. The second parameter is a spatial reference of the coordinates, which you should always pass in to ensure that the precision of your data is maintained.

# Make a new empty array
array = arcpy.Array()

# Make some points
point1 = arcpy.Point(-121.34,47.1)
point2 = arcpy.Point(-121.29,47.32)
point3 = arcpy.Point(-121.31,47.02)

# Put the points in the array
array.add(point1)
array.add(point2)
array.add(point3)

# Make a polyline out of the now-complete array
polyline = arcpy.Polyline(array, spatialRef)

Of course, you usually won't create points manually in your code like this with hard-coded coordinates. It's more likely that you'll parse out the coordinates from a file or capture them from some external source, such as a series of mouse clicks on the screen.

Creating a polyline from a GPS track

Here's how you could parse out coordinates from a GPS-created text file like the one in the previous section of the lesson. This code reads all the points captured by the GPS and adds them to one long polyline. The polyline is then written to an empty, pre-existing polyline shapefile with a geographic coordinate system named tracklines.shp. If you didn't have a shapefile already on disk, you could use the Create Feature Class tool to create one with your script.

# This script reads a GPS track in CSV format and
#  writes geometries from the list of coordinate pairs
import csv
import arcpy

# Set up input and output variables for the script
gpsTrack = open("C:\\data\\Geog485\\gps_track.txt", "r")
polylineFC = "C:\\data\\Geog485\\tracklines.shp"
spatialRef = arcpy.Describe(polylineFC).spatialReference

# Set up CSV reader and process the header
csvReader = csv.reader(gpsTrack)
header = csvReader.next()
latIndex = header.index("lat")
lonIndex = header.index("long")

# Create an empty array object
vertexArray = arcpy.Array()

# Loop through the lines in the file and get each coordinate
for row in csvReader:
    lat = row[latIndex]
    lon = row[lonIndex]

    # Make a point from the coordinate and add it to the array
    vertex = arcpy.Point(lon,lat)
    vertexArray.add(vertex)

# Write the array to the feature class as a polyline feature
with arcpy.da.InsertCursor(polylineFC, ("SHAPE@",)) as cursor:
    polyline = arcpy.Polyline(vertexArray, spatialRef)
    cursor.insertRow((polyline,))   

The above script starts out the same as the one in the previous section of the lesson. First, it parses the header line of the file to determine the position of the latitude and longitude coordinates in each reading. But then, notice that an array is created to hold the points for the polyline:

vertexArray = arcpy.Array()

After that, a loop is initiated that reads each line and creates a point object from the latitude and longitude values. At the end of the loop, the point is added to the array.

for row in csvReader:
    lat = row[latIndex]
    lon = row[lonIndex]

    # Make a point from the coordinate and add it to the array
    vertex = arcpy.Point(lon,lat)
    vertexArray.add(vertex)

Once all the lines have been read, the loop exits and an insert cursor is created using "SHAPE@" as the only element in the tuple of affected fields. Then a Polyline object is created and used as the first (and only) element of the tuple given to insertRow():

# Create an insert cursor
with arcpy.da.InsertCursor(polylineFC, ("SHAPE@",)) as cursor:

    # Put the array in a polyline and write it to the feature class
    polyline = arcpy.Polyline(vertexArray, spatialRef)

    cursor.insertRow((polyline,))  

Remember that the cursor places a lock on your dataset, so this script doesn't create the cursor until absolutely necessary (in other words, after the loop).

Extending the example for multiple polylines

Just for fun, suppose your GPS allows you to mark the start and stop of different tracks. How would you handle this in the code? You can download this modified text file with multiple tracks if you want to try out the following example.

Notice that in the GPS text file, there is an entry new_seg:

type,ident,lat,long,y_proj,x_proj,new_seg,display,color,altitude,depth,temp,time,model,filename,ltime

new_seg is a boolean property that determines whether the reading begins a new track. If new_seg = true, you need to write the existing polyline to the shapefile and start creating a new one. Take a close look at this code example and notice how it differs from the previous one in order to handle multiple polylines:

# This script reads a GPS track in CSV format and
#  writes geometries from the list of coordinate pairs
#  Handles multiple polylines

# Function to add a polyline
def addPolyline(cursor, array, sr):
    polyline = arcpy.Polyline(array, sr)
    cursor.insertRow((polyline,))
    array.removeAll()

# Main script body
import csv
import arcpy

# Set up input and output variables for the script
gpsTrack = open("C:\\data\\Geog485\\gps_track_multiple.txt", "r")
polylineFC = "C:\\data\\Geog485\\tracklines_sept25.shp"
spatialRef = arcpy.Describe(polylineFC).spatialReference

# Set up CSV reader and process the header
csvReader = csv.reader(gpsTrack)
header = csvReader.next()
latIndex = header.index("lat")
lonIndex = header.index("long")
newIndex = header.index("new_seg")

# Write the array to the feature class as a polyline feature
with arcpy.da.InsertCursor(polylineFC, ("SHAPE@",)) as cursor:

    # Create an empty array object
    vertexArray = arcpy.Array()

    # Loop through the lines in the file and get each coordinate
    for row in csvReader:
        
        isNew = row[newIndex].upper()

        # If about to start a new line, add the completed line to the
        #  feature class
        if isNew == "TRUE":
            if vertexArray.count > 0:
                addPolyline(cursor, vertexArray, spatialRef)

        # Get the lat/lon values of the current GPS reading
        lat = row[latIndex]
        lon = row[lonIndex]

        # Make a point from the coordinate and add it to the array
        vertex = arcpy.Point(lon,lat)
        vertexArray.add(vertex)

    # Add the final polyline to the shapefile
    addPolyline(cursor, vertexArray, spatialRef)

The first thing you should notice is that this script uses a function. The addPolyline function adds a polyline to a feature class, given three parameters: (1) an existing insert cursor, (2) an array, and (3) a spatial reference. This function cuts down on repeated code and makes the script more readable.

Here's a look at the addPolyline function:

# Function to add a polyline
def addPolyline(cursor, array, sr):
    polyline = arcpy.Polyline(array, sr)
    cursor.insertRow((polyline,))
    array.removeAll()

Notice it's okay to use arcpy in the above function, since it is going inside the body of a script that imports arcpy. However, you want to avoid using variables in the function that are not defined within the function or passed in as parameters.

The addPolyline function is called twice in the script: once within the loop, which we would expect, and once at the end to make sure the final polyline is added to the shapefile. This is where writing a function cuts down on repeated code.

As you read each line of the text file, how do you determine whether it begins a new track? First of all, notice that we've added one more value to look for in this script:

  newIndex = header.index("new_seg")

The variable newTrackIndex shows us which position in the line is held by the boolean new_seg property that tells us whether a new polyline is beginning. If you have sharp eyes, you'll notice we check for this later in the code:

  isNew = row[newIndex].upper()

  # If about to start a new line, add the completed line to the
  #  feature class
  if isNew == "TRUE":

In the above code, the upper() method converts the string into all upper-case, so we don't have to worry about whether the line says "true," "True," or "TRUE." But there's another situation we have to handle: What about the first line of the file? This line should read "true," but we can't add the existing polyline to the file at that time, because there isn't one yet. Notice that a second check is performed to make sure there are more than zero points in the array before the array is written to the shapefile:

        
  if vertexArray.count > 0:
      addPolyline(cursor, vertexArray, spatialRef)

The above code checks to make sure there's at least one point in the array, then it calls the addPolyline function, passing in the cursor and the array.

Here's another question to consider: How did we know that the Array object has a count property that tells us how many items are in it? This comes from the ArcGIS Desktop Help topic describing the Array class. In this section of the help there are topics describing each class in arcpy, and you'll come here often if you work with ArcGIS geometries in Python.

In the above-linked Array topic, find the Properties table in this topic and notice that Array has a read-only count property. If we were working with a Python list, we could use len(vertexArray), but in our case vertexArray is an Array object that is native to the ArcGIS geoprocessing programming model. This means it is a specialized object designed by Esri, and you can only learn its methods and properties by examining the documentation. Bookmark these pages!

Summary

You can write geometries to ArcGIS feature classes using a combination of geometry objects included with ArcGIS. The common workflow is to create Point objects, which you add to an Array object. You can use the Array object plus a spatial reference to create Polyline and Polygon objects. You then use an insert cursor to assign the geometry in the array to the feature class's geometry field (usually called "shape").

You may be wondering how you might create a multi-part feature (such as the state of Hawaii containing multiple islands), or a polygon with a "hole" in it. There are special rules for ordering and nesting Points and Arrays to create these types of geometries. These are covered in the course textbook, which brings us to...

Readings

Read Zandbergen 8.1 - 8.6, which contains a good summary of how to read and write Esri geometries.