Assignment 1 – Vector Analysis

Assignment

This assignment is based on Lab 6 – Vector Operations from GUS 5062 – Fundamentals of GIS. That lab exercises uses vector operations that you should be familiar with such as buffering, intersection, erasing, etc., to identify an area in Philadelphia suitable for a new store selling healthy foods. Many of you will have taken GUS 5062 and completed Lab 6 in ArcGIS. Now you will see how to automate it using ArcPy. For those of you unfamiliar with the exercise, reading the assignment at the end will give you a little more context.

The criteria for the target area, and the data sources necessary for the analysis, are:

  1. In a Philadelphia Empowerment Zone. These are distressed urban areas targeted for reinvestment.
  2. Within 2000 feet of a SEPTA High Speed Station or SEPTA Regional Rail Station.
  3. NOT within 1200 feet of an establishment that already sells fresh produce, including Farmers Markets and stores participating in the Healthy Corner Stores Initiative.
  4. A larger contiguous area is preferred to a smaller area.

Preliminaries

Begin by downloading all five shapefiles and unzipping them into a project folder. Unzip them so that all the shapefiles are at the same folder level. Do not have each shapefile in its own subfolder, as that will make referring to the file in code more verbose.

Your Python script should be saved and run in the same folder as the unzipped data files. When you run the script using Spyder, Spyder should set the working directory to the script directory. If it doesn’t, you can use the “Run configuration per file” dialog box (CTRL+F6) in Spyder to set the working directory to “The directory of the file being executed”. To make the ArcGIS workspace match your Python working directory, your script should being with:

arcpy.env.workspace = os.getcwd()

With your workspace set correctly, each shapefile can be referred to in code using a bare filename rather than a fully qualified path, e.g. “my_shapefile.shp”, not “C:/folder1/folder2/folder3/my_shapefile.shp”.

Analysis Plan

  1. Begin by creating an inclusion layer. This should be the intersection of the areas that you want to include.
  2. Then create an exclusion layer. This is the area that you do NOT want to include.
  3. Use the exclusion layer to erase the non-desired areas from the inclusion layer.
  4. In order to calculate the area of each candidate zone, first use the arcpy.MultipartToSinglepart_management tool to explode the resulting multipolygon; then use the arcpy.AddGeometryAttributes_management tool to calculate the area.
  5. Find the area of the largest candidate zone. You can do this by opening the attribute table in ArcGIS. Put your answer in a comment at the end of your script. If you would like to do this step in code, you can try to do so by using a search cursor (PfA 17) or by using the arcpy.Statistics_analysis tool.

ArcPy Toolbox Functions That You Will Use

You will use the following ArcPy Toolbox functions in your analysis. Help for all these functions can be viewed within Spyder using CTRL+I, or by searching the documentation at ArcGIS.com. The names are specific enough that a Google search of, e.g., “ArcPy Buffer_analysis” will usually take you directly to the ArcGIS Pro or Desktop documentation for the tool. If the tool is discussed in the textbook, a page reference is given.

All of the following functions take as their first parameter a filename or list of filenames as inputs, and take as their second parameter an output filename that will get created on disk. They may have additional required or optional parameters.

  • arcpy.Project_management.
  • arcpy.Merge_management: PfA p. 101.
  • arcpy.Intersect_analysis: PfA p. 102.
  • arcpy.Buffer_analysis: First introduced in PfA on page 7 and used repeatedly for examples throughout the textbook.
  • arcpy.Dissolve_management.
  • arcpy.Erase_analysis.
  • arcpy.MultipartToSinglepart_management

Additionally you will need the arcpy.AddGeometryAttributes_managementtool. This function only takes an input file, which will be modified, as the first parameter (that is, no output file is specified). The required second parameter is the geometry attribute (or list of attributes) that you want to add. Use the following line of code to create a new field with the area of each feature. The new field will automatically be named POLY_AREA:

arcpy.AddGeometryAttributes_management("your_file.shp", "AREA")

Since the coordinate system units are feet, the area will be square feet unless you use specify a different unit.

Note that if you search the ArcPy documentation you will find that there are other ways of accomplishing the same thing, including CalculateGeometryAttributes_management and CalculateField_management. The tool I am suggesting is much easier to use, at the cost of being less flexible. For example, you have no control over the name of the created field, which is determined by the geometric measure that you calculate (area, length, centroid, point count, etc.)

Tips for Completing the Assignment

  • All of the data must be in the same projection. The downloaded shapefiles are all already in State Plane Pennsylvania South except the SEPTA Regional Rails Stations. You have to specify the output coordinate system when you call arcpy.Project_management, but keep in mind that you can pass the name of a .prj file to the out_coor_system parameter. All of the other shapefiles are in the desired coordinate system. Use the .prj from one of these files as input to this parameter. (Pay close attention to the filenames because the SEPTA High Speed Rail Stations are already in the correct coordinate system.)
  • In each case where you have to create a buffer, you need to buffer points from two different sources. It will be easiest if you merge the point layers first. Then you only have to run the buffer once. It will also be useful to use dissolve_option = "ALL" to dissolve the buffers.
  • Each toolbox function will create a new shapefile on disk. You can inspect these in ArcCatalog or add them to an ArcMap document to see your progress. If you add a layer to ArcMap, make sure to remove it before proceeding. ArcGIS file locking may prevent Python from being able to work with the file.

Deliverables

  • A script which completes all of these tasks.
  • The area of the largest contiguous zone that satisfies all criteria, either as a comment in the script, or printed to the console if you calculate it using ArcPy functions.
  • A map of the final layer of candidate zones (all candidate areas, not just the largest one). This map can be created in ArcGIS and exported to PDF. You do not have to make it pretty.

Grading

This assignment is worth 10 points, awarded 1 point for each of the following:

  • Create a script that runs to completion with no errors.
  • Set the ArcPy workspace correctly so that the script will execute successfully in a directory with the analysis data.
  • Use course coding conventions throughout.
  • Correctly reproject SEPTA Regional Rail Stations shapefile.
  • Correctly buffer Healthy Corner Stores and Farmers Markets shapefiles.
  • Correctly buffer SEPTA High Speed and Regional Rail Stations shapefiles.
  • Construct correct inclusion zone.
  • Correctly combine inclusion and exclusion zones.
  • Correctly add area field to singlepart result layer.
  • Provide adequate comments throughout your code.