#PanSharpen QuarterQuad Clip

#Created by:  Chris Curlis
#geographer_4@hotmail.com
#Creation Date:  4/20/2010
# Revised Date:  5/4/2010

#This script takes an existing Landsat 7 image and creates a PanSharpened image from 4 bands of the imagery along with the panchromatic band
#It then creates shapefile masks from an existing shapefile of quarterquads, and clips the Pansharpened image for each of the quarterquad shapefiles.
#This script requires the Spatial Analyst Extension to be loaded

#Portions of the script setup and except came from Nathan Jennings
# nate@jenningsplanet.com
# www.jenningsplanet.com

#Portion for directory creation was obtained from Douglas Mayle, (Douglas.Mayle.org) ,answer #7 at:
#(http://stackoverflow.com/questions/273192/python-best-way-to-create-directory-if-it-doesnt-exist-for-file-write)

#Portion was obtained from Robert Claypool, robert.claypool@amec.com  from his response to an inquiry at "Clip a raster dataset against a fishnet":
#http://groups.google.com/group/esri/browse_thread/thread/a70a918d7adc96fb?utoken=ybJ2FicAAAAFGaYcnKQ8OO3624bCIEfpwg-L0Ck4OY0txtIHwwRAqAdp9zfgqaDRkjVr7gxBuhM


#The general steps are:

#Setup and load all workspaces and needed modules, licenses etc.
#Check for existance of files from previous attempts and delete them
#Define variables
#Create PanSharpened image
#Create search cursor for the predefined shapefile of quarterquads of interest (those without gaps in a Landsat 7 scene)
#Delete the directory "outputs" created from previous attempts
#Create a new "outputs" directory
#Select and name the output quarterquad shapefiles then save them to new shapefiles
#Use new shapefiles as masks to clip the newly created PanSharpened image by quarterquad


#Import standard/necessary library modules
import arcgisscripting, sys, os, string, traceback
  
#Create the Geoprocessor object

gp=arcgisscripting.create(9.3)

#Set the input/output workspace
gp.workspace = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\"

#Set image path (may not be required)

image_path = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\"

# Check out any necessary licenses
gp.CheckOutExtension("spatial")

# Load required toolboxes...
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")

# Set the Geoprocessing environment...

# Not absolutely needed

gp.outputCoordinateSystem = ""
gp.rasterStatistics = "STATISTICS 1 1"
gp.pyramid = "PYRAMIDS -1 NEAREST"

#Check existance of and delete pansharpened image

if gp.exists("C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\44033_03320071102_all_PS.img"):
    gp.delete_management("C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\44033_03320071102_all_PS.img")
print "Previous PanSharpened image deleted"

# Local variables...this comes from a modelbuilder derived script
v44033_03320071102_all_PS_img = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\44033_03320071102_all_PS.img"
v44033_03320071102_all_img = "44033_03320071102_all.img"
l72044033_03320071102_b80_img = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\l72044033_03320071102_b80.img"

# Process: Create Pan-sharpened Raster Dataset...set variables for R, G, B and IR bands
gp.CreatePansharpenedRasterDataset_management(v44033_03320071102_all_img, "3", "2", "1", "4", v44033_03320071102_all_PS_img, l72044033_03320071102_b80_img, "ESRI", "0.166", "0.167", "0.167", "0.5")
print "Landsat Image has been PanSharpened"

# Create a search cursor for the QuarterQuad subset shapefile which will loop through all the records in the shapefile
Quad_Polygons = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\NoGaps.shp"
rows = gp.SearchCursor(Quad_Polygons)
row = rows.Next()

# Delete the "outputs" directory if it exists
if gp.exists("C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\Outputs\\"):
    gp.delete_management("C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\Outputs\\")
print "Previous outputs directory deleted"

# If directory "Outputs" doesn’t exist, create it  (this came from Douglas Mayle on stackusers.com (see above))

import os
try:
    os.makedirs('C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\Outputs')
except OSError:
    pass

print "New outputs directory created"

#Select and Name the quarterquads by their "USGSID" combined with their "Qdrnt"
#this is modified from Robert Claypool script (see above)
try:
    while row:
        # Find the FID field
        FID = row.GetValue("FID")

        # Find the USGSID and Qdrnt fields.  Get values for these rows. These values will be used to select the feature and name output files.
        Name = row.GetValue("USGSID")+ row.GetValue("Qdrnt")

        # Generate an expression to select the feature
        Expression = '\"FID\" = ' + str(FID) # Where "FID"=FID

        # Execute Process: Select the feature for this row, save it to a new shapefile
        Output_Path_1 = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\Outputs\\" + Name + "QQ.shp"
        gp.Select_analysis(Quad_Polygons, Output_Path_1, Expression)
        print "Shapefile for QQuad created" + Name
        
        # Execute Process: Clip Image using the new shapefile as a mask
        PanImg = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\44033_03320071102_all_PS.img"
        Output_Path_2 = "C:\\Documents and Settings\\Owner\\My Documents\\Geog376\\FinalProject\\Outputs\\" + Name + ".img"
        gp.ExtractByMask_sa(PanImg, Output_Path_1, Output_Path_2)
        print "PanSharpened QQuad clipped" + Name
        
        # Advance the cursor for next iteration of the loop
        row = rows.Next()

    print "Script completed...all quarterquads have been clipped"   

except:
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n     " +        str(sys.exc_type) + ": " + str(sys.exc_value) + "\n"
    msgs = "gp ERRORS:\n" + gp.GetMessages(2) + "\n"

    gp.AddError(msgs)
    gp.AddError(pymsg)

    print msgs
    print pymsg
    
    gp.AddMessage(gp.GetMessages(1))
    print gp.GetMessages(1)