As
I had previously mentioned in my proposal and presentation, the script I set
out to create was designed mainly for archaeology contract work. Basically, it
would allow the user to obtain quick geometry calculations pertaining to the
amount of various disturbances within an archaeological study area (or APE).
The script would work in several distinct steps, first by creating a clip
polygon of all input polygons within the given APE, create an area column,
calculate the acreage of the overlap for each intersection (while at the same
time creating an acreage of an entire APE), printing the results, creating
another column in each of the clipped polygons (and writing the total Study Area
acreage within it) as well as a percentage of coverage column, performing a
division calculation to determine the percentage of the APE that each clip
occupies and finally printing the second set of results. If all went smoothly
in PythonWin, my intention was to revisit the later chapters and create a
toolbox which would allow the user to use the script as a tool with
manually-defined parameters.
The
first step of operation was actually done outside of Python. I created a new
MXD, imported a basemap and zoomed in on a random area in Pennsylvania. After
setting the coordinate system of the MXD to ‘State Plane – Pennsylvania South
US feet’ I began creating hypothetical survey data which I would use to test my
script. I drew an APE polygon around an intersection, covering the area of
several residential yards. Next I created three classes of disturbance which
stretched outside of/and intersected with my study area. I called them ‘Wetland.shp’,
‘Excessive Slope.shp’, and ‘Modern Disturbance.shp’. (All shapefiles above
mentioned were polygons. The Modern Disturbance file covered the asphalt road
and well as all driveways and houses, the Wetland file contained two separate
parts and the slope was made up of just one large area to the south. Once the
shapefiles were finished, I exited arcmap. My intention was to turn it on and
off while testing my script, in order to avoid errors when making edits to the
same file on two different programs. The script I am creating is, as of now,
limited to three separate disturbance inputs. This was done simply to conserve
time in the beginning, but time-permitting I would be able to go back and add
more if necessary.
I
began the script by simply importing all modules I thought may be necessary to
run it. I imported arcpy, os, fileinput and string, and from arcpy I imported
env. I set the workspace to my shapefile folder and enabled overwrite output. Within
the same part of the script, I decided to define all of my most common and most
necessary inputs. I defined my mxd file with a filepath, and the study area was
defined shortened to “clip_features”. I also defined “clip_feature1”,
“clip_feature2”, and “clip_feature3” as my disturbance shapefiles. Because my
first operation would be to clip the three, I defined the path of the initial
clipped polygons and shortened them to “clipped_polygon_1”, “clipped_polygon2”,
and “clipped_polygon3”. I thought that the best way to troubleshoot the script,
would be to print messages at the beginning of each portion of the script, as
well as printing “complete” at the end of each operation. That way, whenever I
would run the code I would know exactly how much of it worked flawlessly and
what part needed fixing.
The
first actual operation performed by the script, would be to clip all of the
disturbance shaefiles in order to isolate the overlap within the APE. I used
the arcpy.Clip_analysis function and named each input feature to be clipped by
the clip_feature. I did this three individual times for each of the polygons,
but added an “if” clause for input 2 and 3. The ‘if’ condition was simply a
presence of the given shapefile as it is defined. If it did not exist, or was
not defined the script would print a blank space and move on. The idea was to
allow the script to be run with less than three disturbance areas whenever
necessary without showing any errors. (I applied the same formula in all of the
following steps, when it came to performing operations on the disturbance
polygons.) After the clip operation was completed, I checked my output folder in
ArcCatalog to make sure that the new shapefiles were contained there. They
were.
Next,
I had to add my area columns to the APE and the clipped disturbances. I began
by defining local variables to make my work easier, and I defined fieldName1 as
“Area” and fieldName2 as “APE_Area”. I added the “APE_Area” column to my APE
shapefile using the arcpy.AddField_management function. The column created was
a Double number with precision set to 15 and field length to 10. This was
necessary when dealing with a decimal point, and because the areas we calculate
at work are always in Acres, a decimal value makes a big difference. I repeated
the step individually (as in the clipping step above) for my disturbance
inputs, but instead of the “APE_Area” column I created the “Area” column. I
made this distinction, because I would have to merge the two shapefiles later
and could not have the same column name twice with different values. My logic
was that this would make my life easier later, when calculating the percentages
of the area is question in a later step. I ran my script to test it, and opened
my mxd to see if the changes had been made. Thus far, it ran smoothly and all
the necessary and the APE file as well as all the initial clips contained a new
column within their attribute table.
Within
the same section of the script, I went back to perform area calculations. As I
had stated before, at work we always deal with acres. I used the
aprcpy.CalculateField_management() function on each of the polygons, populating
the “APE_Area” column in the APE shapefile with its acreage, and the “Area”
column for each of the disturbance polygons. I made sure that the operation is
being done on the clipped polygons instead of the originals, and added each of
the operations directly below the corresponding arcpy.AddField_management
operations from the previous paragraph. I saved my script and ran it, then
re-opened my mxd and went to check my results. Thus far, all seemed to be in
place and the columns were populated by the acreage total with a decimal.
The
fourth step of the script was to simply print the area results. This would be
the first part of the code that would allow the user to obtain new information
about the data they are working with. In order to display the results in a user-friendly
manner I used the search cursor function. For the APE, I defined the cursor to
the “APE_Area” field, and printed a message joined by two string texts: "The
APE has a total area of", (row), "acres." This way, python would
print out the acreage of each separate polygon feature within a shapefile in
question. For the next three input files, I redefined the cursor to focus on
the “Area” field instead, in order to obtain the acreage of the clipped
polygons. Because the names, and the order of disturbances would vary depending
on project the message I printed was: "The first disturbance listed
includes an area the size of", (row), "acres within the APE." I
saved and ran my script, and all was well thus far. I was happy to see that
both rows of the first disturbance were listed separately and contained
correct, individual acreage.
Next
came the time to import the “APE_Area” field into my disturbance polygons. I
had begun running into a bit of trouble here, as I had initially hoped that I
would be able to perform these calculations by calling values from various
polygons. For example, I wanted to input a simple equation where the value of
the “Area” field for each of the polygons could be divided by the value of the
“APE_Area” and multiplied by 100. I intended to use the result of this equation
to populate a new field in the attribute table of the disturbance polygons in
order to display the percentage of APE that is taken up by each of these
features. But after a very long time on Esri forums, and many attempts at
troubleshooting I saw that I was not able to figure out how to do this
properly. It was then that I decided to use the Spatial Join function to create
yet another set of shapefiles, which would contain both columns in their
Attribute table (meaning: each shapefile would contain as “Area” field that
would display the acreage of its own features, as well as an “APE_Area” field
that would display the number representing the acreage of the entire APE. Before
doing this, however I defined a few local attributes once more and created “PerCofAPE”.
(This field was meant to display the result of the percentage equation I had
mentioned previously.)
I
executed the arcpy.SpatialJoin_analysis() function on all three of my clipped
disturbance polygons (combining them with my APE polygon) and called them "OFC1_join.shp",
"OFC2_join.shp" and "OFC3_join.shp" respectively. I ran
another arcpy.AddField_management() operation to add the “PerCofAPE” column to
each of the joined polygons as a Double, saved and ran my script. When I opened
up ArcMap, I saw that the created polygons were in their correct place and
after viewing their attributes I saw that they contained individual columns
with each of the desired numbers as well as the new, blank “PerCofAPE” field.
The last (and most difficult) operation in my script was to populate the “PerCofAPE” field with the result of the percentage equation. (As I had mentioned above, this would be (“Area”)/(“APE_Area”)*100). I figured that the easiest way to do this was to simply perform a field calculation on all the joined shapes using the arcpy.CalculateField_management function. It looked as though the equation would yield a result (in a form of a decimal, for now) but no matter what I did I would receive an error message stating that the field in question was not nullable. This was extremely confusing, as I had made sure that all created fields throughout the script were designated as nullable. I spent a long time trying to troubleshoot the problem, but ultimately I was not able to figure it out on time. The last remaining step in the script, was supposed be to simply print these percentages in python in a legible manner. Unfortunately, no matter what I had tried the “PerCofAPE” field remained blank and not editable.
Ultimately, this script proved to be a lot more challenging and time consuming. I had previously thought that simple calculations would be easy to figure out and troubleshoot but in this particular case it proved to be too much for me. I am proud that I was able to achieve most of what I set out to do, and if I were able to get past this particular issue, I think the script would have worked well. I had also toyed with the idea of setting the script up as a user-friendly tool and get parameters that could be set in ArcMap (like we did in our later modules) directly. But unfortunately this was impossible as I had gotten stuck on a crucial part of the script and could not move on to the final steps.
No comments:
Post a Comment