# -*- coding: utf-8 -*-
""" ---------------------------------------------------------------------------
Name:       SSM_OnLand.py
Version:    2.10    12/18/2018
Author:     Thomas Follett
            Whale Telemetry Group
            Marine Mammal Institute
            Oregon State University

Description: Estimated locations of tracked marine species (e.g. from telemetry devices)
	        should be confined to being over water. However, for various reasons, locations
	        can sometimes be estimated over land.
	        This Python script uses the uncertainty associated with each estimated location
	        to correct locations erroneously estimated over land, in this case from the 
	        output of a state-space model (SSM) applied to Argos satellite locations. 
	        The 95% credible limits in longitude and latitude provided by the SSM are used 
	        as the semi-major and semi-minor axes of a polyline ellipse around a <<single>> 
	        selected location. The polyline is then converted to a polygon and then the 
	        portions over land are erased. The centroid of the remaining ocean-based portion 
	        of the polygon is then written to a new point feature class with the attributes
	        of the original input point preserved.

Requirements: ArcMap v10.3+

Usage:  Rename this text file to "SSM_OnLand.py" (remove the .txt extension)
        Create a new script tool in ArcMap. Go to the tool Properties,
        click the Source tab and browse to this file's location.
        Parameters:
        Display Name   Data Type               Type       Direction
        -----------------------------------------------------------
        Input Points | Point Feature Layer   | Required | Input
        Land Polygon | Polygon Feature Layer | Required | Input
        Output Layer | Feature Dataset       | Derived  | Output

NOTE:   Tool users must select a single point before running the tool.

Reference: Jiménez-López, M.E., D.M. Palacios, A. Jaramillo, J. Urbán,
           and B.R. Mate. Accepted. Fin whale movements in the Gulf of
           California, Mexico, from satellite telemetry. PLoS ONE 2019
---------------------------------------------------------------------------"""

import arcpy
from arcpy import env
import collections

in_point = arcpy.GetParameterAsText(0)
in_land =  arcpy.GetParameterAsText(1)

sr = arcpy.SpatialReference(4326)
env.overwriteOutput = True

# trap for single point selection
if int(arcpy.GetCount_management(in_point).getOutput(0)) > 1:
    arcpy.AddMessage("You must select ONE point before tool use.")
    quit()
# add major/minor fields to input, if they don't already exist
if arcpy.ListFields(in_point, "major"):
    pass
else:
    arcpy.AddField_management(in_point,"major","DOUBLE")
if arcpy.ListFields(in_point, "minor"):
    pass
else:
    arcpy.AddField_management(in_point,"minor","DOUBLE")
# cursor field list
fields = ["FID",   #0
        "long025", #1
        "lat025",  #2
        "long",    #3
        "lat",     #4
        "long975", #5
        "lat975",  #6
        "major",   #7
        "minor"]   #8
# update the selected point with major & minor axes lengths
with arcpy.da.UpdateCursor(in_point, fields) as uCur:
    row = uCur.next()
    # X axis of error distribution = semimajor
    row[7] = arcpy.Polyline(arcpy.Array([arcpy.Point(row[3],row[2]),  #long, lat025
                                         arcpy.Point(row[3],row[6])]),#long, lat975
                                         sr).getLength("GEODESIC", "METERS")
    # Y axis of error distribution = semiminor
    row[8] = arcpy.Polyline(arcpy.Array([arcpy.Point(row[1],row[4]), #long025, lat
                                         arcpy.Point(row[5],row[4])]),#long975, lat
                                         sr).getLength("GEODESIC", "METERS")

# optionally you could choose to switch the axes ----------------------------------------
##    # Add a new Boolean parameter to the tool at the third position
##    # Move to line 48: invert =  arcpy.GetParameterAsText(2)
##    if invert:
##        # X axis of error distribution = semiminor
##        row[8] = arcpy.Polyline(arcpy.Array([arcpy.Point(row[3],row[2]),  #long, lat025
##                                             arcpy.Point(row[3],row[6])]),#long, lat975
##                                             sr).getLength("GEODESIC", "METERS")
##        # Y axis of error distribution = semimajor
##        row[7] = arcpy.Polyline(arcpy.Array([arcpy.Point(row[1],row[4]), #long025, lat
##                                           arcpy.Point(row[5],row[4])]),#long975, lat
##                                             sr).getLength("GEODESIC", "METERS")
##      # modify line 119: arcpy.SetParameterAsText(3, candidate)
#-----------------------------------------------------------------------------------------
    uCur.updateRow(row)

# create ellipse from major& minor axes
ellipse_name = arcpy.CreateUniqueName("ellipse_line")
ell_line = arcpy.TableToEllipse_management (in_point, ellipse_name,
            "long", "lat", "major", "minor", "METERS", "", "", "FID", sr)
# create a polygon from polyline ellipse
poly_name = arcpy.CreateUniqueName("in_memory\\ellipse_poly")
polygon = arcpy.FeatureToPolygon_management(ell_line, poly_name, None, "ATTRIBUTES", in_point)
# erase land overlay from polygon ellipse
erase_name = arcpy.CreateUniqueName("ellipse_erase")
erased = arcpy.Erase_analysis(polygon, in_land, erase_name)
# write centroid of remainder to a new point featureclass
pt_name = arcpy.CreateUniqueName("centroid")
candidate = arcpy.FeatureToPoint_management(erased, pt_name)
# append lat & lon to point
arcpy.AddXY_management(candidate)

arcpy.SetParameterAsText(2, candidate)