ERROR 000989 Issue with Python Script and Field Calculating distance between two long lat points

Question

very new to scripting and I have been struggling with my script for a while. I'm not sure why it doesn't work in field calculator. I am trying to calculate a list of distance between two long lat points. The following is the code that I have implemented in the Pre-Logic Script Code for FieldName1. I am currently using ArcGIS 10.1

def Distance(Latitude, Latitude2, Longitude, Longitude2)
 Longitude, Latitude, Longitude2, Latitude2 = map(math.radians, [Longitude, Latitude, Longitude2, Latitude2])
 dlong = Longitude2-Longitude
 dlat = Latitude2 - Latitude
 a = math.sin(dlat/2)**2 + math.cos(Latitude * math.cos(Latitude2) * math.sin(dlong/2)**2
 c = 2 * math.asin(math.sqrt(a))
 Distance = 6371 * C * 1000

The following are sample points for reference: Latitude :48.5761 Longitude : -124.31111 Latitude 2: 48.57611 Longitude 2: -124.31111

I kept getting error 000989: Python syntax error: Parsing error SyntaxError:invalid syntax(line1). It also adds that the parameter are not valid.

Answer 1

According to this tutorial the map function returns a list:

>>> items = [1, 2, 3, 4, 5]
>>> 
>>> def sqr(x): return x ** 2
>>> list(map(sqr, items))
[1, 4, 9, 16, 25]
>>> 

So your code should be something like this:

def Distance(Latitude, Latitude2, Longitude, Longitude2)
 pList = map(math.radians, [Longitude, Latitude, Longitude2, Latitude2])
 Longitude = pList[0]
 Latitude = pList[1]
 Longitude2 = pList[2]
 Latitude2 = pList[3]
 dlong = Longitude2-Longitude
 dlat = Latitude2 - Latitude
 a = math.sin(dlat/2)**2 + math.cos(Latitude * math.cos(Latitude2) * math.sin(dlong/2)**2
 c = 2 * math.asin(math.sqrt(a))
 Distance = 6371 * C * 1000

Which would be the same as writing:

def Distance(Latitude, Latitude2, Longitude, Longitude2)
 Longitude = math.radians(Longitude)
 Latitude = math.radians(Latitude)
 Longitude2 = math.radians(Longitude2)
 Latitude2 = math.radians(Latitude2)
 dlong = Longitude2-Longitude
 dlat = Latitude2 - Latitude
 a = math.sin(dlat/2)**2 + math.cos(Latitude * math.cos(Latitude2) * math.sin(dlong/2)**2
 c = 2 * math.asin(math.sqrt(a))
 Distance = 6371 * C * 1000

Answer 2

There are 4 problems with this function, one causing the SyntaxError:invalid syntax(line1) message.

1. The function definition needs to end with a colon :
2. The call to math.cos(Latitude needs a close parenthesis )
3. The variable name c needs to be lower case or upper case every time you use it.
4. After calculating the result, the function needs to return it.

def Distance(Latitude, Latitude2, Longitude, Longitude2):
 Longitude, Latitude, Longitude2, Latitude2 = map(math.radians, [Longitude, Latitude, Longitude2, Latitude2])
 dlong = Longitude2-Longitude
 dlat = Latitude2 - Latitude
 a = math.sin(dlat/2)**2 + math.cos(Latitude) * math.cos(Latitude2) * math.sin(dlong/2)**2
 c = 2 * math.asin(math.sqrt(a))
 Distance = 6371 * c * 1000
 return Distance

While this now works, there are 2 other things that would be good practice when you write longer Python programs:

1. Don't re-use the names of your arguments (Latitude, Latitude2, Longitude, Longitude2) with newly assigned values. Instead give new names that describe how you have modified them, such as (LatitudeRadians, Latitude2Radians, LongitudeRadians, Longitude2Radians)
2. Don't create a variable Distance with the same name as your function, even if your variable name has to be something arbitrary like DistanceMeters.

While this answers your question about getting this script to run, there is probably a still better way to get the distance, but I will put that in a separate answer.

Answer 3

Calculating Distance with ArcPy PointGeometry

I posted a separate answer on how to tidy your existing syntax. But geometry on the Earth's surface is hard, and ArcMap is there to help us. It provides objects to store geometry, with methods to compare them. To use them to solve your problem, you need to:

1. Construct a SpatialReference object. This must have the datum your coordinates are relative to.
2. Construct Point objects. This just stores XY (and possibly Z) values and has no notion of coordinate reference systems.
3. Construct PointGeometry objects using your Point and SpatialReference objects.
4. Call the angleAndDistanceTo method to get the distance.

This will handle matters like the Earth being a bit squashed at the poles, instead of treating it as a fixed 6,371km radius. What difference does it make? From (0°N,0°E) to (0°N,90°E), angleAndDistanceTo returns a distance 11 kilometers longer. From (0°N,0°E) to (90°N,0°E), angleAndDistanceTo returns a distance 5 kilometers shorter.

Easy Traps

1. If you don't construct a SpatialReference object appropriate to your coordinates, your distances will be meaningless.
2. The constructor for arcpy.Point takes the X value (latitude), then the Y (longitude). To make the arguments clear, the sample code below uses keyword arguments.
3. The angleAndDistanceTo method has an argument that can be "GEODESIC", "GREAT_ELLIPTIC", "LOXODROME", "PLANAR" or "PRESERVE_SHAPE". Make sure you know what these mean and choose the one appropriate to your task. If using "PLANAR" you will need to reproject your PointGeometry.
4. There is a method called distanceTo, except it always returns distances in the coordinate system of your PointGeometry. If you are in WGS84, your "distances" will be in degrees. You cannot specify a calculation method the way that angleAndDistanceTo does.

Place some points on your map and use the ArcMap Measure Tool to measure the distance between them. You should find that, as long as you specify the same method, the distances it gives you are the same as from your Field Calculator function.

The Code

def Distance(Latitude1, Latitude2, Longitude1, Longitude2):
 # Need to define a coordinate reference system appropriate to our inputs.
 # Here I assume they are coming from GPS, so WGS 1984.
 # These can also be constructed using their EPSG numbers (4326 for WGS 1984)
 spatialReferenceWGS84 = arcpy.SpatialReference("WGS 1984")
 # Make sure to specify which argument to the Point constructor is X (Longitude) and Y (Latitude).
 point1 = arcpy.Point(X = Longitude1, Y = Latitude1)
 pointGeometry1 = arcpy.PointGeometry(point1, spatialReferenceWGS84)
 point2 = arcpy.Point(X = Longitude2, Y = Latitude2)
 pointGeometry2 = arcpy.PointGeometry(point2, spatialReferenceWGS84)
 # Do not use the pointGeometry.distanceTo function. It will return answers in degrees.
 # For the different method values, see:
 # https://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-classes/pointgeometry.htm#M2_GUID-5853F68B-6164-4FDA-81C9-C0DEF50F0830
 # They are "GEODESIC", "GREAT_ELLIPTIC", "LOXODROME", "PLANAR" and "PRESERVE_SHAPE"
 angleBetweenPoints,distanceBetweenPoints = pointGeometry1.angleAndDistanceTo(pointGeometry2, "GEODESIC")
 # Distance between the two points is known, return it.
 return distanceBetweensPoints

Thank you to this answer that got me on the right track here.