# ---------------------------------------------------------------------------

# MakeGridCells.py

# August 4,2009

# Usage: MakeGridCells <Workspace> <GDB_Name> <FDS_Name> <Input_Grid> <PolysFC_Name>

# ---------------------------------------------------------------------------

#Import system modules
import sys, string, os, arcgisscripting, math

def Fishnet(file_array,Workspace,GDB_Name,Poly_Names,Cell_Width,Cell_Height):
    #Create folder to store files if DNE
    try:
        os.makedirs(Workspace + "/Fishnet")
        fishnet_Workspace = Workspace + "/Fishnet/" 
        
    except OSError:
        fishnet_Workspace = Workspace + "/Fishnet/"
    
    #Create the Geoprocessor and set overwrite to true
    gp = arcgisscripting.create(9.3)
    gp.CheckOutExtension("Spatial")
    gp.OverwriteOutput = True
    
    #Process script arguments and derive some variables
    GDB = GDB_Name + ".gdb"
    FDS_Name = "AggregateCells" #Provide a default value if unspecified
    FDS = GDB + "/" + FDS_Name
    Input_Grid = file_array[0]

    #Set the Workspace
    gp.Workspace = fishnet_Workspace
    print Input_Grid
    #Determine spatial reference of Input_Grid for feature data set
    desc = gp.Describe(Input_Grid)
    SpRef = desc.SpatialReference
    
    #Verify that the Input_Grid is in projected space.
    if SpRef.Type != "Projected":
        print "Input raster data is not a projected data set."
        sys.exit()
    
    #Create the Geodatabase and empty Feature Data Set.
    print "Creating Geodatabase " + GDB_Name + " with Feature Data Set " + FDS_Name + "..."
    gp.CreateFileGDB_management(fishnet_Workspace,GDB_Name)
    gp.CreateFeatureDataset_management(GDB,FDS_Name,SpRef)
    
    #Specify names for feature classes of Fishnet Lines and Lables.
    FishnetFC = FDS +"/Fishnet"
    LabelsFC = FDS + "/Fishnet_label"
    
    #Describe Input_Grid and derive extents.
    desc = gp.Describe(Input_Grid)
    Extent = desc.Extent
    OriginPt = str(Extent.XMin) + " " + str(Extent.YMin)
    AxisPt = str(Extent.XMin) + " " + str(Extent.YMin + 10)
    
    #Calculate rows and columns needed for fishnet
    Rows = int(math.ceil((Extent.YMax - Extent.YMin)/int(Cell_Height)))
    Cols = int(math.ceil((Extent.XMax - Extent.XMin)/int(Cell_Width)))
    
    #Create fishnet lines and labels
    print "Creating Fishnet with " + str(Rows) + " rows and " + str(Cols) + " columns at Origin Point:" + str(OriginPt) + "..."
    gp.CreateFishnet_management(FishnetFC, OriginPt, AxisPt, Cell_Width, Cell_Height, Rows, Cols, "#", "labels", Input_Grid)
    
    #Make feature layers of the Fishnet and Labels feature classes.
    FishnetLyr = "FishnetLayer"
    LabelsLyr = "LabelsLayer"
    gp.MakeFeatureLayer_management(FishnetFC, FishnetLyr)
    gp.MakeFeatureLayer_management(LabelsFC, LabelsLyr)
    
    #Add X,Y coordinates to labels
    print "Getting coordinates for cells..."
    gp.AddXY_management(LabelsLyr)
    
    #Create Polygons feature class from fishnet lines and labels
    iter = 0
    for raster in file_array:
        PolysFC_Name = Poly_Names[iter]
        PolysFC = FDS + "/" + PolysFC_Name
        Input_Grid = raster
        print "Creating " + PolysFC_Name
        gp.FeatureToPolygon_management(FishnetLyr,PolysFC,"#","Attributes",LabelsLyr)
        
        #Add ZONE_ID number and set it equal to Object ID
        print "Assigning Zone IDs..."
        gp.AddField_management(PolysFC,"ZONE_ID","long")
        gp.CalculateField_management(PolysFC,"ZONE_ID","[OBJECTID]")
        
        #Calculate zonal statistics as a table.
        print "Calculating zonal statistics..."
        Stat_Table = GDB + "\ZonalStats"
        gp.ZonalStatisticsAsTable_sa(PolysFC,"ZONE_ID",Input_Grid,Stat_Table,"DATA")
        
        #Join ZonalStats table to PolysFC
        print "Joining tables"
        gp.JoinField_management(PolysFC,"ZONE_ID",Stat_Table,"ZONE_ID")
        
        iter += 1
        
    #Create Polygons for Head,Error, and Extraction
    calcs=['Head','Error','Extraction']
    for name in calcs:
        PolysFC = FDS + "/" + name
        print "Creating " + name
        gp.FeatureToPolygon_management(FishnetLyr,PolysFC,"#","Attributes",LabelsLyr)
        gp.AddField_management(PolysFC,str(name),"float")
        
    
    #Delete variables, releasing memory.
    del gp, Workspace, GDB_Name, GDB, FDS_Name, FDS, PolysFC_Name, PolysFC, Input_Grid, Cell_Width, Cell_Height
    del desc, SpRef, FishnetFC, LabelsFC, Extent, OriginPt, AxisPt, Rows, Cols, FishnetLyr, LabelsLyr, Stat_Table
    
    print "Finished!"

