Wednesday, July 13, 2011

Creating Extent Polygons Using ArcPy

Awhile ago I wrote an article about creating extent polygons using the old arcgisscripting geoprocessing framework at version 931.  I've decided to update the posting because ArcPy makes it easy because you can easily get the extent object of a geometry and directly grab the extent point objects from the extent object.

To create an extent of the whole dataset use the describe method:

import arcpy
from arcpy import env
import os

inFC = r"c:\temp\USA.gdb\States"
extentPoly = str(env.scratchWorkspace) + os.sep + "extent.shp"
# Feature extent
#
extent = arcpy.Describe(inFC).extent
# Array to hold points
#
array = arcpy.Array()
# Create the bounding box
array.add(extent.lowerLeft)
array.add(extent.lowerRight)
array.add(extent.upperRight)
array.add(extent.upperLeft)
# ensure the polygon is closed
array.add(extent.lowerLeft)
# Create the polygon object
polygon = arcpy.Polygon(array)
array.removeAll()
# save to disk
arcpy.CopyFeatures_management(polygon, extentPoly)
del polygon

That should create a box around all your data in a feature class. Very simple!

Next, let's say you need the bounding box for each feature. In the example, I will continue with using a states polygon:

import arcpy
from arcpy import env
import os


env.overwriteOutput = True

inFC = r"c:\temp\USA.gdb\States"
extentPoly = str(env.scratchWorkspace) + os.sep + "extent.shp"
ptList = []
rows = arcpy.SearchCursor(inFC)

for row in rows:
   extent = row.Shape.extent
   array.add(extent.lowerLeft)
   array.add(extent.lowerRight)
   array.add(extent.upperRight)
   array.add(extent.upperLeft)
   array.add(extent.lowerLeft)
   polygon = arcpy.Polygon(array)
   ptList.append(polygon)
   array.removeAll()
del rows
arcpy.CopyFeatures_management(ptList, extentPoly)

In this code sample, it loops through the feature class and write the bounding box geometry to a python list. That list is then given as an input to CopyFeatures() and it writes it to disk. Notice that you can get the extent right from the Shape value from the row object. This makes life easy.

Here is a image of what the lower 48 would look like as extents:
The white lines are the extents.

Enjoy

14 comments:

Anonymous said...

This code leaves the resulting file with no projection. How can I add it? I have tried setting the psatialrefrence property to the same as the file it came from but no luck.. I will not know what the prj it ahead of time.

Andrew said...

You can add a spatial reference to your polygon object. See http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v000000n1000000 for all the details.

To create a spatial reference say for WGS-1984 do the following:
sr = arcpy.SpatialReference()
sr.factoryCode = 4326
sr.create()
then you would pass the spatial reference along with your coordinates to the polygon object.

polygon = arcpy.Polygon(array,sr)

Enjoy

Jon said...

I read your heading and got all excited that my searching had come to an end. But it was a teaser.

I'm trying to find an arcpy script that does a data export to the data frame extents. I don't think it exists.

What I want my end result to do, is clip all (or selected) layer's features in the current mxd to the data frame extents and put them in a gdb with the same name (or maybe a '_clipped' suffix. I'm new to python scripting and have scoured the forums looking for an answer to what i thought would be simple.

All help is appreciated. Thanks.

Jon said...

Ok, so since I posted yesterday, i've been putting in some work. My script is now creating a .shp based on data frame extents, but it's not clipping layers in my data frame. Any ideas?

import arcpy
import arcpy.mapping
import os
# Set scratch workspace
arcpy.env.scratchWorkspace = "C:/Users/jsommerville/Documents/ArcGIS/scratch"
scratchpath = arcpy.env.scratchWorkspace

#Set the current map
mxd = arcpy.mapping.MapDocument("CURRENT")

# Set the dataframe
df = arcpy.mapping.ListDataFrames(mxd)[0]

# Get list of Layers from user imput
layers = arcpy.mapping.ListLayers(mxd,"",df)

#Set the output workspace
outWorkspace = r"C:\Users\jsommerville\Documents\ArcGIS\California\NewFolder"


# The following code uses the extents of the first data frame in the TOC of current map
# to create a polygon shapefile.

# Set the dataframe extents
dfx = df.extent

# Name and path for extents polygon
extentPoly = str(scratchpath) + os.sep + "extent.shp"

# Array to hold points
array = arcpy.Array()

# Create the bounding box
array.add(dfx.lowerLeft)
array.add(dfx.lowerRight)
array.add(dfx.upperRight)
array.add(dfx.upperLeft)
array.add(dfx.lowerLeft)

# Create the polygon object
polygon = arcpy.Polygon(array)
array.removeAll()

# save to scratch folder
arcpy.CopyFeatures_management(polygon, extentPoly)
del polygon


for lyr in layers:
outFeatureClass = os.path.join(outWorkspace, lyr)
arcpy.Clip_analysis(lyr, extentPoly, outFeatureClass)

Jake said...

I believe you're missing

array = arcpy.Array()

from the second example.

Thanks though, saved me writing it all out myself.

Anonymous said...

Thanks so much! I got a dbf of extents for a grid and could not figure out how to get it into a geometry file. I modified yours...


import arcpy, os
from arcpy import env

env.overwirteOutput = True

inTbl=r"dbf file path here of records of extents"

outShp=r"already made empty shapefile path here"

#Array to hold points.
array=arcpy.Array()

#Point List to generate polygon
ptList= []

#For table use SearchCursor, for another GIS file use da.SearchCurstor
rows=arcpy.SearchCursor(inTbl)


for row in rows:
att1 = row.getValue("TILENM_400")
att2 = row.getValue("SUBBLK_200")
att3 = row.getValue("TILENM_200")
xmin = row.getValue("X_Min")
xmax = row.getValue("X_Max")
ymin = row.getValue("Y_Min")
ymax = row.getValue("Y_Max")
extent = arcpy.Extent(xmin,ymin,xmax,ymax)
array.add(extent.lowerLeft)
array.add(extent.lowerRight)
array.add(extent.upperRight)
array.add(extent.upperLeft)
#ensure Polygon is closed
array.add(extent.lowerLeft)
polygon = arcpy.Polygon(array)
ptList.append(polygon)
array.removeAll()

del rows
arcpy.CopyFeatures_management(ptList, outShp)

Anonymous said...

Hey Andrew,
I've been working with your python script trying to build bounding boxes for GPS locations. Basically I have 24,000 multipart points and I want each of the 24,000 to have a bounding box.

The easy answer is just use ET Geowizards, unfortunately I can't because it's 3rd party software and that's a big no no at my company.

Your script has been huge so far but I'm getting returns on the second part that I just have no answer for. The two problems are:
1) Some of the bounding boxes return with 0 value in the shape_length and shape_area columns, and have no presence in the map itself. They are empty features. Maybe the area is too small to be calculated? Should I run the extent again?
2)The bounding boxes don't have the same attributes as my multipoint shapefile. I can fix this without python, but is this a normal problem?

Thanks for any help,
Alex

Andrew said...

@alex, so you are asking a lot. If you have a multi-point group with only 1 feature, then the extent is most likely none, or you have empty geometry feature classes. An easy test of this, is to copy your data (very important) to a another location and run repair geometry on the data set with erase NULL geometries checked.

To get the attributes from your current feature would require python most likely. You can do it two ways. Using a searchcursor from the source dataset, replace the geometry object with the new polygon extent. Then you would insert it into a polygon feature class. Way 2, you might be able to do a spatial join on your polygon extents with the multi-point feature classes.

I haven't tested any of these, but I'm just throwing it out there.

Good luck and post your results!

Tim said...

Dear Andrew,
Would the same rectangle shapes be possible, with area of the respective polygons (instead of extent)?! And it should be a predefined width of the rectangle.. This would be completely amazing for my work...

Andrew said...

You mean create the polygon rectangle but put the area of the shape as an attribute?

Tim said...

hm..I want ArcPy to create rectangles according to the area (instead of extent) of the respective polygons. Means create a rectangle which is length*width=area of the respective polygons..
Furthemore.. I want to define a fixed length of the respective resulting rectangular so that the surface area varies only in width (and still has the same surface area as the resoective basis Polygon)
Thanks, Tim!

Andrew said...

Tim,

How would you know what the width and length values are?
Unless you wanted squares, then it is easy.

1TIGA11 said...

Dear Andrew,
Width and length should actually come from an excel table (I know its possible to import excel tables as dbf to arcgis)
Do you have an idea?!

Best Tim

Andrew said...

Take your width and divide by the area to get your Y length. Use this to build all your polygons.