# ------------------------------------------------------------------------------
#
# Copyright 2011, 2012, 2013 Brent L. Brock and the Craighead Institute
#
# This file is part of Wild Planner.
#
# Wild Planner is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Wild Planner is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Wild Planner in the file named LICENSE.TXT.  If not, see <http://www.gnu.org/licenses/>.
#
# ------------------------------------------------------------------------------
#
# evaluate_corridor.py
# Created on: Wed Dec 24 2008 10:35:19 AM
#   (generated by ArcGIS/ModelBuilder)
# Usage: evaluate_corridor <Structures_090208> 
# ---------------------------------------------------------------------------

# Import system modules
import sys, os, arcpy #arcgisscripting
import functionlib as fl

# Create the Geoprocessor object
#gp = arcgisscripting.create()

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

# Make sure outputs can be overwritten
arcpy.env.overwriteOutputs = 1

# Script arguments...
pointLayer = sys.argv[1]
minCorridorWidth = sys.argv[2]
distDistance = sys.argv[3]
outRaster = sys.argv[4]
cellSize = sys.argv[5]
extractMask = sys.argv[6]
aExtent = sys.argv[7]
viewShed = sys.argv[8]
elevLayer = sys.argv[9]

if aExtent == "" or aExtent == "#":
    if extractMask == "" or extractMask == "#":
        arcpy.AddMessage('Setting analysis extent to ' + pointLayer + '...')
        aExtent = arcpy.Describe(pointLayer).Extent
    else:
        arcpy.AddMessage('Setting analysis extent to ' + extractMask + '...')
        aExtent = arcpy.Describe(extractMask).Extent       
else:
    arcpy.AddMessage('Setting analysis extent to user specified extent...')
    
arcpy.env.extent = aExtent

# Run Repair Geometry on pointLayer

arcpy.RepairGeometry_management (pointLayer)


#n = os.path.dirname(outRaster)
desc = arcpy.Describe(os.path.dirname(outRaster))
if not desc.workspaceType == u'FileSystem':
    if len(os.path.splitext(outRaster)[1]) > 0:
        arcpy.AddWarning("Output grid format is not supported for writing to geodatabase. \nRemoving extension from output name...")
    outRaster = os.path.splitext(outRaster)[0]
    if len(os.path.basename(outRaster)) > 13:
        arcpy.AddWarning("ESRI grid name exceeds 13 charcters. Truncating output raster name...")
        bn = os.path.basename(outRaster)[:13]
        outRaster = os.path.join(os.path.dirname(outRaster), bn)
        arcpy.AddWarning("Output will be written to: " + outRaster)
        


                      
# Local variables...

# Process distance inputs...
arcpy.AddMessage("Converting input units to meters...")
width = minCorridorWidth.split(" ")
minCorridorWidth = str(fl.ConvertDistanceToMeters(width[0], width[1]))
del width

distance = distDistance.split(" ")
distDistance = str(fl.ConvertDistanceToMeters(distance[0], distance[1]))
del distance

cell = cellSize.split(" ")
if cellSize == "" or cellSize == "#":
    cellSize = round(float(minCorridorWidth) / 30)
    arcpy.AddMessage("\tCell size set to " + str(cellSize) + " meters...")
else:
    cellSize = str(fl.ConvertDistanceToMeters(cell[0], cell[1]))
del cell

radius = float(minCorridorWidth)/2
InNeighborhood = "CIRCLE " + str(radius) + " MAP"

# Process: Euclidean Distance...
arcpy.AddMessage("Calculating euclidian distance from input structures...")
Euc_Dist = arcpy.sa.EucDistance(pointLayer, "", cellSize, "")
arcpy.AddMessage("Identifying disturbance areas...")
# Reclassify to define disturbed areas...
reclassRaster = arcpy.sa.Reclassify(Euc_Dist, "Value", "0 " + distDistance + " 1", "NODATA")
arcpy.AddMessage("Identifying potential travel areas...")

# If extraction option selected, process extract mask...
if not extractMask == "" and not extractMask == "#":
    arcpy.AddMessage("Extract mask selected. Processing " + arcpy.Describe(extractMask).baseName + " for corridor calculations...")
    dType = arcpy.Describe(extractMask).dataType
    extractProjName = arcpy.Describe(extractMask).spatialreference.name
    outProjName = arcpy.Describe(pointLayer).spatialreference.name

    if not extractProjName == outProjName:
        arcpy.AddMessage("\tReprojecting extractMask layer...")
        if dType == "ShapeFile" or dType == "FeatureLayer":
            arcpy.AddMessage("extractMask: " + extractMask)
            arcpy.AddMessage("pointLayer: " + pointLayer)
            output = arcpy.Project_management (extractMask, arcpy.env.scratchFolder + os.path.sep + "reproj.shp", pointLayer)
        else:
            output = arcpy.ProjectRaster_management (extractMask, "in_memory\\reproj", pointLayer)
    if dType == "RasterLayer":
        arcpy.AddMessage("\tConverting extract mask raster to polygons...")
        extractPoly = arcpy.RasterToPolygon_conversion (extractMask, "in_memory\\extractPoly", "NO_SIMPLIFY")
    else:
        extractPoly = extractMask

#Convert polygon boundaries to lines...
    arcpy.AddMessage("\tExtracting mask polygon edges to polylines...")
    extractLine = arcpy.FeatureToLine_management (extractPoly, "in_memory\\extractLine")
    
#Convert polygon boundary lines to raster...
    arcpy.AddMessage("\tConverting extract mask polylines to raster...")
    arcpy.env.snapRaster = arcpy.Describe(reclassRaster).catalogPath
    extractEdge = arcpy.FeatureToRaster_conversion (extractLine, "FID", "in_memory\\extractEdge", cellSize)
    arcpy.env.snapRaster = ""
    MAX = arcpy.GetRasterProperties_management (extractEdge, "Maximum")
    extractEdge = arcpy.sa.Reclassify(extractEdge, "Value", "0 " + str(MAX) + " 1 ; NODATA 0", "NODATA")
    arcpy.AddMessage("\tCombining extraction polygons with disturbance areas...")
    reclassRaster = arcpy.sa.Con(reclassRaster, 1, extractEdge, "Value = 1")

# Calculate euclidian distance from disturbed areas...
Euc_Dist2 = arcpy.sa.EucDistance(reclassRaster, "", cellSize, "")

# Get focal max within radius equal one half minimum corridor width...
maxRaster = arcpy.sa.FocalStatistics(Euc_Dist2, InNeighborhood, "MAXIMUM", "DATA")
MAX = arcpy.GetRasterProperties_management (Euc_Dist2, "Maximum")
# Reclassify final output. Movement area has focal max >= 1/2 min corridor width...
Reclassification = "0 " + str(radius) + " 0;" + str(radius) + " " + str(MAX) + " 1"
output = arcpy.sa.Reclassify(maxRaster, "VALUE", Reclassification, "DATA")

if viewShed == 'true':
    arcpy.AddMessage("Initiating viewshed analysis for travel areas...")
##    vName = outWorkspace + os.sep + raster[1][:7]
#    outVShed = tWorkspace + os.sep + "xxvshd"
    arcpy.AddMessage("\tFinding areas hidden from view...")
    outVShed = arcpy.sa.Viewshed(elevLayer, pointLayer)
    output = arcpy.sa.Con(outVShed, "1", output, "Value = 0")

if not extractMask == "" and not extractMask == "#":
    arcpy.AddMessage('Extracting output raster by ' + extractMask + '...')
    output = arcpy.sa.ExtractByMask (output, extractMask)
    
output.save(outRaster)

try:
    arcpy.Delete_management(os.path.join(arcpy.env.scratchWorkspace, "reproj.shp"))
except:
    pass
    
##else:
##    output = rclss.save(outRaster)
    
# Set output raster as parameter
arcpy.SetParameterAsText(3, outRaster)

# Clean up temporary workspace
try:
    arcpy.Delete_management("in_memory")
    arcpy.Delete_management(arcpy.env.scratchFolder)
except:
    pass
