# SSURGO_ExportMuRaster.py
#
# Convert MUPOLYGON featureclass to raster for the specified SSURGO geodatabase.
# By default any small NoData areas (< 5000 sq meters) will be filled using
# the Majority value.
#
# Input mupolygon featureclass must have a projected coordinate system or it will skip.
# Input databases and featureclasses must use naming convention established by the
# 'SDM Export By State' tool.
#
# For geographic regions that have USGS NLCD available, the tool wil automatically
# align the coordinate system and raster grid to match.
#
# 10-31-2013 Added gap fill method
#
# 11-05-2014
# 11-22-2013
# 12-10-2013  Problem with using non-unique cellvalues for raster. Going back to
#             creating an integer version of MUKEY in the mapunit polygon layer.
# 12-13-2013 Occasionally see error messages related to temporary GRIDs (g_g*) created
#            under "C:\Users\steve.peaslee\AppData\Local\Temp\a subfolder". These
#            are probably caused by orphaned INFO tables.
# 01-08-2014 Added basic raster metadata (still need process steps)
# 01-12-2014 Restricted conversion to use only input MUPOLYGON featureclass having
#            a projected coordinate system with linear units=Meter
# 01-31-2014 Added progressor bar to 'Saving MUKEY values..'. Seems to be a hangup at this
#            point when processing CONUS geodatabase
# 02-14-2014 Changed FeatureToLayer (CELL_CENTER) to PolygonToRaster (MAXIMUM_COMBINED_AREA)
#            and removed the Gap Fill option.
# 2014-09-27 Added ISO metadata import
#
# 2014-10-18 Noticed that failure to create raster seemed to be related to long
# file names or non-alphanumeric characters such as a dash in the name.
#
# 2014-10-29 Removed ORDER BY MUKEY sql clause because some computers were failing on that line.
#            Don't understand why.
#
# 2014-10-31 Added error message if the MUKEY column is not populated in the MUPOLYGON featureclass
#
# 2014-11-04 Problems occur when the user's gp environment points to Default.gdb for the scratchWorkpace.
#            Added a fatal error message when that occurs.
#
# 2015-01-15 Hopefully fixed some of the issues that caused the raster conversion to crash at the end.
#            Cleaned up some of the current workspace settings and moved the renaming of the final raster.
#
# 2015-02-26 Adding option for tiling raster conversion by areasymbol and then mosaicing. Slower and takes
#            more disk space, but gets the job done when otherwise PolygonToRaster fails on big datasets.

# 2015-02-27 Make bTiling variable an integer (0, 2, 5) that can be used to slice the areasymbol value. This will
#            give the user an option to tile by state (2) or by survey area (5)
# 2015-03-10 Moved sequence of CheckInExtension. It was at the beginning which seems wrong.
#
# 2015-03-11 Switched tiled raster format from geodatabase raster to TIFF. This should allow the entire
#            temporary folder to be deleted instead of deleting rasters one-at-a-time (slow).
# 2015-03-11 Added attribute index (mukey) to raster attribute table
# 2015-03-13 Modified output raster name by incorporating the geodatabase name (after '_' and before ".gdb")

## ===================================================================================
class MyError(Exception):
    pass

## ===================================================================================
def PrintMsg(msg, severity=0):
    # prints message to screen if run as a python script
    # Adds tool message to the geoprocessor
    #
    #Split the message on \n first, so that if it's multiple lines, a GPMessage will be added for each line
    try:
        for string in msg.split('\n'):
            #Add a geoprocessing message (in case this is run as a tool)
            if severity == 0:
                arcpy.AddMessage(string)

            elif severity == 1:
                arcpy.AddWarning(string)

            elif severity == 2:
                arcpy.AddMessage("    ")
                arcpy.AddError(string)

    except:
        pass

## ===================================================================================
def errorMsg():
    try:
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]
        theMsg = tbinfo + "\n" + str(sys.exc_type)+ ": " + str(sys.exc_value)
        PrintMsg(theMsg, 2)

    except:
        PrintMsg("Unhandled error in errorMsg method", 2)
        pass

## ===================================================================================
def setScratchWorkspace():
    # This function will set the scratchGDB environment based on the scratchWorkspace environment.
    #
    # The scratch workspace will typically not be defined by the user but the scratchGDB will
    # always be defined.  The default location for the scratchGDB is at C:\Users\<user>\AppData\Local\Temp
    # on Windows 7 or C:\Documents and Settings\<user>\Localsystem\Temp on Windows XP.  Inside this
    # directory, scratch.gdb will be created.
    #
    # If scratchWorkspace is set to something other than a GDB or Folder than the scratchWorkspace
    # will be set to C:\temp.  If C:\temp doesn't exist than the ESRI scratchWorkspace locations will be used.
    #
    # If scratchWorkspace is an SDE GDB than the scratchWorkspace will be set to C:\temp.  If
    # C:\temp doesn't exist than the ESRI scratchWorkspace locations will be used.

    try:
        # -----------------------------------------------
        # Problems can occur (especially raster gp) when environment ScratchWorkspace = Default.gdb
        if os.path.basename(env.scratchGDB) == "Default.gdb":
            raise MyError, "Please set geoprocessing environment 'Scratch Workspace' to another database. Default.gdb not allowed"

        if arcpy.env.scratchWorkspace is not None:

            # describe scratch workspace
            scratchWK = arcpy.env.scratchWorkspace
            descSW = arcpy.Describe(scratchWK)
            descDT = descSW.dataType.upper()

            # scratch workspace is geodatabase
            if descDT == "WORKSPACE":
                progID = descSW.workspaceFactoryProgID

                # scratch workspace is a FGDB
                if  progID == "esriDataSourcesGDB.FileGDBWorkspaceFactory.1":
                    arcpy.env.scratchWorkspace = os.path.dirname(scratchWK)
                    arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

                # scratch workspace is a Personal GDB -- set scratchWS to folder of .mdb
                elif progID == "esriDataSourcesGDB.AccessWorkspaceFactory.1":
                    arcpy.env.scratchWorkspace = os.path.dirname(scratchWK)
                    arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

                # scratch workspace is an SDE GDB.
                elif progID == "esriDataSourcesGDB.SdeWorkspaceFactory.1":

                    # set scratch workspace to C:\Temp; avoid the server
                    if os.path.exists(r'C:\Temp'):

                        arcpy.env.scratchWorkspace = r'C:\Temp'
                        arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

                    # set scratch workspace to default ESRI location
                    else:
                        arcpy.env.scratchWorkspace = arcpy.env.scratchFolder
                        arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

            # scratch workspace is simply a folder
            elif descDT == "FOLDER":
                arcpy.env.scratchWorkspace = scratchWK
                arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

            # scratch workspace is set to something other than a GDB or folder; set to C:\Temp
            else:
                # set scratch workspace to C:\Temp
                if os.path.exists(r'C:\Temp'):

                    arcpy.env.scratchWorkspace = r'C:\Temp'
                    arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

                # set scratch workspace to default ESRI location
                else:
                    arcpy.env.scratchWorkspace = arcpy.env.scratchFolder
                    arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

        # -----------------------------------------------
        # Scratch Workspace is not defined. Attempt to set scratch to C:\temp
        elif os.path.exists(r'C:\Temp'):
            arcpy.env.scratchWorkspace = r'C:\Temp'
            arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

        # set scratch workspace to default ESRI location
        else:
            arcpy.env.scratchWorkspace = arcpy.env.scratchFolder
            arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

        return True

    except MyError, e:
        # Example: raise MyError, "This is an error message"
        PrintMsg(str(e) + " \n ", 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def SnapToNLCD(inputFC, iRaster):
    # This function will set an output extent that matches the NLCD raster dataset.
    # In effect this is like using NLCD as a snapraster as long as the projections are the same,
    # which is USA_Contiguous_Albers_Equal_Area_Conic_USGS_version
    #
    # Returns empty string if linear units are not 'foot' or 'meter'

    try:
        theDesc = arcpy.Describe(inputFC)
        sr = theDesc.spatialReference
        inputSRName = sr.name
        theUnits = sr.linearUnitName
        pExtent = theDesc.extent

        PrintMsg(" \nCoordinate system: " + inputSRName + " (" + theUnits.lower() + ")", 0)

        if pExtent is None:
            raise MyError, "Failed to get extent from " + inputFC

        x1 = float(pExtent.XMin)
        y1 = float(pExtent.YMin)
        x2 = float(pExtent.XMax)
        y2 = float(pExtent.YMax)

        if 'foot' in theUnits.lower():
            theUnits = "feet"

        elif theUnits.lower() == "meter":
            theUnits = "meters"

        # USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (NAD83)
        xNLCD = 532695
        yNLCD = 1550295

        # Hawaii_Albers_Equal_Area_Conic  -345945, 1753875
        # Western_Pacific_Albers_Equal_Area_Conic  -2390975, -703265 est.
        # NAD_1983_Alaska_Albers  -2232345, 344805
        # WGS_1984_Alaska_Albers  Upper Left Corner:  -366405.000000 meters(X),  2380125.000000 meters(Y)
        # WGS_1984_Alaska_Albers  Lower Right Corner: 517425.000000 meters(X),  2032455.000000 meters(Y)
        # Puerto Rico 3092415, -78975 (CONUS works for both)

        if theUnits != "meters":
            PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
            raise MyError, "Unable to align raster output with this coordinate system"

        elif inputSRName == "USA_Contiguous_Albers_Equal_Area_Conic_USGS_version":
            xNLCD = 532695
            yNLCD = 1550295

        elif inputSRName == "Hawaii_Albers_Equal_Area_Conic":
            xNLCD = -29805
            yNLCD = 839235

        elif inputSRName == "NAD_1983_Alaska_Albers":
            xNLCD = -368805
            yNLCD = 1362465

        elif inputSRName == "WGS_1984_Albers":
            # New WGS 1984 based coordinate system matching USGS 2001 NLCD for Alaska
            xNLCD = -366405
            yNLCD = 2032455

        elif inputSRName == "Western_Pacific_Albers_Equal_Area_Conic":
            # WGS 1984 Albers for PAC Basin area
            xNLCD = -2390975
            yNLCD = -703265

        else:
            PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
            raise MyError, "Unable to align raster output with this coordinate system"

        pExtent = theDesc.extent
        x1 = float(pExtent.XMin)
        y1 = float(pExtent.YMin)
        x2 = float(pExtent.XMax)
        y2 = float(pExtent.YMax)

        # Round off coordinates to integer values based upon raster resolution
        # Use +- 5 meters to align with NLCD
        # Calculate snapgrid using 30 meter Kansas NLCD Lower Left coordinates = -532,695 X 1,550,295
        #
        xNLCD = 532695
        yNLCD = 1550295
        iRaster = int(iRaster)

        # Calculate number of columns difference between KS NLCD and the input extent
        # Align with NLCD CONUS
        iCol = int((x1 - xNLCD) / 30)
        iRow = int((y1 - yNLCD) / 30)
        x1 = (30 * iCol) + xNLCD - 30
        y1 = (30 * iRow) + yNLCD - 30

        numCols = int(round(abs(x2 - x1) / 30))
        numRows = int(round(abs(y2 - y1) / 30))

        x2 = numCols * 30 + x1
        y2 = numRows * 30 + y1

        theExtent = str(x1) + " " + str(y1) + " " + str(x2) + " " + str(y2)
        # Format coordinate pairs as string
        sX1 = Number_Format(x1, 0, True)
        sY1 = Number_Format(y1, 0, True)
        sX2 = Number_Format(x2, 0, True)
        sY2 = Number_Format(y2, 0, True)
        sLen = 11
        sX1 = ((sLen - len(sX1)) * " ") + sX1
        sY1 = " X " + ((sLen - len(sY1)) * " ") + sY1
        sX2 = ((sLen - len(sX2)) * " ") + sX2
        sY2 = " X " + ((sLen - len(sY2)) * " ") + sY2

        PrintMsg(" \nAligning output raster to match NLCD:", 0)
        PrintMsg("\tUR: " + sX2 + sY2 + " " + theUnits.lower(), 0)
        PrintMsg("\tLL: " + sX1 + sY1 + " " + theUnits.lower(), 0)
        PrintMsg(" \n\tNumber of rows =    \t" + Number_Format(numRows * 30 / iRaster), 0)
        PrintMsg("\tNumber of columns = \t" + Number_Format(numCols * 30 / iRaster), 0)

        return theExtent

    except MyError, e:
        # Example: raise MyError, "This is an error message"
        PrintMsg(str(e) + " \n ", 2)
        return ""

    except:
        errorMsg()
        return ""

## ===================================================================================
def WriteToLog(theMsg, theRptFile):
    # prints message to screen if run as a python script
    # Adds tool message to the geoprocessor
    #print msg
    #
    try:
        fh = open(theRptFile, "a")
        theMsg = "\n" + theMsg
        fh.write(theMsg)
        fh.close()

    except:
        errorMsg()
        pass

## ===================================================================================
def elapsedTime(start):
    # Calculate amount of time since "start" and return time string
    try:
        # Stop timer
        #
        end = time.time()

        # Calculate total elapsed seconds
        eTotal = end - start

        # day = 86400 seconds
        # hour = 3600 seconds
        # minute = 60 seconds

        eMsg = ""

        # calculate elapsed days
        eDay1 = eTotal / 86400
        eDay2 = math.modf(eDay1)
        eDay = int(eDay2[1])
        eDayR = eDay2[0]

        if eDay > 1:
          eMsg = eMsg + str(eDay) + " days "
        elif eDay == 1:
          eMsg = eMsg + str(eDay) + " day "

        # Calculated elapsed hours
        eHour1 = eDayR * 24
        eHour2 = math.modf(eHour1)
        eHour = int(eHour2[1])
        eHourR = eHour2[0]

        if eDay > 0 or eHour > 0:
            if eHour > 1:
                eMsg = eMsg + str(eHour) + " hours "
            else:
                eMsg = eMsg + str(eHour) + " hour "

        # Calculate elapsed minutes
        eMinute1 = eHourR * 60
        eMinute2 = math.modf(eMinute1)
        eMinute = int(eMinute2[1])
        eMinuteR = eMinute2[0]

        if eDay > 0 or eHour > 0 or eMinute > 0:
            if eMinute > 1:
                eMsg = eMsg + str(eMinute) + " minutes "
            else:
                eMsg = eMsg + str(eMinute) + " minute "

        # Calculate elapsed secons
        eSeconds = "%.1f" % (eMinuteR * 60)

        if eSeconds == "1.00":
            eMsg = eMsg + eSeconds + " second "
        else:
            eMsg = eMsg + eSeconds + " seconds "

        return eMsg

    except:
        errorMsg()
        return ""

## ===================================================================================
def Number_Format(num, places=0, bCommas=True):
    try:
    # Format a number according to locality and given places
        #locale.setlocale(locale.LC_ALL, "")
        if bCommas:
            theNumber = locale.format("%.*f", (places, num), True)

        else:
            theNumber = locale.format("%.*f", (places, num), False)
        return theNumber

    except:
        errorMsg()
        return False

## ===================================================================================
def ListEnv():
    # List geoprocessing environment settings
    try:
        environments = arcpy.ListEnvironments()

        # Sort the environment list, disregarding capitalization
        #
        environments.sort(key=string.lower)

        for environment in environments:
            # As the environment is passed as a variable, use Python's getattr
            #   to evaluate the environment's value
            #
            envSetting = getattr(arcpy.env, environment)

            # Format and print each environment and its current setting
            #
            #print "{0:<30}: {1}".format(environment, envSetting)
            PrintMsg("\t" + environment + ": " + str(envSetting), 0)

    except:
        errorMsg()

## ===================================================================================
def StateNames():
    # Create dictionary object containing list of state abbreviations and their names that
    # will be used to name the file geodatabase.
    # For some areas such as Puerto Rico, U.S. Virgin Islands, Pacific Islands Area the
    # abbrevation is

    # NEED TO UPDATE THIS FUNCTION TO USE THE LAOVERLAP TABLE AREANAME. AREASYMBOL IS STATE ABBREV

    try:
        stDict = dict()
        stDict["AL"] = "Alabama"
        stDict["AK"] = "Alaska"
        stDict["AS"] = "American Samoa"
        stDict["AZ"] = "Arizona"
        stDict["AR"] = "Arkansas"
        stDict["CA"] = "California"
        stDict["CO"] = "Colorado"
        stDict["CT"] = "Connecticut"
        stDict["DC"] = "District of Columbia"
        stDict["DE"] = "Delaware"
        stDict["FL"] = "Florida"
        stDict["GA"] = "Georgia"
        stDict["HI"] = "Hawaii"
        stDict["ID"] = "Idaho"
        stDict["IL"] = "Illinois"
        stDict["IN"] = "Indiana"
        stDict["IA"] = "Iowa"
        stDict["KS"] = "Kansas"
        stDict["KY"] = "Kentucky"
        stDict["LA"] = "Louisiana"
        stDict["ME"] = "Maine"
        stDict["MD"] = "Maryland"
        stDict["MA"] = "Massachusetts"
        stDict["MI"] = "Michigan"
        stDict["MN"] = "Minnesota"
        stDict["MS"] = "Mississippi"
        stDict["MO"] = "Missouri"
        stDict["MT"] = "Montana"
        stDict["NE"] = "Nebraska"
        stDict["NV"] = "Nevada"
        stDict["NH"] = "New Hampshire"
        stDict["NJ"] = "New Jersey"
        stDict["NM"] = "New Mexico"
        stDict["NY"] = "New York"
        stDict["NC"] = "North Carolina"
        stDict["ND"] = "North Dakota"
        stDict["OH"] = "Ohio"
        stDict["OK"] = "Oklahoma"
        stDict["OR"] = "Oregon"
        stDict["PA"] = "Pennsylvania"
        stDict["PRUSVI"] = "Puerto Rico and U.S. Virgin Islands"
        stDict["RI"] = "Rhode Island"
        stDict["Sc"] = "South Carolina"
        stDict["SD"] ="South Dakota"
        stDict["TN"] = "Tennessee"
        stDict["TX"] = "Texas"
        stDict["UT"] = "Utah"
        stDict["VT"] = "Vermont"
        stDict["VA"] = "Virginia"
        stDict["WA"] = "Washington"
        stDict["WV"] = "West Virginia"
        stDict["WI"] = "Wisconsin"
        stDict["WY"] = "Wyoming"
        return stDict

    except:
        PrintMsg("\tFailed to create list of state abbreviations (CreateStateList)", 2)
        return stDict

## ===================================================================================
def CheckStatistics(outputRaster):
    # For no apparent reason, sometimes ArcGIS fails to build statistics. Might work one
    # time and then the next time it may fail without any error message.
    #
    try:
        #PrintMsg(" \n\tChecking raster statistics", 0)

        for propType in ['MINIMUM', 'MAXIMUM', 'MEAN', 'STD']:
            statVal = arcpy.GetRasterProperties_management (outputRaster, propType).getOutput(0)
            #PrintMsg("\t\t" + propType + ": " + statVal, 1)

        return True

    except:
        return False

## ===================================================================================
def UpdateMetadata(outputWS, target, surveyInfo, iRaster):
    #
    # Used for non-ISO metadata
    #
    # Search words:  xxSTATExx, xxSURVEYSxx, xxTODAYxx, xxFYxx
    #
    try:
        PrintMsg("\tUpdating metadata...")
        arcpy.SetProgressor("default", "Updating metadata")

        # Set metadata translator file
        dInstall = arcpy.GetInstallInfo()
        installPath = dInstall["InstallDir"]
        prod = r"Metadata/Translator/ARCGIS2FGDC.xml"
        mdTranslator = os.path.join(installPath, prod)

        # Define input and output XML files
        mdImport = os.path.join(env.scratchFolder, "xxImport.xml")  # the metadata xml that will provide the updated info
        xmlPath = os.path.dirname(sys.argv[0])
        mdExport = os.path.join(xmlPath, "gSSURGO_MapunitRaster.xml") # original template metadata in script directory

        # Cleanup output XML files from previous runs
        if os.path.isfile(mdImport):
            os.remove(mdImport)

        # Get replacement value for the search words
        #
        stDict = StateNames()
        st = os.path.basename(outputWS)[8:-4]

        if st in stDict:
            # Get state name from the geodatabase
            mdState = stDict[st]

        else:
            # Leave state name blank. In the future it would be nice to include a tile name when appropriate
            mdState = ""

        # Set date strings for metadata, based upon today's date
        #
        d = datetime.date.today()
        today = str(d.isoformat().replace("-",""))

        # Set fiscal year according to the current month. If run during January thru September,
        # set it to the current calendar year. Otherwise set it to the next calendar year.
        #
        if d.month > 9:
            fy = "FY" + str(d.year + 1)

        else:
            fy = "FY" + str(d.year)

        # Convert XML to tree format
        tree = ET.parse(mdExport)
        root = tree.getroot()

        # new citeInfo has title.text, edition.text, serinfo/issue.text
        citeInfo = root.findall('idinfo/citation/citeinfo/')

        if not citeInfo is None:
            # Process citation elements
            # title, edition, issue
            #
            for child in citeInfo:
                #PrintMsg("\t\t" + str(child.tag), 0)
                if child.tag == "title":
                    if child.text.find('xxSTATExx') >= 0:
                        child.text = child.text.replace('xxSTATExx', mdState)

                    elif mdState != "":
                        child.text = child.text + " - " + mdState

                elif child.tag == "edition":
                    if child.text == 'xxFYxx':
                        child.text = fy

                elif child.tag == "serinfo":
                    for subchild in child.iter('issue'):
                        if subchild.text == "xxFYxx":
                            subchild.text = fy

        # Update place keywords
        ePlace = root.find('idinfo/keywords/place')

        if not ePlace is None:
            #PrintMsg("\t\tplace keywords", 0)

            for child in ePlace.iter('placekey'):
                if child.text == "xxSTATExx":
                    child.text = mdState

                elif child.text == "xxSURVEYSxx":
                    child.text = surveyInfo

        # Update credits
        eIdInfo = root.find('idinfo')
        if not eIdInfo is None:
            #PrintMsg("\t\tcredits", 0)

            for child in eIdInfo.iter('datacred'):
                sCreds = child.text

                if sCreds.find("xxSTATExx") >= 0:
                    #PrintMsg("\t\tcredits " + mdState, 0)
                    child.text = child.text.replace("xxSTATExx", mdState)

                if sCreds.find("xxFYxx") >= 0:
                    #PrintMsg("\t\tcredits " + fy, 0)
                    child.text = child.text.replace("xxFYxx", fy)

                if sCreds.find("xxTODAYxx") >= 0:
                    #PrintMsg("\t\tcredits " + today, 0)
                    child.text = child.text.replace("xxTODAYxx", today)

        idPurpose = root.find('idinfo/descript/purpose')
        if not idPurpose is None:
            ip = idPurpose.text

            if ip.find("xxFYxx") >= 0:
                idPurpose.text = ip.replace("xxFYxx", fy)
                #PrintMsg("\t\tpurpose", 0)

        #  create new xml file which will be imported, thereby updating the table's metadata
        tree.write(mdImport, encoding="utf-8", xml_declaration=None, default_namespace=None, method="xml")

        # import updated metadata to the geodatabase table
        # Using three different methods with the same XML file works for ArcGIS 10.1
        #
        #PrintMsg("\t\tApplying metadata translators...")
        arcpy.MetadataImporter_conversion (mdImport, target)
        arcpy.ImportMetadata_conversion(mdImport, "FROM_FGDC", target, "DISABLED")
        #arcpy.ImportMetadata_conversion(mdImport, "FROM_FGDC", target, "DISABLED")

        # delete the temporary xml metadata file
        if os.path.isfile(mdImport):
            os.remove(mdImport)
            pass

        # delete metadata tool logs
        logFolder = os.path.dirname(env.scratchFolder)
        logFile = os.path.basename(mdImport).split(".")[0] + "*"


        currentWS = env.workspace
        env.workspace = logFolder
        logList = arcpy.ListFiles(logFile)

        for lg in logList:
            arcpy.Delete_management(lg)

        env.workspace = currentWS

        return True

    except:
        errorMsg()
        False

## ===================================================================================
def GetExtent(tmpPolys, iRaster):
    # Set output extent to that of the temporary 'tile' featurelayer, snapping
    # it to the appropriate NLCD coordinates
    try:

        # User searchcursor on polygon geometry to get extent
        #
        # Set default extent values
        xMin =  1000000000
        yMin =  1000000000
        xMax = -1000000000
        yMax = -1000000000

        # Assuming that input cooordinate system is Albers
        with arcpy.da.SearchCursor(tmpPolys, ["SHAPE@"]) as srcCursor:
            for row in srcCursor:
                pExtent = row[0].extent
                a1 = pExtent.XMin
                b1 = pExtent.YMin
                a2 = pExtent.XMax
                b2 = pExtent.YMax

                if a1 < xMin:
                    xMin = a1

                if b1 < yMin:
                    yMin = b1

                if a2 > xMax:
                    xMax = a2

                if b2 > yMax:
                    yMax = b2


        x1 = float(xMin)
        y1 = float(yMin)
        x2 = float(xMax)
        y2 = float(yMax)

        sExtent = str(x1) + " " + str(y1) + " " + str(x2) + " " + str(y2)

        # Begin SnapToNLCD Code
        theDesc = arcpy.Describe(tmpPolys)
        sr = theDesc.spatialReference
        inputSRName = sr.name
        theUnits = sr.linearUnitName

        if 'foot' in theUnits.lower():
            theUnits = "feet"

        elif theUnits.lower() == "meter":
            theUnits = "meters"

        # USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (NAD83)
        #xNLCD = 532695
        #yNLCD = 1550295

        # Hawaii_Albers_Equal_Area_Conic  -345945, 1753875
        # Western_Pacific_Albers_Equal_Area_Conic  -2390975, -703265 est.
        # NAD_1983_Alaska_Albers  -2232345, 344805
        # WGS_1984_Alaska_Albers  Upper Left Corner:  -366405.000000 meters(X),  2380125.000000 meters(Y)
        # WGS_1984_Alaska_Albers  Lower Right Corner: 517425.000000 meters(X),  2032455.000000 meters(Y)
        # Puerto Rico 3092415, -78975 (CONUS works for both)

        if theUnits != "meters":
            PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
            raise MyError, "Unable to align raster output with this coordinate system"

        elif inputSRName == "USA_Contiguous_Albers_Equal_Area_Conic_USGS_version":
            xNLCD = 532695
            yNLCD = 1550295

        elif inputSRName == "Hawaii_Albers_Equal_Area_Conic":
            xNLCD = -29805
            yNLCD = 839235

        elif inputSRName == "NAD_1983_Alaska_Albers":
            xNLCD = -368805
            yNLCD = 1362465

        elif inputSRName == "WGS_1984_Albers":
            # New WGS 1984 based coordinate system matching USGS 2001 NLCD for Alaska
            xNLCD = -366405
            yNLCD = 2032455

        elif inputSRName == "Western_Pacific_Albers_Equal_Area_Conic":
            # WGS 1984 Albers for PAC Basin area
            xNLCD = -2390975
            yNLCD = -703265

        else:
            PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
            raise MyError, "Unable to align raster output with this coordinate system"

        # Calculate number of columns difference between KS NLCD and the input extent
        # Align with NLCD CONUS
        iCol = int((x1 - xNLCD) / 30)
        iRow = int((y1 - yNLCD) / 30)

        x1 = (30 * iCol) + xNLCD - 30
        y1 = (30 * iRow) + yNLCD - 30

        numCols = int(round(abs(x2 - x1) / 30))
        numRows = int(round(abs(y2 - y1) / 30))

        x2 = numCols * 30 + x1
        y2 = numRows * 30 + y1

        theExtent = str(x1) + " " + str(y1) + " " + str(x2) + " " + str(y2)
        # Format coordinate pairs as string
        sX1 = Number_Format(x1, 0, True)
        sY1 = Number_Format(y1, 0, True)
        sX2 = Number_Format(x2, 0, True)
        sY2 = Number_Format(y2, 0, True)
        sLen = 11
        sX1 = ((sLen - len(sX1)) * " ") + sX1
        sY1 = " X " + ((sLen - len(sY1)) * " ") + sY1
        sX2 = ((sLen - len(sX2)) * " ") + sX2
        sY2 = " X " + ((sLen - len(sY2)) * " ") + sY2

        env.extent = theExtent
        return theExtent
        # End of SnapToNLCD code

    except MyError, e:
        # Example: raise MyError, "This is an error message"
        PrintMsg(str(e), 2)
        return ""

    except:
        errorMsg()
        return ""

## ===================================================================================
def CheckSpatialReference(inputFC):
    # Make sure that the coordinate system is projected and units are meters
    try:
        desc = arcpy.Describe(inputFC)
        inputSR = desc.spatialReference

        if inputSR.type.upper() == "PROJECTED":
            if inputSR.linearUnitName.upper() == "METER":
                env.outputCoordinateSystem = inputSR
                return True

            else:
                raise MyError, os.path.basename(theGDB) + ": Input soil polygon layer does not have a valid coordinate system for gSSURGO"

        else:
            raise MyError, os.path.basename(theGDB) + ": Input soil polygon layer must have a projected coordinate system"

    except MyError, e:
        # Example: raise MyError, "This is an error message"
        PrintMsg(str(e), 2)
        return False

    except:
        errorMsg()
        return False

## ===================================================================================
def ConvertToRaster(theGDB, outputWS, theSnapRaster, iRaster, bTiled):
    # main function
    try:
        # Set geoprocessing environment
        #
        #bScratch = setScratchWorkspace()
        env.overwriteOutput = True
        env.workspace = theGDB
        arcpy.CheckOutExtension("Spatial")
        arcpy.env.compression = "LZ77"
        import shutil

        # Make sure that the env.scratchGDB is NOT Default.gdb. This causes problems for
        # some unknown reason.
        if os.path.basename(env.scratchGDB).lower() == "default.gdb" or os.path.basename(env.scratchWorkspace).lower() == "default.gdb":
            raise MyError, "Invalid scratch workspace setting (" + env.scratchWorkspace + ")"

        # turn off automatic Pyramid creation and Statistics calculation
        env.rasterStatistics = "NONE"
        env.pyramid = "PYRAMIDS 0"

        if bTiled:
            # Try creating a temporary folder for holding temporary rasters
            # This may allow the entire folder to be deleted at once instead of one raster at-a-time
            tmpFolder = os.path.join(env.scratchFolder, "TmpRasters")

            if arcpy.Exists(tmpFolder):
                shutil.rmtree(tmpFolder)

        if outputWS != "":
            # determine type of output workspace and use it to set the output raster format
            if outputWS[-4:].lower() in ['.gdb','.mdb']:
                rasterExt = ""

            else:
                rasterExt = ".img"

            if arcpy.CheckExtension("Spatial") == "Available":
                arcpy.CheckOutExtension("Spatial")
                tileName = os.path.basename(theGDB)[os.path.basename(theGDB).rfind("_") + 1:os.path.basename(theGDB).rfind(".gdb")]
                outputRaster = os.path.join(outputWS, "MapunitRaster_" + tileName + "_" + str(iRaster) + "m" + rasterExt)

            else:
                raise MyError, "Required Spatial Analyst extension is not available"

        else:
            # create final raster in same workspace as input polygon featureclass
            outputWS = theGDB
            tileName = os.path.basename(theGDB)[os.path.basename(theGDB).rfind("_") + 1:os.path.basename(theGDB).rfind(".gdb")]
            outputRaster = os.path.join(theGDB, "MapunitRaster_"    + tileName + "_" + str(iRaster) + "m")

        inputFC = os.path.join(theGDB, "MUPOLYGON")

        if not arcpy.Exists(inputFC):
            raise MyError, "Could not find input featureclass: " + inputFC

        # Check input layer's coordinate system to make sure horizontal units are meters
        # set the output coordinate system for the raster (neccessary for PolygonToRaster)
        if CheckSpatialReference(inputFC) == False:
            return False

        # For rasters named using an attribute value, some attribute characters can result in
        # 'illegal' names.
        outputRaster = outputRaster.replace("-", "")

        if arcpy.Exists(outputRaster):
            arcpy.Delete_management(outputRaster)
            time.sleep(1)

        if arcpy.Exists(outputRaster):
            err = "Output raster (" + os.path.basename(outputRaster) + ") already exists"
            raise MyError, err

        start = time.time()   # start clock to measure total processing time
        #begin = time.time()   # start clock to measure set up time
        time.sleep(2)

        PrintMsg(" \nBeginning raster conversion process", 0)

        # Create Lookup table for storing MUKEY values
        #
        if bTiled:
            lu = os.path.join(theGDB, "Lookup")

        else:
            lu = os.path.join(env.scratchGDB, "Lookup")

        if arcpy.Exists(lu):
            arcpy.Delete_management(lu)

        # worked
        arcpy.CreateTable_management(os.path.dirname(lu), "Lookup")
        arcpy.AddField_management(lu, "CELLVALUE", "LONG")
        arcpy.AddField_management(lu, "MUKEY", "TEXT", "#", "#", "30")

        # Create list of MUKEY values from the MUPOLYGON featureclass
        #
        if not bTiled:
            # Create a list of map unit keys present in the MUPOLYGON featureclass
            #
            PrintMsg("\tGetting feature count and list of mukeys from input soil polygon layer...", 0)
            arcpy.SetProgressor("default", "Getting inventory of map units...")
            tmpPolys = "SoilPolygons"

            with arcpy.da.SearchCursor(inputFC, ("MUKEY",)) as srcCursor:
                # Create a unique, sorted list of MUKEY values in the MUPOLYGON featureclass
                mukeyList = [row[0] for row in srcCursor]
                polyCnt = len(mukeyList)
                mukeySet = set(mukeyList)
                mukeyList = sorted(list(mukeySet))
                del mukeySet

            if len(mukeyList) == 0:
                raise MyError, "Failed to get MUKEY values from " + inputFC

        # Create list of map unit keys AND areasymbols present in the MUPOLYGON featureclass
        #
        if bTiled:
            PrintMsg("\tGetting feature count and list of mukeys from input soil polygon layer...", 0)
            arcpy.SetProgressor("default", "Creating map unit inventory...")
            tileList = list()
            mukeyList = list()

            with arcpy.da.SearchCursor(inputFC, ["MUKEY", "AREASYMBOL"]) as cur:
                # Create a unique, sorted list of MUKEY values in the MUPOLYGON featureclass

                for rec in cur:
                    # get mukey
                    mukeyList.append(rec[0])

                    # get areasymbol for tile value
                    if not rec[1] in tileList:
                        tileList.append(rec[1])

                # create sorted and unique list of mukeys
                polyCnt = len(mukeyList)
                mukeySet = set(mukeyList)
                mukeyList = sorted(list(mukeySet))
                del mukeySet

            if len(mukeyList) == 0:
                raise MyError, "Failed to get MUKEY values from " + inputFC

            if len(tileList) == 1:
                # Only one tile in list
                # Switch to the 'un-tiled' mode
                bTiled = 0

        muCnt = len(mukeyList)

        # Load MUKEY values into Lookup table
        #
        PrintMsg("\tSaving " + Number_Format(muCnt, 0, True) + " MUKEY values for " + Number_Format(polyCnt, 0, True) + " polygons"  , 0)
        arcpy.SetProgressorLabel("Creating lookup table...")

        with arcpy.da.InsertCursor(lu, ("CELLVALUE", "MUKEY") ) as inCursor:
            for mukey in mukeyList:
                rec = mukey, mukey
                inCursor.insertRow(rec)

        # Add attribute index here.
        arcpy.AddIndex_management(lu, ["mukey"], "Indx_LU")
        #
        # End of Lookup table code

        # Match NLCD raster (snapraster)
        if theSnapRaster == "":
            fullExtent = SnapToNLCD(inputFC, iRaster)

            if fullExtent == "":
                err = "Please specify a Snap Raster"
                raise MyError, err

        else:
            PrintMsg(" \nSetting snap raster environment to: " + os.path.basename(theSnapRaster), 0)
            arcpy.snapRaster = theSnapRaster

        # Raster conversion process...
        #
        if bTiled:
            # Tiled raster process...
            #
            # Reset stopwatch for measuring the raster tile conversion plus mosaic

            #theMsg = "\tRaster conversion setup: " + elapsedTime(begin)
            #PrintMsg(theMsg, 1)

            #begin = time.time()
            PrintMsg(" \n\tCreating " + Number_Format(len(tileList), 0, True) + " raster tiles from " + os.path.join(os.path.basename(os.path.dirname(inputFC)), os.path.basename(inputFC)) + " featureclass", 0)

            # Create output folder
            arcpy.CreateFolder_management(env.scratchFolder, os.path.basename(tmpFolder))
            rasterList = list()
            i = 0

            for tile in tileList:
                i += 1
                tmpPolys = "poly_" + tile
                #tileRaster = os.path.join(env.scratchGDB, "r" + tile.lower())
                tileRaster = os.path.join(tmpFolder, tile.lower() + ".tif")
                rasterList.append(tileRaster)
                arcpy.SetProgressor("default", "Performing raster conversion for " + tile + "  (tile " + str(i) + " of " + str(len(tileList)) + ")...")
                wc = "AREASYMBOL = '" + tile + "'"
                arcpy.MakeFeatureLayer_management (inputFC, tmpPolys, wc)
                tileExtent = GetExtent(tmpPolys, iRaster)

                if tileExtent == "":
                    raise MyError, "Bad extent for " + tile

                arcpy.AddJoin_management (tmpPolys, "MUKEY", lu, "MUKEY", "KEEP_ALL")
                arcpy.PolygonToRaster_conversion(tmpPolys, "Lookup.CELLVALUE", tileRaster, "MAXIMUM_COMBINED_AREA", "", iRaster)
                arcpy.Delete_management(tmpPolys)
                del tmpPolys

            del tileRaster
            PrintMsg("\tMosaicing tiles to a single raster...", 0)
            env.extent = fullExtent
            arcpy.MosaicToNewRaster_management (rasterList, os.path.dirname(outputRaster), os.path.basename(outputRaster), "", "32_BIT_UNSIGNED", iRaster, 1, "MAXIMUM")

            # Cleanup the temporary tiled rasters
            # Should probably compact the scratchGDB afterwards...
            #
            arcpy.SetProgressor("step", "Cleaning up temporary raster tiles...", 1, len(rasterList), 1)
            #PrintMsg("\tCleaning up raster tiles...", 0)

            #for ras in rasterList:
            #    arcpy.Delete_management(ras)
            #    arcpy.SetProgressorPosition()


            del rasterList
            # Compact the scratch geodatabase after deleting all the rasters
            arcpy.Compact_management(env.scratchGDB)
            #theMsg = "\tPolygonToRaster conversion and mosaic: " + elapsedTime(begin)
            #PrintMsg(theMsg, 1)

        else:
            # Create a single raster, no tiles
            #
            # Reset stopwatch for measuring just the single raster conversion
            #begin = time.time()
            #PrintMsg(" \n" + (65 * "*"), 0)
            PrintMsg(" \nConverting featureclass " + os.path.join(os.path.basename(os.path.dirname(inputFC)), os.path.basename(inputFC)) + " to raster (" + str(iRaster) + " meter)", 0)
            #PrintMsg(" \n" + (65 * "*"), 0)
            tmpPolys = "poly_tmp"
            arcpy.MakeFeatureLayer_management (inputFC, tmpPolys)
            arcpy.AddJoin_management (tmpPolys, "MUKEY", lu, "MUKEY", "KEEP_ALL")

            arcpy.SetProgressor("default", "Running PolygonToRaster conversion...")
            #theMsg = "\tRaster conversion setup: " + elapsedTime(begin)
            #PrintMsg(theMsg, 1)

            env.extent = fullExtent
            arcpy.PolygonToRaster_conversion(tmpPolys, "Lookup.CELLVALUE", outputRaster, "MAXIMUM_COMBINED_AREA", "", iRaster)
            #theMsg = " \n\tPolygonToRaster conversion: " + elapsedTime(begin)
            #PrintMsg(theMsg, 1)

            # immediately delete temporary polygon layer to free up memory for the rest of the process
            arcpy.Delete_management(tmpPolys)

            # End of single raster process

        # Now finish up the single temporary raster
        #
        PrintMsg(" \nCompleting conversion process:", 0)
        # Reset the stopwatch for the raster post-processing
        #begin = time.time()

        # Remove lookup table
        if arcpy.Exists(lu):
            arcpy.Delete_management(lu)

        # ****************************************************
        # Build pyramids and statistics
        # ****************************************************
        if arcpy.Exists(outputRaster):
        #if arcpy.Exists(tmpRaster):
            time.sleep(3)
            arcpy.SetProgressor("default", "Calculating raster statistics...")
            PrintMsg("\tCalculating raster statistics...", 0)
            env.pyramid = "PYRAMIDS -1 NEAREST"
            arcpy.env.rasterStatistics = 'STATISTICS 100 100'
            #arcpy.BuildPyramidsandStatistics_management(os.path.dirname(outputRaster), "NONE", "BUILD_PYRAMIDS", "CALCULATE_STATISTICS")
            arcpy.CalculateStatistics_management (outputRaster, 1, 1, "", "OVERWRITE" )

            if CheckStatistics(outputRaster) == False:
                # For some reason the BuildPyramidsandStatistics command failed to build statistics for this raster.
                #
                # Try using CalculateStatistics while setting an AOI
                PrintMsg("\tInitial attempt to create statistics failed, trying another method...", 0)
                time.sleep(3)

                if arcpy.Exists(os.path.join(outputWS, "SAPOLYGON")):
                    # Try running CalculateStatistics with an AOI to limit the area that is processed
                    # if we have to use SAPOLYGON as an AOI, this will be REALLY slow
                    #arcpy.CalculateStatistics_management (outputRaster, 1, 1, "", "OVERWRITE", os.path.join(outputWS, "SAPOLYGON") )
                    arcpy.CalculateStatistics_management (outputRaster, 1, 1, "", "OVERWRITE" )

                if CheckStatistics(outputRaster) == False:
                    time.sleep(3)
                    PrintMsg("\tFailed in both attempts to create statistics for raster layer", 1)

            arcpy.SetProgressor("default", "Building pyramids...")
            PrintMsg("\tBuilding pyramids...", 0)
            arcpy.BuildPyramids_management(outputRaster, "-1", "NONE", "NEAREST", "DEFAULT", "", "SKIP_EXISTING")

            # ****************************************************
            # Add MUKEY to final raster
            # ****************************************************
            # Build attribute table for final output raster. Sometimes it fails to automatically build.
            PrintMsg("\tBuilding raster attribute table and updating MUKEY values", )
            arcpy.SetProgressor("default", "Building raster attrribute table...")
            arcpy.BuildRasterAttributeTable_management(outputRaster)

            # Add MUKEY values to final mapunit raster
            #
            arcpy.SetProgressor("default", "Adding MUKEY attribute to raster...")
            arcpy.AddField_management(outputRaster, "MUKEY", "TEXT", "#", "#", "30")
            with arcpy.da.UpdateCursor(outputRaster, ["VALUE", "MUKEY"]) as cur:
                for rec in cur:
                    rec[1] = rec[0]
                    cur.updateRow(rec)

            # Add attribute index (MUKEY) for raster
            arcpy.AddIndex_management(outputRaster, ["mukey"], "Indx_RasterMukey")

        else:
            err = "Missing output raster (" + outputRaster + ")"
            raise MyError, err

        # Compare list of original mukeys with the list of raster mukeys
        # Report discrepancies. These are usually thin polygons along survey boundaries,
        # added to facilitate a line-join.
        #
        arcpy.SetProgressor("default", "Looking for missing map units...")
        rCnt = int(arcpy.GetRasterProperties_management (outputRaster, "UNIQUEVALUECOUNT").getOutput(0))

        if rCnt <> muCnt:
            missingList = list()
            rList = list()

            # Create list of raster mukeys...
            with arcpy.da.SearchCursor(outputRaster, ("MUKEY",)) as rcur:
                for rec in rcur:
                    mukey = rec[0]
                    rList.append(mukey)

            missingList = list(set(mukeyList) - set(rList))
            queryList = list()
            for mukey in missingList:
                queryList.append("'" + mukey + "'")

            PrintMsg("\tDiscrepancy in mapunit count for new raster", 1)
            PrintMsg("\t\tInput polygon mapunits: " + Number_Format(muCnt, 0, True), 0)
            PrintMsg("\t\tOutput raster mapunits: " + Number_Format(rCnt, 0, True), 0)
            PrintMsg("\t\tMUKEY IN (" + ", ".join(queryList) + ") \n ", 0)

        # Update metadata file for the geodatabase
        #
        # Query the output SACATALOG table to get list of surveys that were exported to the gSSURGO
        #
        saTbl = os.path.join(theGDB, "sacatalog")
        expList = list()

        with arcpy.da.SearchCursor(saTbl, ("AREASYMBOL", "SAVEREST")) as srcCursor:
            for rec in srcCursor:
                expList.append(rec[0] + " (" + str(rec[1]).split()[0] + ")")

        surveyInfo = ", ".join(expList)
        time.sleep(2)
        arcpy.SetProgressorLabel("Updating metadata...")

        bMetaData = UpdateMetadata(outputWS, outputRaster, surveyInfo, iRaster)

        arcpy.SetProgressorLabel("Compacting metadata...")
        PrintMsg("\tCompacting geodatabase...", 0)
        arcpy.Compact_management(theGDB)
        arcpy.RefreshCatalog(os.path.dirname(outputRaster))

        if bTiled:

            if arcpy.Exists(tmpFolder):
                PrintMsg("\tCleaning up raster tiles...", 0)
                shutil.rmtree(tmpFolder)

        #theMsg = "\tPost-processing: " + elapsedTime(begin)
        #PrintMsg(theMsg, 1)

        theMsg = " \nProcessing time for " + outputRaster + ": " + elapsedTime(start) + " \n "
        PrintMsg(theMsg, 0)

        del outputRaster
        del inputFC

        #WriteToLog(theMsg, theRptFile)
        arcpy.CheckInExtension("Spatial")

        return True

    except MyError, e:
        # Example: raise MyError, "This is an error message"
        PrintMsg(str(e), 2)
        return False
        arcpy.CheckInExtension("Spatial")

    except:
        errorMsg()
        return False
        arcpy.CheckInExtension("Spatial")

## ===================================================================================
## ===================================================================================
## MAIN
## ===================================================================================

# Import system modules
import sys, string, os, arcpy, locale, traceback, math, time, datetime
import xml.etree.cElementTree as ET
from arcpy import env
from arcpy.sa import *

# Create the Geoprocessor object
try:
    if __name__ == "__main__":
        # get parameters
        theGDB = arcpy.GetParameterAsText(0)                   # required geodatabase containing MUPOLYGON featureclass
        outputWS = arcpy.GetParameterAsText(1)                # optional output workspace. If not set, output raster will be created in the same workspace
        theSnapRaster = arcpy.GetParameterAsText(2)           # optional snap raster. If not set, uses NLCD Albers USGS NAD83 for CONUS
        iRaster = arcpy.GetParameter(3)                       # output raster resolution
        bTiled = arcpy.GetParameter(4)                                         # boolean - split raster into survey-tiles and then mosaic

        env.overwriteOutput= True
        # Call function that does all of the work
        bRaster = ConvertToRaster(theGDB, outputWS, theSnapRaster, iRaster, bTiled)

except MyError, e:
    # Example: raise MyError, "This is an error message"
    PrintMsg(str(e), 2)

except:
    errorMsg()
