here the thing-
in the same layer, I have POINT_X, POINT_Y and POINT_X_1, POINT_Y_1, as decimal degrees.
I´m trying to find the linear distance in meters from the 2 points.
I´m trying to get it with pyhton, using the below code, but I get always the same syntax error
Any help appreciated
Thanks
Joao
pre-logic code
ef Haversine(Lat0, Lat1, Lon0, Lon1):
[INDENT]dlat = Lat0 - Lat1
dlon = Lon0 - Lon1
a = (math.sin(math.radians(dlat)/2 )*math.sin(math.radians(dlat)/2 )) + math.cos(math.radians(Lat1) )*math.cos(math.radians( Lat0 ) )*(math.sin(math.radians(dlon)/2 )*math.sin(math.radians(dlon)/2 ))
c = 2*math.atan2( math.sqrt(a ),math.sqrt(1-a ))
H = c * 6371
return H
[/INDENT]
Code:
Haversine( !POINT_Y! , !POINT_Y_1! , !POINT_X! , !POINT_X_1! )
Solved! Go to Solution.
At the moment this process is comparing apples and oranges and it will give incorrect results.
This said, the process needed is to build in a reproject function to take the coordinates from Geographic Coordinates (decimal degrees) into Map units (meters).
#Todo Get distance between 2 point
Step 1: Build up representation of 1st point in memory as PointGeometry Object
Step 2; Reproject point in map units
Step 3: Build up representation of 2nd point in memory as PointGeometryObject
Step 4: Reproject this into Map Units
Step 5: Compare the 2 points using Pythagoras Formula.
Attached is a function to perform the Reproject using the arcpy module
#projectSR
def projectSR (pt, sr): #sr = spatial reference, pt is geometry in memory
#constantts
gcs = arcpy.SpatialReference(2193) #2193 = factory code for NZ Transverse Mercator
gt = 'NZGD_2000_To_WGS_1984_1' #geographic transformation, if application
coords = []
ptgeom = arcpy.PointGeometry(pt, nztm).projectAs(gcs, gt)
coords.insert(0,(ptgeom.centroid.X))
coords.insert(1,(ptgeom.centroid.Y))
del ptgeom
return coords
Call projectSR to get the coordinates in map units meters and then you'll be able to use the Pythagoras theorem.
Susan
At the moment this process is comparing apples and oranges and it will give incorrect results.
This said, the process needed is to build in a reproject function to take the coordinates from Geographic Coordinates (decimal degrees) into Map units (meters).
#Todo Get distance between 2 point
Step 1: Build up representation of 1st point in memory as PointGeometry Object
Step 2; Reproject point in map units
Step 3: Build up representation of 2nd point in memory as PointGeometryObject
Step 4: Reproject this into Map Units
Step 5: Compare the 2 points using Pythagoras Formula.
Attached is a function to perform the Reproject using the arcpy module
#projectSR
def projectSR (pt, sr): #sr = spatial reference, pt is geometry in memory
#constantts
gcs = arcpy.SpatialReference(2193) #2193 = factory code for NZ Transverse Mercator
gt = 'NZGD_2000_To_WGS_1984_1' #geographic transformation, if application
coords = []
ptgeom = arcpy.PointGeometry(pt, nztm).projectAs(gcs, gt)
coords.insert(0,(ptgeom.centroid.X))
coords.insert(1,(ptgeom.centroid.Y))
del ptgeom
return coords
Call projectSR to get the coordinates in map units meters and then you'll be able to use the Pythagoras theorem.
Susan
If you still want to use Haversine, there are many examples on the web in python, for example