[GIS] Get extent of Georeferenced Rasters in Python and output to polygon shapefile

arcpyerror-000725error-999999gdalraster

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.

Best Answer

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.

Related Solutions

[GIS] Inserting DATAFRAME_ELEMENT value into dynamically generated shapefile containing its extent

I'm a bit lost here, so my answer may be off base. When you say you want the ExtentID, that appears to be the name of the Dataframe you're exporting the extent of. If this is the case (which is what your code is doing to my best guess), you don't need the [x] index when printing it. :

 for el in arcpy.mapping.ListLayoutElements(mxd2, "DATAFRAME_ELEMENT"):
     ExtentIDStore = str(el.name)
     #print x, y, z <--- You don't need to print x y and z. They are unrelated at this point.
     print ExtentIDStore #used to see which item to add to the shape file

     print ExtentIDStore[x] # <- This will just print an indexed character in the string

By the time you reach the loop above, your code is only concerned with the Dataframe Element in question, so there is no need to do any indexing for storage. ExtendIDStore should be what you add to the shapefile row.

Your error which results at this line:

row.strField = str(x)

Is probably because you are not referencing a field that exists. What comes after the row reference is the name of the Field. In this case, it should be:

row.ExtentID = ExtentIDStore

File name you can get via the filePath property of the MapDocument class and the os.path.basename method, i.e.:

filepath = mxd2.filePath # Gets full path of MXD used in your dataframe loop
filename = os.path.basename(filepath) # Extracts the file name and extension
row.FileRef = filename

If you don't want the extension, use filename.strip(".mxd") to remove it.

Scale would be:

scale = el.scale # Scale of the dataframe referenced in the loop
row.Scale = scale

And that should fill in the rows you're looking for.

For the new problem, I think you are misusing the insert cursor. The second argument should be a spatial reference, which cannot be gleaned via an MXD. I would recommend creating the cursor with only one argument (the output file or dfStore). You also forgot the last (and most important) line:

rows.insertRow(row)

Without that your data never gets put into the table. Put that down at the end of the row loop so that the data is actually inserted after all the variables are set.

For more information about accessing Dataframe, layer, and map properties (scale, name, current position), check out the mapping module. The scale property of the Dataframe class can be seen on the classes page of the same name.

[GIS] ogr2ogr DXF to PostgreSQL some LINESTRING coordinates mirrored to the negative longitude

I can't find the origin of this problem but I finally reached the final objective through another tool, with the GIS:

1) open the DXF file with QGIS

2) convert the DXF to Shapefile (right click on the layer > "Save as...")

3) use the Shapefiles to PostgreSQL converter in the menu option: Database > Split > "Import Shapefiles to PostgreSQL"

Trough this tool I imported the DXF to PostgreSQL without the wrong coordinates problem. I suppose that in the DXF to Shapefiles conversion the origin of the problem is solved.

If anyone have more information about what can origin that situation when I use the ogr2ogr tool it would be nice to know, because right know I cannot make that conversion (DXF to PostgreSQL) automatic (with an script that would execute the ogr2ogr command).

Best Answer

Related Solutions

[GIS] Inserting DATAFRAME_ELEMENT value into dynamically generated shapefile containing its extent

[GIS] ogr2ogr DXF to PostgreSQL some LINESTRING coordinates mirrored to the negative longitude

Related Question