Defining Raster Extents using Python Geoprocessing Script

4398
2
03-30-2011 04:13 PM
KatrinaBennett
New Contributor II
Hello, I'm trying to import and define the extents of raster GeoTIFF images. For some reason, after the update of ArcGIS 10 SP1 I am not able to define extents of these images such that the extents remain in the raster images after I convert them.

I've specified my extents several times through my code, but I don't understand why this doesn't 'take' in the final product.

Any ideas? Thanks for your thoughts in advance... My Python code is posted below.

Katrina



import arcgisscripting, sys, os, shutil, datetime, re, arcpy

gp = arcgisscripting.create()
gp.overwriteoutput = 1
gp.CheckOutExtension("spatial")
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")

sourceDirRas = r"G:\modis\updated_example\my_own\albers"
extent = "-1805197.82136111 374274.164059672 1712802.17863889 2569774.16405967"
i = 1

for root, dirs, files in os.walk(sourceDirRas):
   for file in files:
      if file[-4:] == ".tif":       
        wholeFile = root + "\\" + file
        outFile = root + "\\out\\" + "ms" + str(i)
        print file

        #Calculate Stats
        print "Calculating Stats!"
        gp.CalculateStatistics_management(wholeFile, "", "", "")

        print "Setting Extent and Define Projection!"
        tempEnvironment0 = arcpy.env.extent
        arcpy.env.extent = "-1805197.82136111 374274.164059672 1712802.17863889 2569774.16405967"
        #gp.Extent = "-1805197.82136111 374274.164059672 1712802.17863889 2569774.16405967"
        arcpy.DefineProjection_management(wholeFile, "PROJCS['AEA_North_American_Datum_1983',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',-154.0],PARAMETER['standard_parallel_1',55.0],PARAMETER['standard_parallel_2',65.0],PARAMETER['latitude_of_origin',50.0],UNIT['Meter',1.0]]")
        arcpy.env.extent = tempEnvironment0

        print "Copy Raster!"
        tempEnvironment0 = arcpy.env.extent
        arcpy.env.extent = "-1805197.82136111 374274.164059672 1712802.17863889 2569774.16405967"
        #gp.CopyRaster_management(wholeFile, outFile, "#","0","#","NONE","NONE","#")
        arcpy.CopyRaster_management(wholeFile, outFile, "", "", "", "NONE", "NONE", "")
        arcpy.env.extent = tempEnvironment0

        i = i + 1
        print "Done!"
2 Replies
KatrinaBennett
New Contributor II
I made a few changes to my code and it works now! I'm unsure what was happening before...

import arcgisscripting, sys, os, shutil, datetime, re

gp = arcgisscripting.create()
gp.overwriteoutput = 1
gp.CheckOutExtension("spatial")
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")

sourceDirRas = r"G:\modis\updated_example\my_own\albers"
extent = "-1805197.82136111 374274.164059672 1712802.17863889 2569774.16405967"
i = 1

for root, dirs, files in os.walk(sourceDirRas):
   for file in files:
      if file[-4:] == ".tif":       
        wholeFile = root + "\\" + file
        outFile = root + "\\out\\" + "ms" + str(i)
        print file

        #Calculate Stats
        print "Calculating Stats!"
        gp.CalculateStatistics_management(wholeFile, "", "", "")

        print "Setting Extent and Define Projection!"
        gp.Extent = extent
        gp.DefineProjection_management(wholeFile, "PROJCS['AEA_North_American_Datum_1983',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',-154.0],PARAMETER['standard_parallel_1',55.0],PARAMETER['standard_parallel_2',65.0],PARAMETER['latitude_of_origin',50.0],UNIT['Meter',1.0]]")

        print "Copy Raster!"
        gp.Extent = extent
        gp.CopyRaster_management(wholeFile, outFile, "", "", "", "NONE", "NONE", "")

        i = i + 1
        print "Done!"
0 Kudos
KarenMaloof1
New Contributor

Thanks for posting your solution. I've been banging my head on this problem all day, so a solution that works for an inexplicable reason is great.

0 Kudos