3

I would like to create a shapefile containing the extents of each of the rasters in a directory. Is it possible to capture the extent of a raster using Python?

I have tried

extent1=arcgisscripting.Raster.extent('stg1_05.jpg') Runtime error : 'getset_descriptor' object is not callable

and I can't seem to find any help on the module.

I also tried

arcpy.RasterToPolygon_conversion(inRaster, outPolygons, "SIMPLIFY", field)

Runtime error : ERROR 010247: The size of the output shapefile exceeds the 2 GB limit. ERROR 010067: Error in executing grid expression. Failed to execute (RasterToPolygon).

Anyway this will be a polygon of the whole raster file and not just the extents -even if this is generated, I guess we could then run a merge/dissolve but the files created are to big.

Another option I was thinking of was to convert the raster to layer

import arcpy, glob, os, sys
from arcpy import env, mapping
path = os.getcwd()
env.workspace = path
RasterType = 'TIF'
FileList = glob.glob(path + "\*." + RasterType)

print 'Reading files from ' + path
os.chdir(path)
#print FileList

x=0
z=1005
File=FileList[x]
LayerWorking=arcpy.mapping.Layer(File)
print File
LayerExtent=LayerWorking.getExtent()
XMAX = LayerExtent.XMax
XMIN = LayerExtent.XMin
YMAX = LayerExtent.YMax
YMIN = LayerExtent.YMin
pnt1 = arcpy.Point(XMIN, YMIN)
pnt2 = arcpy.Point(XMIN, YMAX)
pnt3 = arcpy.Point(XMAX, YMAX)
pnt4 = arcpy.Point(XMAX, YMIN)
array = arcpy.Array()
array.add(pnt1)
array.add(pnt2)
array.add(pnt3)
array.add(pnt4)
array.add(pnt1)
polygon = arcpy.Polygon(array)
ShapeFile = path + "\Polygon_Extent" + "_" + str(z) + ".shp"
print ShapeFile
print arcpy.GetMessages()
arcpy.CreateFeatureclass_management(path, ShapeFile, "POLYGON")
print arcpy.GetMessages()
arcpy.CopyFeatures_management(polygon, ShapeFile)

Gives following output

Reading files from C:\Python26\script_tests\rectify C:\Python26\script_tests\rectify\Stage1_01a.TIF C:\Python26\script_tests\rectify\Polygon_Extent_1005.shp Executing: CopyFeatures in_memory\fB4DC2172_7D02_44B9_B55C_9E71427AE96E C:\Python26\script_tests\rectify\Polygon_Extent_1004.shp # 0 0 0 Start Time: Thu Jul 14 16:11:38 2011 Failed to execute. Parameters are not valid. ERROR 000725: Output Feature Class: Dataset C:\Python26\script_tests\rectify\Polygon_Extent_1004.shp already exists. Failed to execute (CopyFeatures). Failed at Thu Jul 14 16:11:38 2011 (Elapsed Time: 0.00 seconds) Traceback (most recent call last): File "C:\Python26\script_tests\rectify\temp.py", line 36, in arcpy.CreateFeatureclass_management(path, ShapeFile, "POLYGON") File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\management.py", line 1539, in CreateFeatureclass raise e ExecuteError: ERROR 999999: Error executing function. Failed to execute (CreateFeatureclass).

using IDLE...Which is wierd as when you run it from the Python window in ArcMap it works fine when the Create/Copy feature commands are run individually.

ShapeFile = "Polygon_Extent" + "_" + str(z) + ".shp" arcpy.CreateFeatureclass_management(path, ShapeFile, "POLYGON")

Result 'C:\Python26\script_tests\rectify\Polygon_Extent_1005.shp'>

arcpy.CopyFeatures_management(polygon, ShapeFile)

Result 'C:\Python26\script_tests\rectify\Polygon_Extent_1005.shp'>

this seems to be a very convoluted process...

UPDATE: I am trying to add the filename to the table and for some reason it inserts a new row into the table and doesn't accept "updateRow(row)"? what am I doing wrong? Also the files don't seem to retain the projection I assign them.

Reading files from C:\Python26\script_tests\rectify C:\Python26\script_tests\rectify\Stage1_01a.TIF Created: Polygon_Extent_1.shp Filled in C:\Python26\script_tests\rectify\Stage1_01a.TIF Traceback (most recent call last): File "C:\Python26\script_tests\rectify\ResterExtent_toSHP.py", line 45, in rows.updateRow(row) File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\arcobjects.py", line 102, in updateRow return convertArcObjectToPythonObject(self._arc_object.UpdateRow(*gp_fixargs(args))) RuntimeError: ERROR 999999: Error executing function.

SCRIPT all working now.

import arcpy, glob, os, sys, arcgisscripting
from arcpy import env, mapping
path = os.getcwd()
env.workspace = path
env.overwriteOutput = True
RasterType = 'TIF'
FileList = glob.glob(path + "\*." + RasterType)
print 'Reading files from ' + path
os.chdir(path)
#print FileList
geometry_type = "POLYGON"
template = "C:\\Python26\\GDA_1994_MGA_Zone_55.shp"
has_m = "DISABLED"
has_z = "DISABLED"
# Creating a spatial reference object
spatial_reference = arcpy.SpatialReference("C:\\Python26\\GDA_1994_MGA_Zone_55.prj")

x=0
z=x+1
for File in FileList:
#File=FileList[x]
    RasterFile = arcgisscripting.Raster(File)
    RasterExtent = RasterFile.extent
    print File
    XMAX = RasterExtent.XMax
    XMIN = RasterExtent.XMin
    YMAX = RasterExtent.YMax
    YMIN = RasterExtent.YMin
    pnt1 = arcpy.Point(XMIN, YMIN)
    pnt2 = arcpy.Point(XMIN, YMAX)
    pnt3 = arcpy.Point(XMAX, YMAX)
    pnt4 = arcpy.Point(XMAX, YMIN)
    array = arcpy.Array()
    array.add(pnt1)
    array.add(pnt2)
    array.add(pnt3)
    array.add(pnt4)
    array.add(pnt1)
    polygon = arcpy.Polygon(array)
    arcpy.CreateFeatureclass_management(path, "Temp_Polygon_Extent" + "_" + str(z), geometry_type, template, has_m, has_z, spatial_reference)
    arcpy.CopyFeatures_management(polygon, "Temp_Polygon_Extent" + "_" + str(z))
    ShapeFile = "Temp_Polygon_Extent" + "_" + str(z) + ".shp"
    print "Created: " + ShapeFile
    arcpy.AddField_management(ShapeFile,'FileName','TEXT')
    desc = arcpy.Describe(ShapeFile)
    print desc, ShapeFile
    #rows = arcpy.InsertCursor(ShapeFile, desc)
    rows = arcpy.UpdateCursor(ShapeFile)
    #row = rows.newRow()
    #row.FileName = str(File)
    #row.FileName = File
    print 'Filled in: ', str(File) 
    #rows.insertRow(row)
    for row in rows:
        row.FileName = str(ShapeFile)
        rows.updateRow(row)
    x=x+1
    z=z+1

#cleaning up
arcpy.CreateFeatureclass_management(path, "Extent.shp", geometry_type, template, has_m, has_z, spatial_reference)
list =[]
lstFCs = arcpy.ListFeatureClasses("Temp_Polygon_Extent*")
print 'Merging Polygon_Extents* to Extent.shp'

for fc in lstFCs:
    print fc
    list.append(fc)

arcpy.Merge_management(list, "Extent.shp")
#print 'Deleting identical entries and temp files'
#arcpy.DeleteIdentical_management("Extent.shp", ["SHAPE"])
print 'Created Extent.shp and exiting'

for item in list:
    arcpy.Delete_management(item) 
del row, rows

GDAL Ver.

import os, gdal
gdaltindex Extent.shp *.tif

File "". line 1 SyntaxError: invalid syntax

BTW FTTools Python window is a pain as you can't copy/past code/errors from it to.

2 Answers 2

6

You almost had it right with your first attempt. What you are doing is trying to call the extent property with a filename as its parameter, when you need to construct a Raster object with that parameter.

In practical terms, that means:

extent1 = arcgisscripting.Raster('stg1_05.jpg').extent

Although it is usually better practice to break it down into two steps:

raster1 = arcgisscripting.Raster('stg1_05.jpg')
extent1 = raster1.extent

Note the lack of parentheses on extent, this is because it is a property rather than a method.

2
  • Thanks mate. That worked great. I have uploaded the new code for anyone to use. Useful to generate the extents of all images in a directory so that they can be easily referenced. I will expand this to os.walk and all imagery types in due course...but I am having an issue with insert/update.Row. Please see updated post.
    – GeorgeC
    Commented Jul 15, 2011 at 1:25
  • CORRECTED IT...all working fine now. Thanks Mersey.
    – GeorgeC
    Commented Jul 15, 2011 at 6:59
2

GDAL has a gdaltindex command-line utility that does the same thing you are after:

gdaltindex Extent.shp *.jpg

The utility was intended for tile indexing for MapServer, but it works just the same for any other purpose. One caveat (which is the same for your script) is that it does not rotate the polygon footprint for rotated rasters.

4
  • thanks. I haven't been able to get GDAL to work in my ArcGIS python install yet so I haven't been able to try any GDAL based solutions.
    – GeorgeC
    Commented Jul 15, 2011 at 5:12
  • Yeah, it's tricky to install older GDAL versions for the same version of Python from ArcGIS. I install OSGeo4W to keep it separate from Esri software, then run the commands in their "OSGeo4W Shell". FWTools also work just the same.
    – Mike T
    Commented Jul 15, 2011 at 5:22
  • Thanks mike. I do have GDAL in my FWTOOL (forgot that I did). I am getting errors though. See the ticket please.
    – GeorgeC
    Commented Jul 15, 2011 at 5:38
  • The commands need to be used in a FWTools Shell .. not a Python console
    – Mike T
    Commented Jul 15, 2011 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.