Geometry: Points in the field calculator
Date: 2015-10-08 Modified: 2017-04-19 new ***
Included:
- Add X Y coordinates to table by distance or percent along poly* features ***
- Distance to a specific point
- Cumulative distance
- Inter-point distance
- Azimuth to a point
- Angle between consecutive points
- Convert Azumith to compass bearing ***
- Convert Degrees decimal minutes to decimal degrees ***
Purpose:
Many posts on GeoNet come from people looking to write a script or produce a tool which can be handled in a much simpler fashion. This is first in a series of posts that will address how to use the field calculator to its fullest. The use, and perhaps the limitations to using, the field calculator should be guided by the following considerations:
- you need to perform some rudimentary task for which there is no existing tool (ie. it is so simple to do...what is the point of having it builtin)
- your datasets are relatively small ...
- of course this is a subjective classification... so if processing takes more than a minute... you have left the realm of small
- you don't have to perform the task that often...
- basically you are trying to fix up an oversight is the initial data requirements, or the requirements have changed since project inception
- you got the data elsewhere and you are trying to get it to meet project standards
- you want to explore
- you have no clue how to script but want to dabble as a prelude to full script development.
- whatever else I left out ....
Notes:
- Make sure you have projected data and a double field or nothing will make sense. I care not about performing great circle calculations or determining geodetic areas. The field calculator...IMHO... is not the place to do that.
- Ensure that the sequence of points is indeed correct. If you have been editing or fiddling around with the data, make sure that nothing is amiss with the data. What you see, is often not what is. You can get erroneous results from any environment
- Sorting will not change the outcome, the results will be in feature order.
- If you want them in another order, then use tools to do that.
- I will be paralleling geometric calculations using NumPy and Arcpy is a separate series which removes this issue.
- Selections can be used but they will be in feature order.
- sorting and selecting can get it to look like you want it ... but it doesn't make it so... want and is, are two separate concepts like project and define projection (but I digress)
- VB anything is not covered.
- It is time to move on to the next realm in the sequence and evolution of languages.
- in ArcGIS Pro, VB isn't even an option, so get used to it
- code is written verbosely, for the most part
- I will try to use the simplest of programming concepts... simplicity trumps coolness and speed.
- parallel, optional and/or optimized solutions will be discussed elsewhere.
- the math module is already imported
Useage:
For all functions, do the following:
- select the Python parser ... unfortunately, it can't be selected by default
- toggle on, Show code block
- Pre-logic Script Code:
- paste the code in this block
- Expression box
- paste the function call in the expression box
- ensure that there is preceding space before the expression...it will generate an error otherwise
- ie dist_between(!Shape!) ... call the dist_between function ... the shape field is required
- field names are surrounded by exclamation marks ( ! ) ... the shapefile and file geodatabase standard
----------------------------------------------------------------------------------------------------------------------------------------
Add X or Y coordinate to table field by distance or percentage
For use in tables, to retrieve the X or Y coordinates in an appropriate field. See the header for requirements.
Pre-logic Script Code:
def pnt_along(shape, value=0.0, use_fraction=False, XorY="X"):
"""Position X or Y coordinate, x/y meters or decimal fraction along a line.
:Requires:
:--------
: shape field: python parser use !Shape!
: value: (distance or decimal fraction, 0-1)
: use_fraction: (True/False)
: XorY: specify X or Y coordinates
:
:Returns: the specified coordinate (X or Y) meters or % along line or boundary
:-------
:
:Useage: pnt_along(!Shape!, 100, False, "X") # X, 100 m from start point
:
"""
XorY = XorY.upper()
if use_fraction and (value > 1.0):
value = value/100.0
if shape.type.lower() == "polygon":
shape = shape.boundary()
pnt = shape.positionAlongLine(value,use_fraction)
if XorY == 'X':
return pnt.centroid.X
else:
return pnt.centroid.Yexpression =
pnt_along(!Shape!, 100, False, "X") # X, 100 m from start point--------------------------------------------------------------------------------------------------------------
Distance to a specific point
Calculate the distance to a specific point within a dataset. For example, you can determine the distance from a point cloud's center to other points. Features not in the file can be obtained by other means. Useful to get summary information on inter-point distances as a prelude to a full clustering or KD-means study.
For example, the potential workflow to get the distance of every point to the point clouds center might include the following steps:
- add X and Y style fields to contain the coordinates
- use Field Statistics to get the mean center (use a pen to record them... not everything is automagic)
- add a Dist_to field to house the calculations
- use the script below
Pre-logic Script Code:
""" dist_to(shape, from_x, from_y)
input: shape field, origin x,y
returns: distance to the specified point
expression: dist_to(!Shape!, x, y)
"""
def dist_to(shape, from_x, from_y):
x = shape.centroid.X
y = shape.centroid.Y
distance = math.sqrt((x - from_x)**2 + (y - from_y)**2)
return distanceexpression =
dist_to(!Shape!, x, y)-----------------------------------------------------------------------------------------------------------------------------------------
Cumulative Distance
To determine the distance from a point in a data table to every other point in sequence. Essentially a sequential perimeter calculation.
Pre-logic Script Code:
""" input shape field: returns cumulative distance between points
dist_cumu(!Shape!) #enter into the expression box"""
x0 = 0.0
y0 = 0.0
distance = 0.0
def dist_cumu(shape):
global x0
global y0
global distance
x = shape.centroid.X
y = shape.centroid.Y
if x0 == 0.0 and y0 == 0.0:
x0 = x
y0 = y
distance += math.sqrt((x - x0)**2 + (y - y0)**2)
x0 = x
y0 = y
return distanceexpression =
dist_cum(!Shape!)----------------------------------------------------------------------------------------------------------------------------------------
Inter-point distance
Determine distances between sequential point pairs (not cumulative).
Pre-logic Script Code:
""" dist_between(shape)
input: shape field
returns: distance between successive points
expression: dist_between(!Shape!)
"""
x0 = 0.0
y0 = 0.0
def dist_between(shape):
global x0
global y0
x = shape.centroid.X
y = shape.centroid.Y
if x0 == 0.0 and y0 == 0.0:
x0 = x
y0 = y
distance = math.sqrt((x - x0)**2 + (y - y0)**2)
x0 = x
y0 = y
return distanceexpression =
dist_between(!Shape!)-----------------------------------------------------------------------------------------------------------------------------------------
Azimuth to a specific point
Determine the azimuth to a specific point, for example, a point cloud's center.
Pre-logic Script Code:
""" azimuth_to(shape, from_x, from_y)
input: shape field, from_x, from_y
returns: angle between 0 and <360 between a specified point and others
expression: azimuth_to(!Shape!, from_x, from_y)
"""
def azimuth_to(shape, from_x, from_y):
radian = math.atan((shape.centroid.X - from_x)/(shape.centroid.Y - from_y))
degrees = math.degrees(radian)
if degrees < 0:
return degrees + 360.0
else:
return degrees
expression =
azimuth_to(!Shape!,from_x, from_y)
-----------------------------------------------------------------------------------------------------------------------------------------
Angle between successive points
Determine the angle between points, for example, angle changes between waypoints.
Pre-logic Script Code:
""" angle_between(shape)
input: shape field
returns: angle between successive points,
NE +ve 0 to 90, NW +ve 90 to 180,
SE -ve <0 to -90, SW -ve <-90 to -180
expression: angle_between(!Shape!)
"""
x0 = 0.0
y0 = 0.0
angle = 0.0
def angle_between(shape):
global x0
global y0
x = shape.centroid.X
y = shape.centroid.Y
if x0 == 0.0 and y0 == 0.0:
x0 = x
y0 = y
return 0.0
radian = math.atan2((shape.centroid.Y - y0),(shape.centroid.X - x0))
angle = math.degrees(radian)
x0 = x
y0 = y
return angle
expression =
angle_between(!Shape!)
-----------------------------------------------------------------------------------------------------------------------------------------
Line direction or Azimuth to Compass Bearing
Can be used to determine the direction/orientation between two points which may or may not be on a polyline. Provide the origin and destination points. The origin may be the 0,0 origin or the beginning of a polyline or a polyline segment.
def line_dir(orig, dest, fromNorth=False):
"""Direction of a line given 2 points
: orig, dest - two points representing the start and end of a line.
: fromNorth - True or False gives angle relative to x-axis)
:
"""
orig = np.asarray(orig)
dest = np.asarray(dest)
dx, dy = dest - orig
ang = np.degrees(np.arctan2(dy, dx))
if fromNorth:
ang = np.mod((450.0 - ang), 360.)
return ang-----------------------------------------------------------------------------------------------------------------------------------------
Convert Azimuth to Compass Bearing
Once the Azimuth to a particular point is determined, this can be converted to a compass direction, centered in 22.5 degree bands. The type of compass can be altered to suit... see script header.
Pre-logic Script Code:
import numpy as np
global a
global c
c = np.array(['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE',
'S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW', 'N'])
a = np.arange(11.25, 360., 22.5)
def compass(angle):
"""Return the compass direction based on supplied angle.
:Requires:
:--------
: angle - angle(s) in degrees, no check made for other formats.
: - a single value, list or np.ndarray can be used as input.
: - angles are assumed to be from compass north, alter to suit.
:
:Returns: The compass direction.
:-------
:
:Notes:
:-----
: Compass ranges can be altered to suit the desired format.
: See various wiki's for other options. This incarnation uses 22.5
: degree ranges with the compass centered on the range.
: ie. N between 348.75 and 11.25 degrees, range equals 22.5)
:
:----------------------------------------------------------------------
"""
if isinstance(angle, (float, int, list, np.ndarray)):
angle = np.atleast_1d(angle)
comp_dir = c[np.digitize(angle, a)]
if len(comp_dir) == 1:
comp_dir[0]
return comp_direxpression =
compass(!Azimuth_Field_Name!) # python parser, field name enclosed in quotes:--------------------------------------------------------------------------------------------------------------------------------------------------------------
Degrees decimal minutes to decimal degrees
def ddm_ddd(a, sep=" "):
""" convert decimal minute string to decimal degrees
: a - degree, decimal minute string
: sep - usually a space, but check
"""
d, m = [abs(float(i)) for i in a.split(sep)]
sign = [-1, 1][d < 0]
dd = sign*(d + m/60.)
return ddFor example ddm_ddd(!YourStringField!, sep=" ") will convert a string/text field into a double in your new field