GEOG 485:
GIS Programming and Software Development

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, if you put the X and Y coordinates in a tuple or list, you can plug it in the tuple given to insertRow() for the geometry field referred to using the "SHAPE@XY" token (see page 4.1).

# Create coordinate tuple
inPoint = (-121.34, 47.1)
...

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

At ArcGIS Desktop v10.6/ArcGIS Pro v2.1, Esri made it possible to create polylines and polygons by putting together a list or tuple of coordinate pairs (vertices) like the one above in sequence.  When you pass that list or tuple to the insertRow() method, arcpy will "connect the dots" to create a polyline or a polygon (depending on the geometry type of the feature class you opened the insert cursor on).  Multi-part and multi-ring geometries are a bit more complicated than that, but that's the basic idea for single-part geometries.  

The code below creates an empty list and adds three points using the list.append() method. Then that list is plugged into an insertRow() statement, where it will result in the creation of a Polyline object.

# Make a new empty list
coords = []

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

# Put the points in the list
coords.append(point1)
coords.append(point2)
coords.append(point3)

# Open an insert cursor on the FC, add a new feature from the coords list
with arcpy.da.InsertCursor(polylineFC, ("SHAPE@")) as cursor:
    cursor.insertRow((coords))

In addition to the requirement that the geometry be single-part, you should also note that this list-of-coordinate-pairs approach to geometry creation requires that the spatial reference of your coordinate data matches the spatial reference of the feature class.  If the coordinates are in a different spatial reference, you can still create the geometry, but you'll need to use the alternate approach covered at the bottom of this page.  

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

polylineFC = r"C:\PSU\geog485\L4\trackLines.shp"
 
# Open the input file 
with open(r"C:\PSU\geog485\L4\gps_track.txt", "r") as gpsTrack:
    # Set up CSV reader and process the header
    csvReader = csv.reader(gpsTrack)
    header = next(csvReader)
    latIndex = header.index("lat")  
    lonIndex = header.index("long")
 
    # Create an empty list
    vertices = []
 
    # Loop through the lines in the file and get each coordinate
    for row in csvReader:
        lat = float(row[latIndex])
        lon = float(row[lonIndex])
 
        # Put the coords into a tuple and add it to the list
        vertex = (lon,lat)
        vertices.append(vertex)
 
    # Write the coordinate list to the feature class as a polyline feature
    with arcpy.da.InsertCursor(polylineFC, ('SHAPE@')) as cursor:
        cursor.insertRow((vertices,))

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.  After that, a loop is initiated that reads each line and creates a tuple containing the longitude and latitude values. At the end of the loop, the tuple is added to the list.

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 the insertRow() method is called, passing it the list of coordinate tuples within a tuple.  It's very important to note that this statement is cursor.insertRow((vertices,)), not cursor.insertRow(vertices,).  Just as the fields supplied when opening the cursor must be in the form of a tuple, even if it's only one field, the values in the insertRow() statement must be a tuple.

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).  Finally, note that the tuple plugged into the insertRow() statement includes a trailing comma.  This odd syntax is needed only in the case where the tuple contains just a single item.  For tuples containing two or more items, the trailing comma is not needed.  Alternatively, the coordinate list could be supplied as a list within a list ([vertices]) rather than within a tuple, in which case a trailing comma is also not needed.

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, coords):
    cursor.insertRow((coords,))
    del coords[:] # clear list for next segment

# Main script body
import csv
import arcpy

polylineFC = "C:\\data\\Geog485\\tracklines.shp" 

# Open the input file 
with open(r"C:\PSU\geog485\L4\gps_track_multiple.txt", "r") as gpsTrack:
    # Set up CSV reader and process the header
    csvReader = csv.reader(gpsTrack)
    header = next(csvReader)
    latIndex = header.index("lat")  
    lonIndex = header.index("long")
    newIndex = header.index("new_seg")

    # Write the coordinates to the feature class as a polyline feature
    with arcpy.da.InsertCursor(polylineFC, ("SHAPE@")) as cursor:
        # Create an empty vertex list
        vertices = []

        # 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 len(vertices) > 0:
                    addPolyline(cursor, vertices)
            # Get the lat/lon values of the current GPS reading
            lat = float(row[latIndex])
            lon = float(row[lonIndex])

            # Add coordinate pair tuple to vertex list
            vertices.append((lon, lat))

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

The first thing you should notice is that this script uses a function. The addPolyline() function adds a polyline to a feature class, given two parameters: (1) an existing insert cursor, and (2) a list of coordinate pairs. 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, coords):
    cursor.insertRow((coords,))
    del coords[:]

The addPolyline function is referred to 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 newIndex 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 list before attempting to add a new polyline:

        
  if len(vertices) > 0:
      addPolyline(cursor, vertices)

Only if there's at least one point in the list does the addPolyline() function get called, passing in the cursor and the list.

Alternate method for creating Polylines and Polygons

Prior to ArcGIS Desktop v10.6/ArcGIS Pro v2.1, the list-of-coordinate-pairs approach to creating geometries described above was not available.  The only way to create Polylines and Polygons was to create an arcpy Point object from each set of coordinates and add that Point to an Array object. Then a Polyline or Polygon object could be constructed from the Array. 

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.

# 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)

The first parameter you pass in when creating a polyline is the array containing the points for the polyline. The second parameter is the spatial reference of the coordinates.  Recall that we didn't have any spatial reference objects in our earlier list-of-coordinate-pairs examples.  That's because you can only use that method when the spatial reference of the coordinates is the same as the feature class.  But when creating a Polyline (or Polygon) object from an Array, you have the option of specifying the spatial reference of the coordinates.  If that spatial reference doesn't match that of the feature class, then arcpy will re-project the geometry into the feature class spatial reference.  If the two spatial references are the same, then no re-projection is needed.  It can't hurt to include the spatial reference, so it's not a bad idea to get in the habit of including it if you find yourself creating geometries with this alternate syntax.  

Here is a version of the GPS track script that uses the Array-of-Points approach:

# 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

polylineFC = "C:\\data\\Geog485\\tracklines_sept25.shp" 
spatialRef = arcpy.Describe(polylineFC).spatialReference

# Open the input file 
with open("C:\\data\\Geog485\\gps_track_multiple.txt", "r") as gpsTrack:
    # Set up CSV reader and process the header
    csvReader = csv.reader(gpsTrack)
    header = next(csvReader)
    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 = float(row[latIndex])
            lon = float(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)

Summary

The ArcGIS Pro documentation does a nice job of summarizing the benefits and limitations of creating geometries from lists of coordinates:

Geometry can also be created from a list of coordinates. This approach can provide performance gains, as it avoids the overhead of creating geometry objects. However, it is limited to only features that are singlepart, and in the case of polygons, without interior rings. All coordinates should be in the units of the feature class's spatial reference.

If you need to create a multi-part feature (such as the state of Hawaii containing multiple islands), or a polygon with a "hole" in it, then you'll need to work with Point and Array objects as described in the Alternate method section of this page.  You would also use this method if your coordinates are in a different spatial reference than the feature class.

Readings

Read the Writing geometries page in the ArcGIS Pro documentation, and pay particular attention to the multipart polygon example if you deal with these sorts of geometries in your work.

ArcGIS Pro edition:

Read Zandbergen 9.1 - 9.7, which contains a good summary of how to read and write Esri geometries using the Points-in-an-Array method.

ArcMap edition:

Read Zandbergen 8.1 - 8.6, which contains a good summary of how to read and write Esri geometries using the Points-in-an-Array method.