#------------------------------------------------------------------------------
# Copyright 2014 Esri
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#------------------------------------------------------------------------------

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# TestRadialLOS.py
# Created on: 2013-05-08 12:23:38.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: TestRadialLOS <observers> <input_surface> <output_rlos> <RADIUS2_to_infinity> <towerHeight>
# Description:
# Shows the areas visible and invisible to one or more observer points.
# ---------------------------------------------------------------------------
# history:
# 4/23/2015 - mf - Updates for Extent properties
# ==================================================

# Import arcpy module
import os, sys, traceback, math, decimal
import arcpy
from arcpy import env
from arcpy import sa

def zfactor(dataset):
    desc = arcpy.Describe(dataset)
    # if it's not geographic return 1.0
    if desc.spatialReference.type != "Geographic":
        return 1.0
    extent = desc.Extent
    extent_split = [extent.XMin,extent.YMin,extent.XMax,extent.YMax]

    top = float(extent_split[3])
    bottom = float(extent_split[1])

    #find the mid-latitude of the dataset
    if (top > bottom):
        height = (top - bottom)
        mid = (height/2) + bottom
    elif (top < bottom):  # Unlikely, but just in case
        height = bottom - top
        mid = (height/2) + top
    else: # top == bottom
        mid = top

    # convert degrees to radians
    mid = math.radians(mid)

    # Find length of degree at equator based on spheroid's semi-major axis
    spatial_reference = desc.SpatialReference
    semi_major_axis = spatial_reference.semiMajorAxis # in meters
    equatorial_length_of_degree = ((2.0 * math.pi * float(semi_major_axis))/360.0)

    # function:
    # Z-Factor = 1.0/(111320 * cos(mid-latitude in radians))
    decimal.getcontext().prec = 28
    decimal.getcontext().rounding = decimal.ROUND_UP
    a = decimal.Decimal("1.0")
    b = decimal.Decimal(str(equatorial_length_of_degree))
    c = decimal.Decimal(str(math.cos(mid)))
    zfactor = a/(b * c)
    zfactor = "%06f" % (zfactor.__abs__())
    return zfactor

def updateValue(fc, field, value):
    cursor = arcpy.UpdateCursor(fc)
    for row in cursor:
        row.setValue(field, value)
        cursor.updateRow(row)
    return


# Script arguments
input_surface = arcpy.GetParameterAsText(0) #Input Surface
RADIUS2_to_infinity = arcpy.GetParameterAsText(1) #Force visibility to infinity
towerClass = arcpy.GetParameterAsText(2) #Defensive Position Feature Class
descField= arcpy.GetParameterAsText(3) #Defensive Position Description Field
towerName = arcpy.GetParameterAsText(4) #Defensive Position Description
towerHeightField = arcpy.GetParameterAsText(5) #Defensive Position Height Field
# The name of the workspace in which the features should be stored
outWorkspace = arcpy.GetParameterAsText(6) #Output Workspace
# The name of the featureclass in which the features should be stored
outFeatureClassName = arcpy.GetParameterAsText(7) #Output Visibility

# Scrub the name
scrubbedFeatureClassName = arcpy.ValidateTableName(outFeatureClassName,outWorkspace)

# Put it all together
output_rlos = os.path.join(outWorkspace,scrubbedFeatureClassName)
if RADIUS2_to_infinity == '#' or not RADIUS2_to_infinity:
    RADIUS2_to_infinity = "false"

terrestrial_refractivity_coefficient = 0.13
polygon_simplify = "SIMPLIFY"


delete_me = []
DEBUG = True
arcpy.AddMessage("Using scratch GDB of: " + arcpy.env.scratchWorkspace)


if RADIUS2_to_infinity == 'true':
    RADIUS2_to_infinity = True
else:
    RADIUS2_to_infinity = False

if DEBUG == True:
    arcpy.AddMessage("Input surface is " + input_surface)
    arcpy.AddMessage("Output is " + output_rlos)
    arcpy.AddMessage("Tower Name is " + towerName)
    arcpy.AddMessage("Tower Height Field is " + towerHeightField)
    arcpy.AddMessage("RADIUS2_to_infinity:" + str(RADIUS2_to_infinity))
    arcpy.AddMessage("outWorkspace:" + str(outWorkspace))
    arcpy.AddMessage("outFeatureClassName:" + outFeatureClassName)
    arcpy.AddMessage("scrubbedFeatureClassName: " + scrubbedFeatureClassName)

###############################
# Set the snap raster to the input surface
###############################
if DEBUG == True:
        arcpy.AddMessage("Setting snap raster to: " + input_surface)
env.snapRaster = input_surface


# mf 11/26/2014: Fix issue #150
arcpy.SelectLayerByAttribute_management(towerClass, "CLEAR_SELECTION")

try:

    # Local variables:
    Output_Feature_Class = output_rlos
    Output_Feature_Class__2_ = Output_Feature_Class
    sr = arcpy.SpatialReference()
    sr.factoryCode = 4326
    sr.create()
    GCS_WGS_1984 = sr
    #GCS_WGS_1984 = arcpy.SpatialReference(r"WGS 1984")
    Feature_Class = ""

    # Process: sourceRLOSscript
     # get/set initial environment
    env.overwriteOutput = True
    installInfo = arcpy.GetInstallInfo("desktop")
    installDirectory = installInfo["InstallDir"]
    
    # Get the towerHeight value from the towerHeightField (assuming that there is only one)
    towerHeight = None
    towerRows = arcpy.da.SearchCursor(towerClass,[descField,towerHeightField])
    for towerRow in towerRows:
        if str(towerRow[0]) == str(towerName):
            towerHeight = float(towerRow[1])
    if DEBUG == True: arcpy.AddMessage("towerHeight: " + str(towerHeight))
    del towerRow
    del towerRows

    # get observer's vibility modifier maximums
    obsMaximums = {'SPOT':None,'OFFSETA':towerHeight, 'RADIUS2':4000, 'REMOVE_SPOT':False}
    #TODO: Why is RADIUS2 hardcoded to be 4000?
    
    removeSPOT = obsMaximums['REMOVE_SPOT']
    if (removeSPOT == True):
        if DEBUG == True:
            arcpy.AddMessage("Observer SPOT is <NULL>, deleteing field ...")
        arcpy.DeleteField_management(observers,"SPOT")

    # Remove special characters and replace spaces with underscores
    scrubbedTowerName = ''.join(e for e in towerName if (e.isalnum() or e == " " or e == "_"))
    scrubbedTowerName = scrubbedTowerName.replace(" ", "_")
    if DEBUG == True:
        arcpy.AddMessage("towerName: " + towerName)
        arcpy.AddMessage("scrubbedTowerName: " + scrubbedTowerName)
        arcpy.AddMessage("DescField: " + descField)
    tower = arcpy.MakeFeatureLayer_management(towerClass, scrubbedTowerName, descField + " = '" + str(towerName) +"'")

    # Do a Minimum Bounding Geometry (MBG) on the tower
    observers_mbg = os.path.join(env.scratchWorkspace,"observers_mbg_towerlos")
    delete_me.append(observers_mbg)
    arcpy.AddMessage("Finding observer's minimum bounding envelope using " + observers_mbg)
    arcpy.MinimumBoundingGeometry_management(tower,observers_mbg)

    # Now find the center of the (MBG)
    arcpy.AddMessage("Finding center of Defensive Position ...")
    mbgCenterPoint = os.path.join(env.scratchWorkspace,"mbgCenterPoint_towerlos")
    mbgExtent = arcpy.Describe(observers_mbg).extent
    mbgSR = arcpy.Describe(observers_mbg).spatialReference
    mbgCenterX = mbgExtent.XMin + (mbgExtent.XMax - mbgExtent.XMin)
    mbgCenterY = mbgExtent.YMin + (mbgExtent.YMax - mbgExtent.YMin)
    arcpy.CreateFeatureclass_management(os.path.dirname(mbgCenterPoint),os.path.basename(mbgCenterPoint),"POINT","#","DISABLED","DISABLED",mbgSR)
    mbgShapeFieldName = arcpy.Describe(mbgCenterPoint).ShapeFieldName
    rows = arcpy.InsertCursor(mbgCenterPoint)
    feat = rows.newRow()
    feat.setValue(mbgShapeFieldName,arcpy.Point(mbgCenterX,mbgCenterY))
    rows.insertRow(feat)
    del rows
    delete_me.append(mbgCenterPoint)

    # Get the maximum radius of the Defensive Position
    maxRad = obsMaximums['RADIUS2']
    maxOffset = obsMaximums['OFFSETA']
    z_factor = float(zfactor(tower))
    horizonDistance = 0.

    if RADIUS2_to_infinity == True:
        # if going to infinity what we really need is the distance to the horizon based on height/elevation
        arcpy.AddMessage("Finding horizon distance ...")
        result = arcpy.GetCellValue_management(input_surface, str(mbgCenterX) + " " + str(mbgCenterY))
        centroid_elev = result.getOutput(0)
        R2 = float(centroid_elev) + float(maxOffset)
        R = 6378137.0 # length, in meters, of semimajor axis of WGS_1984 spheroid.
        horizonDistance = math.sqrt(math.pow((R + R2),2) - math.pow(R,2))

        horizonExtent = str(mbgCenterX - horizonDistance) + " " + str(mbgCenterY - horizonDistance) + " " + str(mbgCenterX + horizonDistance) + " " + str(mbgCenterY + horizonDistance)
    else:
        pass

    # reset center of AZED using Lat/Lon of MBG center point
    # Project point to WGS 84
    arcpy.AddMessage("Recentering Azimuthal Equidistant to centroid ...")
    mbgCenterWGS84 = os.path.join(env.scratchWorkspace,"mbgCenterWGS84_towerlos")
    arcpy.Project_management(mbgCenterPoint,mbgCenterWGS84,GCS_WGS_1984)
    arcpy.AddXY_management(mbgCenterWGS84)
    pointx = 0.0
    pointy = 0.0
    shapeField = arcpy.Describe(mbgCenterWGS84).ShapeFieldName
    rows = arcpy.SearchCursor(mbgCenterWGS84)
    for row in rows:
        feat = row.getValue(shapeField)
        pnt = feat.getPart()
        pointx = pnt.X
        pointy = pnt.Y
    del row
    del rows
    # write new central meridian and latitude of origin...
    strAZED = 'PROJCS["World_Azimuthal_Equidistant",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION    ["Azimuthal_Equidistant"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",' + str(pointx) + '],PARAMETER["Latitude_Of_Origin",' + str(pointy) + '],UNIT["Meter",1.0],AUTHORITY    ["ESRI",54032]]'
    if DEBUG == True:
        arcpy.AddMessage("strAZED: " + strAZED)
        arcpy.AddMessage("pointx: " + str(pointx))
        arcpy.AddMessage("pointy: " + str(pointy))
    delete_me.append(mbgCenterWGS84)

    # Determine the proper buffer distance based on whether visibility should be generated to Infinity
    bufferDistance = obsMaximums['RADIUS2']
    if RADIUS2_to_infinity == True:
        arcpy.AddMessage("Running viewshed to infinity...")
        bufferDistance = horizonDistance
        if DEBUG == True:
            arcpy.AddMessage("Creating a buffer from the observer points to the horizion distance.  infinity (horizon): " + str(bufferDistance))
    else:
        arcpy.AddMessage("Running viewshed to the maximum observer radius...")
        if DEBUG == True:
            arcpy.AddMessage("Creating a buffer from the observer points to the horizion distance.  maximum observer radius: " + str(bufferDistance))


    # Create the processing buffer to the appropriate distance
    mbgBuffer = os.path.join(env.scratchWorkspace,"mbgBuffer_towerlos")
    arcpy.Buffer_analysis(observers_mbg,mbgBuffer,obsMaximums['RADIUS2'])
    delete_me.append(mbgBuffer)
    if DEBUG == True:
        arcpy.AddMessage("Projecting the buffer to Azimuthal Equidistant")
    mbgBufferPrj = os.path.join(env.scratchWorkspace,"mbgBufferPrj_towerlos")
    arcpy.Project_management(mbgBuffer,mbgBufferPrj,strAZED)
    delete_me.append(mbgBufferPrj)
    mbgBufferPrjExtent = arcpy.Describe(mbgBufferPrj).extent
    if DEBUG == True:
        arcpy.AddMessage("Setting procesing extent to: " + str(mbgBufferPrjExtent))
    env.extent = mbgBufferPrjExtent

    # Project surface to the new AZED
    extract_prj = os.path.join(env.scratchWorkspace,"input_surface_prj_towerlos")
    arcpy.AddMessage("Projecting surface ...")
    arcpy.ProjectRaster_management(input_surface,extract_prj,strAZED)
    delete_me.append(extract_prj)

    # Project observers to the new AZED
    obs_prj = os.path.join(env.scratchWorkspace,"obs_prj_towerlos")
    arcpy.AddMessage("Projecting observers ...")
    arcpy.Project_management(mbgCenterPoint,obs_prj,strAZED)
    
    #Add viewshed-utilized fields
    if DEBUG == True: arcpy.AddMessage("Adding OFFSETA field to: " + str(obs_prj))
    arcpy.AddField_management(obs_prj, "OFFSETA", "DOUBLE", "", "", "", "Observer Offset", "NULLABLE", "NON_REQUIRED", "")
    arcpy.CalculateField_management(obs_prj, "OFFSETA", maxOffset, "PYTHON_9.3", "")
    
    if DEBUG == True: arcpy.AddMessage("Adding RADIUS2 field to: " + str(obs_prj))
    arcpy.AddField_management(obs_prj, "RADIUS2", "DOUBLE", "", "", "", "Farthest distance", "NULLABLE", "NON_REQUIRED", "")
    arcpy.CalculateField_management(obs_prj, "RADIUS2", maxRad, "PYTHON_9.3", "")
    delete_me.append(obs_prj)

    # Project the MBG buffer to AZED
    obs_buf = os.path.join(env.scratchWorkspace,"obs_buf_towerlos")
    arcpy.Project_management(mbgBufferPrj,obs_buf,strAZED)
    delete_me.append(obs_buf)

    # determine the proper z factor
    arcpy.AddMessage("Calculating Z Factor ...")
    z_factor = float(zfactor(extract_prj))
    if DEBUG == True:
        arcpy.AddMessage("z_factor: " + str(z_factor))

    # Finally ... run Viewshed
    arcpy.AddMessage("Calculating Viewshed ...")
    if DEBUG == True:
        arcpy.AddMessage("terrestrial_refractivity_coefficient" + str(terrestrial_refractivity_coefficient))
    vshed = os.path.join(env.scratchWorkspace,"vshed_towerlos")
    delete_me.append(vshed)
    #Use the visibility tool instead of the viewshed. This will allow is to set the observation height offset for the tower
    #outVshed = sa.Visibility(extract_prj,mbgCenterPoint,z_factor=z_factor,curvature_correction="CURVED_EARTH",refractivity_coefficient=terrestrial_refractivity_coefficient,observer_offset=towerHeight)
    outVshed = sa.Viewshed(extract_prj,obs_prj,z_factor,"CURVED_EARTH",terrestrial_refractivity_coefficient)
    outVshed.save(vshed)

    # Raster To Polygon
    arcpy.AddMessage("Converting to polygons ...")
    ras_poly = os.path.join(env.scratchWorkspace,"ras_poly_towerlos")
    arcpy.RasterToPolygon_conversion(vshed,ras_poly,polygon_simplify)
    delete_me.append(ras_poly)

    # clip output polys to buffer
    if RADIUS2_to_infinity != True:
        if DEBUG == True:
            arcpy.AddMessage("Clipping output polygon ...")
        out_buf = os.path.join(env.scratchWorkspace,"out_buf_towerlos")
        arcpy.Buffer_analysis(obs_prj,out_buf,"RADIUS2")
        delete_me.append(out_buf)
        arcpy.Clip_analysis(ras_poly,out_buf,output_rlos)
    else:
        if DEBUG == True:
            arcpy.AddMessage("Copying output polygon ...")
        arcpy.CopyFeatures_management(ras_poly, output_rlos)

    arcpy.AddMessage("Adding visibility fields ...")
    arcpy.AddField_management(output_rlos, "visibility", "DOUBLE", "", "", "", "Observer Visibility", "NULLABLE", "NON_REQUIRED", "")
    arcpy.CalculateField_management(output_rlos, "visibility", "!gridcode!", "PYTHON", "")

    # Add the layer to the map
    #UPDATE
    #layerSymLocation = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', 'layers'))
    #mxd = arcpy.mapping.MapDocument('CURRENT')
    #df = arcpy.mapping.ListDataFrames(mxd)[0]
    #layerToAdd = arcpy.mapping.Layer(output_rlos)
    #arcpy.ApplySymbologyFromLayer_management(layerToAdd, layerSymLocation + "\Radial Line Of Sight Output.lyr")
    #arcpy.mapping.AddLayer(df, layerToAdd, "AUTO_ARRANGE")
        
    #Set the tool output
    arcpy.SetParameter(8,output_rlos)
    

except arcpy.ExecuteError:
    # Get the tool error messages
    msgs = arcpy.GetMessages()
    arcpy.AddError(msgs)
    #print msgs #UPDATE
    print(msgs)

except:
    # Get the traceback object
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]

    # Concatenate information together concerning the error into a message string
    pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
    msgs = "\nArcPy ERRORS:\n" + arcpy.GetMessages() + "\n"

    # Return python error messages for use in script tool or Python Window
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)

    # Print Python error messages for use in Python / Python Window
    #print pymsg + "\n" #UPDATE
    print(pymsg + "\n")
    #print msgs #UPDATE
    print(msgs)

finally:
    # cleanup intermediate datasets
    if DEBUG == True: arcpy.AddMessage("Removing intermediate datasets...")
    for i in delete_me:
        if DEBUG == True: arcpy.AddMessage("Removing: " + str(i))
        #arcpy.Delete_management(i)
    if DEBUG == True: arcpy.AddMessage("Done")

