GEOG 485:
GIS Programming and Software Development

Error message

Notice: unserialize(): Error at offset 14 of 14 bytes in variable_initialize() (line 1255 of /net/www/www.e-education.psu.edu/htdocs/drupal7/includes/bootstrap.inc).

Lesson 4 Practice Exercise A Solution

PrintPrint

Here's one way you could approach Lesson 4 Practice Exercise A with comments to explain what is going on.  If you find a more efficient way to code a solution, please share it through the discussion forums.

# Reads coordinates from a text file and writes a polygon

import arcpy
import csv

arcpy.env.workspace = r"C:\PSU\geog485\L4\Lesson4PracticeExerciseA"

shapefile = "MysteryState.shp"
pointFilePath = arcpy.env.workspace + r"\MysteryStatePoints.txt"

# Open the file and read the first (and only) line
with open(pointFilePath, "r") as pointFile:
    csvReader = csv.reader(pointFile)

    # This list will hold a clockwise "ring" of coordinate pairs
    #  that will form a polygon 
    ptList = []

    # Loop through each coordinate pair 
    for coords in csvReader:
        # Append coords to list as a tuple
        ptList.append((float(coords[0]),float(coords[1])))

    # Create an insert cursor and apply the Polygon to a new row
    with arcpy.da.InsertCursor(shapefile, ("SHAPE@")) as cursor:
        cursor.insertRow((ptList,))

Alternatively, an arcpy Array containing a sequence of Point objects could be passed to the insertRow() method rather than a list of coordinate pairs.  Here is how that approach might look:

# Reads coordinates from a text file and writes a polygon

import arcpy
import csv

arcpy.env.workspace = r"C:\PSU\geog485\L4\Lesson4PracticeExerciseA"

shapefile = "MysteryState.shp"
pointFilePath = arcpy.env.workspace + r"\MysteryStatePoints.txt"
spatialRef = arcpy.Describe(shapefile).spatialReference

# Open the file 
with open(pointFilePath, "r") as pointFile:
    csvReader = csv.reader(pointFile)

    # This Array object will hold a clockwise "ring" of Point
    #  objects, thereby making a polygon.
    polygonArray = arcpy.Array()

    # Loop through each coordinate pair and make a Point object
    for coords in csvReader:

        # Create a point, assigning the X and Y values from your list    
        currentPoint = arcpy.Point(float(coords[0]),float(coords[1]))

        # Add the newly-created Point to your Array    
        polygonArray.add(currentPoint)

    # Create a Polygon from your Array
    polygon = arcpy.Polygon(polygonArray, spatialRef)

    # Create an insert cursor and apply the Polygon to a new row
    with arcpy.da.InsertCursor(shapefile, ("SHAPE@")) as cursor:
        cursor.insertRow((polygon,))


Below is a video offering some line-by-line commentary on the structure of these solutions.  Note that the scripts shown in the video don't use the float() function as shown in the solutions above.  It's likely you can run the script successfully without float(), but we've found that on some systems the coordinates are incorrectly treated as strings unless explicitly cast as floats.  

Video: Solution to Lesson 4 Practice Exercise A (11:02)

Click here for transcript of the Solution to Lesson 4, Practice Exercise A.

This video walks through the solution to Lesson 4, practice exercise A where you have a text file with some coordinates that looks like this. And you need to loop through them, parse them out, and then create a polygon out of it to create the shape of a US state.

This is a very basic list of coordinates. There's no header here. So, we're just going to blast right through this as if it were a csv file with two columns and no header.

In order to do this, we're going to use arcpy, which we import in line 3, and the csv module for Python, which we import in line 4. In lines 8 and 9, we set up some variables referencing the files that we're going to work with. So, we're working with the mystery state shapefile that was provided for you.

This shapefile is in the WGS84 geographic coordinate system. However, there is no shape or polygon within this shapefile yet. You will add that. And then, line 9 references the path to the text file containing all the coordinate points.

In line 12, we open up that text file that has the points. The first parameter in the open method is pointFilePath, which we created in line 9. And then the second parameter is r in quotes. And that means we're opening it in read mode. We're just going to read this file. We're not going to write anything into it.

And in line 13, we then pass that file into the reader() method so that we can create this csvReader object that will allow us to iterate through all the rows in the file in an easy fashion. Before we start doing that though, in line 17, we create an empty list which we’ll use to store the coordinate pairs we read in from the text file.

We know, in this case, we're just going to create one geometry. It's going to be one polygon created from one list of coordinates. So, it's OK to create that list before we start looping. And we, actually, don't want to create it inside the loop because we would be creating multiple lists, and we just need one of them in this case.

In line 20, we actually start looping through the rows in the file. Each row is going to be represented here by a variable called coords. Inside the loop, we’re going to get the X coordinate and the Y coordinate and append them to the list as a tuple. It’s very important that the first value in this tuple be the X and the second value the Y.

Now, I think a lot of us typically say “latitude/longitude” when talking about coordinates, with the word latitude coming first. It’s important to remember that latitude is the Y value and longitude is the X. So in defining this tuple, we want to make sure we’re supplying the longitude first, then the latitude.

Looking at our text file, we see that the first column is the longitude and the second column is the latitude. (We know that because latitudes don’t go over 90.) So the coordinates are in the correct X/Y order, we just need to keep them that way.

Back to our script, as we iterate through the CSV Reader object, remember that on each iteration, the Reader object simply gives us a list of the values from the current line of the file, which we’re storing in the coords variable. So if we want the longitude value, that’ll be the item from the list at position 0. If we want the latitude, that’ll be the item at position 1. So, we use coords[0] to get the longitude and coords[1] to get the latitude.

Now, going back to our file, you'll see that in the first column or column 0, that's the longitude. So, we have coords[0] in our code. And then, the next piece is the latitude. So, we use coords[1] for that.

So we’re only doing one thing in this loop: appending a new X/Y coordinate pair to the list on each iteration.

Now, line 25 is the first statement after the loop. We know that because that line of code is not indented. At this point, we’ve got the coordinates of the mystery state in a list, so we’re ready to add a polygon to the shapefile. To add new features to a feature class, we need to use an InsertCursor. And that’s what we’re doing there on line 25: opening up a new insert cursor.

That cursor needs a couple of parameters in order to get started. It needs the feature class that we're going to modify, which is referenced by the variable shapefile. Remember, we created that back in line 8.

And then, the second parameter we need to supply a tuple of fields that we're going to modify. In our case, we're only going to modify the geometry. We're not going to modify any attributes, so the way that we supply this tuple is it just has one item, or token, SHAPE with the @ sign.

And then, our final line 26 actually inserts the row. And it inserts just the list of points we populated in our loop. Because the shapefile was created to hold polygons, a polygon geometry will be created from our list of points. When this is done, we should be able to go into ArcGIS Pro and verify that we indeed have a polygon inside of that shapefile.

One little quirky thing to note about this code is in the last line, where I supplied the tuple of values corresponding to the tuple of fields that was associated with the cursor. When you’re supplying just a single value as we were here, it’s necessary to follow that value up with a comma. If my cursor was defined with more than one field, say 2 fields, then my values tuple would not need to end in a comma like this. It’s only when specifying just a single value that the values tuple needs this trailing comma.

This is a good place to remind you that you’re not limited to setting the value of a single field. For example, if there were a Name field in this shapefile, and we happened to know the name of the state we were adding, we could add Name to the fields tuple and then when doing insertRow() we could supply the name of the state we’re adding in addition to the point list.

Now, the script we just walked through uses the approach to geometry creation that’s recommended to use whenever possible because it offers the best performance. That approach works when you’re dealing with simple, single-part features and the coordinates are in the same spatial reference as the feature class you’re putting the geometry into. And when you’re working with ArcGIS Desktop 10.6 or later or ArcGIS Pro 2.1 or later because that’s when support for this approach was added.

So for example, this approach would *not* work for adding the state of Hawaii as a single feature since that would require inserting a multi-part polygon. It also wouldn’t work if our empty shapefile was in, say, the State Plane Coordinate System rather than WGS84 geographic coordinates.

In those cases, you’d have no choice but to use this alternate solution that I'm going to talk through now.

In this version of the script, note that in line 10, we get the spatial reference of the MysteryState shapefile using the Describe method. We’ll go on to use this SpatialReference object later in the script when we create the Polygon object. We didn’t bother to get the spatial reference of the shapefile in the other version of the script because that method of creating geometries doesn’t support supplying a spatial reference.

The next difference in this second script comes on line 18. Instead of creating a new empty Python list, we create an arcpy Array object.

Then inside the loop, instead of putting the coordinates into a tuple and appending the tuple to the list, we create an arcpy Point object, again in X/Y order. Then, we use the Array’s add() method to add the Point object to the Array.

After the loop is complete and we’ve populated the Array, we’re ready to create an object of the Polygon class. The Polygon constructor requires that we specify an Array of Points, which is in our polygonArray variable. Optionally, we can supply a SpatialReference object, which we’re doing here with our spatialRef variable.

Supplying that spatial reference doesn’t really have any effect in this instance because our coordinates and feature class are in the same spatial reference, but you could see the difference between these two if you created an empty feature class that was defined to use a different spatial reference, say one of the State Plane coordinate systems for New Mexico. Running this version of the script would result in the data being properly aligned because we’re telling arcpy how our data is projected, so it can re-project it behind the scenes into the State Plane coordinates. With the other version of the script, the data can't be properly aligned because that version doesn’t specify a coordinate system.

In any case, this version of the script finishes up by passing the Polygon object to the insertRow() method rather than a list of coordinate pairs.

So, those are two ways you can solve this exercise.

Credit: S. Quinn and J. Detwiler © Penn State is licensed under CC BY-NC-SA 4.0.