Hi all, I am creating a survey/sampling design using ArcGIS Pro 3.0. In order to randomize study plot locations, I am trying to create lines with angles that are randomized. Each line will go through a centerpoint in a grid and the plots will be centered on each end of the line. Any help with how to create the lines or angle them randomly would be much appreciated.
Solved! Go to Solution.
See another version below. This one will save your points too. It also creates your sampling boxes. There are four different methods you can use (see the "box_rotation" variable in the INPUTS section). Remember to update all the "INPUTS" based on your needs. Should be able to paste the code into the Python window again.
Only one version of the transect lines are shown (lines for the 'ANGLE_DIAMOND'-style boxes), but you can see all four styles of the sampling boxes. Some different sizes and distances from the origin points are also shown.
#-----------------------------------------------------------------------------#
# IMPORTS
import math
import os
import random
import sys
from typing import Literal, Union
import arcpy
#-----------------------------------------------------------------------------#
# INPUTS
# Full path to the geodatabase where the sampling points are, and where outputs will go.
workspace = r"C:TEMP\Random_Sample_Lines.gdb"
# Name of the point feature class in the above geodatabase that houses the
# sampling points.
in_pnt_fc = "_SAMPLE_POINT"
# Distance between sampling points (along x and y, not diagonals).
# Units match coordinate system.
point_spacing = 400 # feet, in my case
box_width = 50
box_rotation = 'ANGLE_BOX' # Literal[NORTH_BOX, NORTH_DIAMOND, ANGLE_BOX, ANGLE_DIAMOND]
# Name of the output feature class to be created.
fc_line = "TRANSECT_LINE"
fc_endpoints = "TRANSECT_LINE_ENDPOINTS"
# Leave this alone if you are not familiar with the f-string notation.
fc_box = f"TRANSECT_{box_rotation}"
# -----------------------------------------------------------------------------#
# SETTINGS
# This just overwrites previous outputs to help with testing.
arcpy.env.overwriteOutput = False
#-----------------------------------------------------------------------------#
# FUNCTIONS
def rotate_point(coordinate: tuple[int, int],
pivot_point: tuple[int, int],
angle: Union[int, float],
angle_format: Literal['DEGREES', 'RADIANS']='DEGREES'
) -> tuple[int, int]:
"""Rotates coordinate values around another point.
:coordinate: (x, y) coordinate pair to be rotated
:pivot_point: (x, y) coordinate pair around which point will be rotated
:angle: angle to rotate the coordinate
:angle_format: angle format, RADIANS or DEGREES
|(x_prime, y_prime)| rotated point coordinates"""
if angle_format == 'DEGREES':
angle = math.radians(angle)
x, y = coordinate
X, Y = pivot_point
x_trans = x - X
y_trans = y - Y
sin = math.sin(angle)
cos = math.cos(angle)
# Counter-clockwise rotation:
x_transprime = cos * x_trans - sin * y_trans
y_transprime = sin * x_trans + cos * y_trans
x_prime = x_transprime + X
y_prime = y_transprime + Y
# print(f"XY TRANS: {x_trans:,.0f}, {y_trans:,.0f}\n"
# f"XY TRANSPRIME: {x_transprime:,.0f}, {y_transprime:,.0f}\n"
# f"XY PRIME: {x_prime:,.0f}, {y_prime:,.0f}\n")
return (x_prime, y_prime)
def create_box(method: Literal['NORTH_BOX', 'NORTH_DIAMOND', 'ANGLE_BOX', 'ANGLE_DIAMOND'],
center_pnt: arcpy.Point,
width: Union[int, float],
angle: Union[int, float],
spat_ref: int) -> arcpy.Polygon:
"""Create a square box around a given point. Optionally rotate it.
:center_pnt: centroid of the square box
:width: width of side of box
:angle: angle the box should be rotated
|box_geom| geometry object representing the box"""
angle_dict = {'NORTH_BOX': 45,
'NORTH_DIAMOND': 0,
'ANGLE_BOX': angle + 45,
'ANGLE_DIAMOND': angle}
# Distance to move point based on known hypotenuse (1/2 width of box).
offset = math.sqrt(2 * ((width / 2) ** 2))
x, y = center_pnt.X, center_pnt.Y
Nx, Ny = x, y + offset
Ex, Ey = x + offset, y
Sx, Sy = x, y - offset
Wx, Wy = x - offset, y
box_corners = []
for coordinate in [(Nx, Ny), (Ex, Ey), (Sx, Sy), (Wx, Wy), (Nx, Ny)]:
box_corners.append(rotate_point(coordinate=coordinate,
pivot_point=(x, y),
angle=angle_dict[method],
angle_format='DEGREES'))
pnt_array = arcpy.Array([arcpy.Point(*corner) for corner in box_corners])
box_geom = arcpy.Polygon(pnt_array, spatial_reference=spat_ref)
return box_geom
#-----------------------------------------------------------------------------#
# MAIN
in_points = os.path.join(workspace, in_pnt_fc)
transect_lines = []
line_enpoint_dict = {}
sample_box_dict = {}
pnt_desc = arcpy.da.Describe(in_points)
pnt_oid_fld = pnt_desc['OIDFieldName']
pnt_spat_ref = pnt_desc['spatialReference'].factoryCode
with arcpy.da.SearchCursor(in_points, [pnt_oid_fld, 'SHAPE@XY']) as scurs:
for oid, (x, y) in scurs:
# Random angle in degrees, 1-360; convert to radians.
angle_deg = random.randint(1, 360)
angle_rad = math.radians(angle_deg)
print(f"POINT {oid}, ANGLE: {angle_deg}")
# Get x/y offset of random point on the imaginary circle around each point.
start_x = (point_spacing/2) * math.cos(angle_rad)
start_y = (point_spacing/2) * math.sin(angle_rad)
# Using the sample point coordinate as the starting point,
# calculate real-world x and y of both start and end of transect line.
start_point = arcpy.Point(x + start_x, y + start_y)
end_point = arcpy.Point(x - start_x, y - start_y)
# -S/-E suffixes differentiate start/edn boxes. Not currently used, but
# documented just in case (also, unique keys are required).
sample_box_dict[f'{oid}-S'] = create_box(box_rotation, start_point,
box_width, angle_deg, pnt_spat_ref)
sample_box_dict[f'{oid}-E'] = create_box(box_rotation, end_point,
box_width, angle_deg, pnt_spat_ref)
# Create a Polyline Geometry Object. Append dict entry to list.
new_line = arcpy.Polyline(arcpy.Array([start_point, end_point]))
transect_lines.append({'OID': oid, 'GEOM': new_line, 'ANGLE': angle_deg})
# Save start and end Point Geometry objects to another list.
line_enpoint_dict[oid] = (arcpy.PointGeometry(start_point),
arcpy.PointGeometry(end_point))
# Create three feature classes.
for fc, geom in [(fc_line, 'POLYLINE'), (fc_endpoints, 'POINT'), (fc_box, 'POLYGON')]:
arcpy.management.CreateFeatureclass(out_path=workspace, out_name=fc,
geometry_type=geom, spatial_reference=pnt_spat_ref)
print(f"CREATED FC {fc}")
# FC paths.
ln_fc = os.path.join(workspace, fc_line)
pt_fc = os.path.join(workspace, fc_endpoints)
box_fc = os.path.join(workspace, fc_box)
# Add fields.
for fc in (ln_fc, pt_fc, box_fc):
arcpy.management.AddField(fc, 'POINT_OID', 'SHORT')
arcpy.management.AddField(ln_fc, 'LINE_ANGLE', 'SHORT')
arcpy.management.AddField(pt_fc, 'POINT_TYPE', 'TEXT', field_length=5)
print("ALL FIELDS ADDED")
# Write transect line features.
with arcpy.da.InsertCursor(ln_fc, ['SHAPE@', 'POINT_OID', 'LINE_ANGLE']) as icurs:
for feature_dict in transect_lines:
icurs.insertRow([feature_dict['GEOM'],
feature_dict['OID'],
feature_dict['ANGLE']])
print("TRANSECTS WRITTEN")
# Write line start and end points.
with arcpy.da.InsertCursor(pt_fc, ['SHAPE@', 'POINT_OID', 'POINT_TYPE']) as icurs:
for oid, (start_geom, end_geom) in line_enpoint_dict.items():
icurs.insertRow([start_geom, oid, 'START'])
icurs.insertRow([end_geom, oid, 'END'])
print("ENDPOINTS WRITTEN")
# Create Sampling Zones----------------------------------#
# Transect Line Oriented
with arcpy.da.InsertCursor(box_fc, ['SHAPE@', 'POINT_OID']) as icurs:
for oid, geom in sample_box_dict.items():
icurs.insertRow([geom, oid.split('-')[0]])
print("BOXES WRITTEN")
Any chance you can us MS Paint or something similar to add a picture of what you are looking for? Not totally clear what you're after from your post. I'd be happy to take a crack at it if I get a better idea of what you are looking for.
Thanks for your response! Hopefully this image makes sense. The basic idea is that I have a grid of points that will remain constant and through each point I want to create lines of a fixed length that each have a unique, random angle. The endpoints of these lines will be the basis for 25mx25m survey plots. Maybe creating lines is the wrong approach but for the end result I want to have two points associated with each center-point that are at opposite, random angles and can be used to make survey plots.
In the image, the circles are the center-point grid and the squares are the plots.
I think your image makes sense, but let me know based on this output.
Here's another random run:
Don't know what your experience with Python is, but that's how I would accomplish this task. More than happy to clarify anything if you have questions.
Couple things:
#-----------------------------------------------------------------------------#
# IMPORTS
import math
import os
import random
import arcpy
#-----------------------------------------------------------------------------#
# INPUTS
# Full path to the geodatabase where the sampling points are, and where outputs will go.
workspace = r"your_geodatabase_path\GDB_NAME.gdb"
spat_ref_epsg = XXXX # This should be an integer code representing the EPSG number of a spatial reference.
# Name of the point feature class in the above geodatabase that houses the
# sampling points.
in_pnt_fc = "SAMPLE_POINT"
# Distance between sampling points (along x and y, not diagonals). Units match coordinate system.
point_spacing = 500 # feet, in my case
# Name of the output feature class to be created.
output_line_fc = "TRANSECT_LINE"
# -----------------------------------------------------------------------------#
# SETTINGS
# This just overwrites previous outputs to help with testing.
arcpy.env.overwriteOutput = False
#-----------------------------------------------------------------------------#
# MAIN
# Gets full path to the input sampling points feature class.
in_points = os.path.join(workspace, in_pnt_fc)
# NOT NEEDED. Compile the origin points here, for testing purposes mostly.
start_pts = []
"""
Using a SearchCursor, iterate over provided sample points.
At the given radius (half the sampling distance) create a random origin point for a line.
Get point on opposite side of the "circle" that is around the given sampling point.
Create a Polyline geometry object from those two points.
"""
transect_lines = []
with arcpy.da.SearchCursor(in_points, ['SHAPE@XY']) as scurs:
for (x, y), in scurs:
# print(f"POINT {oid}")
# Random angle in degrees, 1-360; convert to radians.
angle_deg = random.randint(1, 360)
# print(angle_deg)
angle_rad = math.radians(angle_deg)
# Get the x and y offset of the random point on the imaginary circle around each
# sampling point.
start_x = (point_spacing/2) * math.cos(angle_rad)
start_y = (point_spacing/2) * math.sin(angle_rad)
# Using the sample point coordinate as the starting point,
# calculate the real-world x and y of both the start and end of the transect line.
start_point = arcpy.Point(x + start_x, y + start_y)
end_point = arcpy.Point(x - start_x, y - start_y)
# Create a Polyline Geometry Object.
new_line = arcpy.Polyline(arcpy.Array([start_point, end_point]))
# Append the line geometry for this sampline point to a list.
transect_lines.append(new_line)
# NOT NEEDED. Testing the line origin points only.
start_pts.append(arcpy.PointGeometry(start_point))
print('\n')
# Create a blank polyline feature class to house the new transect lines.
ln_fc = arcpy.management.CreateFeatureclass(workspace, 'LINE', 'POLYLINE', spatial_reference=spat_ref_epsg)
with arcpy.da.InsertCursor(ln_fc, ['SHAPE@']) as icurs:
for geom in transect_lines:
icurs.insertRow([geom])
# NOT NEEDED. This creates a feature class for the origin points of the line.
pt_fc = arcpy.management.CreateFeatureclass(workspace, 'ORIGIN_POINT', 'POINT', spatial_reference=spat_ref_epsg)
with arcpy.da.InsertCursor(pt_fc, ['SHAPE@']) as icurs:
for geom in start_pts:
icurs.insertRow([geom])
Thank you for this documentation, it solves my problem perfectly! The only problem I have run into is that my ArcGIS Pro license doesn't include the "Feature Verticies to Points" tool, do you have any other suggestions for how to create those endpoints using other tools?
See another version below. This one will save your points too. It also creates your sampling boxes. There are four different methods you can use (see the "box_rotation" variable in the INPUTS section). Remember to update all the "INPUTS" based on your needs. Should be able to paste the code into the Python window again.
Only one version of the transect lines are shown (lines for the 'ANGLE_DIAMOND'-style boxes), but you can see all four styles of the sampling boxes. Some different sizes and distances from the origin points are also shown.
#-----------------------------------------------------------------------------#
# IMPORTS
import math
import os
import random
import sys
from typing import Literal, Union
import arcpy
#-----------------------------------------------------------------------------#
# INPUTS
# Full path to the geodatabase where the sampling points are, and where outputs will go.
workspace = r"C:TEMP\Random_Sample_Lines.gdb"
# Name of the point feature class in the above geodatabase that houses the
# sampling points.
in_pnt_fc = "_SAMPLE_POINT"
# Distance between sampling points (along x and y, not diagonals).
# Units match coordinate system.
point_spacing = 400 # feet, in my case
box_width = 50
box_rotation = 'ANGLE_BOX' # Literal[NORTH_BOX, NORTH_DIAMOND, ANGLE_BOX, ANGLE_DIAMOND]
# Name of the output feature class to be created.
fc_line = "TRANSECT_LINE"
fc_endpoints = "TRANSECT_LINE_ENDPOINTS"
# Leave this alone if you are not familiar with the f-string notation.
fc_box = f"TRANSECT_{box_rotation}"
# -----------------------------------------------------------------------------#
# SETTINGS
# This just overwrites previous outputs to help with testing.
arcpy.env.overwriteOutput = False
#-----------------------------------------------------------------------------#
# FUNCTIONS
def rotate_point(coordinate: tuple[int, int],
pivot_point: tuple[int, int],
angle: Union[int, float],
angle_format: Literal['DEGREES', 'RADIANS']='DEGREES'
) -> tuple[int, int]:
"""Rotates coordinate values around another point.
:coordinate: (x, y) coordinate pair to be rotated
:pivot_point: (x, y) coordinate pair around which point will be rotated
:angle: angle to rotate the coordinate
:angle_format: angle format, RADIANS or DEGREES
|(x_prime, y_prime)| rotated point coordinates"""
if angle_format == 'DEGREES':
angle = math.radians(angle)
x, y = coordinate
X, Y = pivot_point
x_trans = x - X
y_trans = y - Y
sin = math.sin(angle)
cos = math.cos(angle)
# Counter-clockwise rotation:
x_transprime = cos * x_trans - sin * y_trans
y_transprime = sin * x_trans + cos * y_trans
x_prime = x_transprime + X
y_prime = y_transprime + Y
# print(f"XY TRANS: {x_trans:,.0f}, {y_trans:,.0f}\n"
# f"XY TRANSPRIME: {x_transprime:,.0f}, {y_transprime:,.0f}\n"
# f"XY PRIME: {x_prime:,.0f}, {y_prime:,.0f}\n")
return (x_prime, y_prime)
def create_box(method: Literal['NORTH_BOX', 'NORTH_DIAMOND', 'ANGLE_BOX', 'ANGLE_DIAMOND'],
center_pnt: arcpy.Point,
width: Union[int, float],
angle: Union[int, float],
spat_ref: int) -> arcpy.Polygon:
"""Create a square box around a given point. Optionally rotate it.
:center_pnt: centroid of the square box
:width: width of side of box
:angle: angle the box should be rotated
|box_geom| geometry object representing the box"""
angle_dict = {'NORTH_BOX': 45,
'NORTH_DIAMOND': 0,
'ANGLE_BOX': angle + 45,
'ANGLE_DIAMOND': angle}
# Distance to move point based on known hypotenuse (1/2 width of box).
offset = math.sqrt(2 * ((width / 2) ** 2))
x, y = center_pnt.X, center_pnt.Y
Nx, Ny = x, y + offset
Ex, Ey = x + offset, y
Sx, Sy = x, y - offset
Wx, Wy = x - offset, y
box_corners = []
for coordinate in [(Nx, Ny), (Ex, Ey), (Sx, Sy), (Wx, Wy), (Nx, Ny)]:
box_corners.append(rotate_point(coordinate=coordinate,
pivot_point=(x, y),
angle=angle_dict[method],
angle_format='DEGREES'))
pnt_array = arcpy.Array([arcpy.Point(*corner) for corner in box_corners])
box_geom = arcpy.Polygon(pnt_array, spatial_reference=spat_ref)
return box_geom
#-----------------------------------------------------------------------------#
# MAIN
in_points = os.path.join(workspace, in_pnt_fc)
transect_lines = []
line_enpoint_dict = {}
sample_box_dict = {}
pnt_desc = arcpy.da.Describe(in_points)
pnt_oid_fld = pnt_desc['OIDFieldName']
pnt_spat_ref = pnt_desc['spatialReference'].factoryCode
with arcpy.da.SearchCursor(in_points, [pnt_oid_fld, 'SHAPE@XY']) as scurs:
for oid, (x, y) in scurs:
# Random angle in degrees, 1-360; convert to radians.
angle_deg = random.randint(1, 360)
angle_rad = math.radians(angle_deg)
print(f"POINT {oid}, ANGLE: {angle_deg}")
# Get x/y offset of random point on the imaginary circle around each point.
start_x = (point_spacing/2) * math.cos(angle_rad)
start_y = (point_spacing/2) * math.sin(angle_rad)
# Using the sample point coordinate as the starting point,
# calculate real-world x and y of both start and end of transect line.
start_point = arcpy.Point(x + start_x, y + start_y)
end_point = arcpy.Point(x - start_x, y - start_y)
# -S/-E suffixes differentiate start/edn boxes. Not currently used, but
# documented just in case (also, unique keys are required).
sample_box_dict[f'{oid}-S'] = create_box(box_rotation, start_point,
box_width, angle_deg, pnt_spat_ref)
sample_box_dict[f'{oid}-E'] = create_box(box_rotation, end_point,
box_width, angle_deg, pnt_spat_ref)
# Create a Polyline Geometry Object. Append dict entry to list.
new_line = arcpy.Polyline(arcpy.Array([start_point, end_point]))
transect_lines.append({'OID': oid, 'GEOM': new_line, 'ANGLE': angle_deg})
# Save start and end Point Geometry objects to another list.
line_enpoint_dict[oid] = (arcpy.PointGeometry(start_point),
arcpy.PointGeometry(end_point))
# Create three feature classes.
for fc, geom in [(fc_line, 'POLYLINE'), (fc_endpoints, 'POINT'), (fc_box, 'POLYGON')]:
arcpy.management.CreateFeatureclass(out_path=workspace, out_name=fc,
geometry_type=geom, spatial_reference=pnt_spat_ref)
print(f"CREATED FC {fc}")
# FC paths.
ln_fc = os.path.join(workspace, fc_line)
pt_fc = os.path.join(workspace, fc_endpoints)
box_fc = os.path.join(workspace, fc_box)
# Add fields.
for fc in (ln_fc, pt_fc, box_fc):
arcpy.management.AddField(fc, 'POINT_OID', 'SHORT')
arcpy.management.AddField(ln_fc, 'LINE_ANGLE', 'SHORT')
arcpy.management.AddField(pt_fc, 'POINT_TYPE', 'TEXT', field_length=5)
print("ALL FIELDS ADDED")
# Write transect line features.
with arcpy.da.InsertCursor(ln_fc, ['SHAPE@', 'POINT_OID', 'LINE_ANGLE']) as icurs:
for feature_dict in transect_lines:
icurs.insertRow([feature_dict['GEOM'],
feature_dict['OID'],
feature_dict['ANGLE']])
print("TRANSECTS WRITTEN")
# Write line start and end points.
with arcpy.da.InsertCursor(pt_fc, ['SHAPE@', 'POINT_OID', 'POINT_TYPE']) as icurs:
for oid, (start_geom, end_geom) in line_enpoint_dict.items():
icurs.insertRow([start_geom, oid, 'START'])
icurs.insertRow([end_geom, oid, 'END'])
print("ENDPOINTS WRITTEN")
# Create Sampling Zones----------------------------------#
# Transect Line Oriented
with arcpy.da.InsertCursor(box_fc, ['SHAPE@', 'POINT_OID']) as icurs:
for oid, geom in sample_box_dict.items():
icurs.insertRow([geom, oid.split('-')[0]])
print("BOXES WRITTEN")
This is fantastic, thank you again for your help, problem solved!
Thanks again for this solution! One additional output feature that we would like to create is a label for which directional quadrant (Northeast, Southwest, etc) that the "TRANSECT_LINE_ENDPOINTS" are in. This is defined in relation to the Origin point or Centerpoint, so naturally each pair of Endpoints will be in directionally opposite quadrants. The way we identify our plots currently is by using the Point_OID of the Endpoints and then manually labeling them with the quadrant when we record the associated data in the field. As shown in the attached example, if the Point_OID for these two Endpoints is 171, then we would record the associated plots as 171NW and 171SE.
What we're looking for is a way to automatically generate these quadrant codes in either the Point_OID field or some other field that we can use to identify plots. It seems to me that it should be a relatively simple to assign each code to an angle range (0-90 being NE for example), but I don't have the Python experience to write this code 🙂
I think I understand what you're after. I reworked some of the code internals to simplify things, so I'd try copying this version, rather than updating anything you're currently using (unless you're comfortable doing that). The outputs should be the same, with the addition of the new label field. I'm using a label expression to show the "point_id" + "quadrant" labels.
This takes the angle of the line (say, 45 degrees) and then assigns the quadrant (NE). It then inverts that to get the opposing quadrant (SE). The code is set up so that if the line is perfectly north-south or east-west (90 degrees, 180 degrees) then those labels will be "N" and "S" or "W" and "E", without the secondary direction. Let me know if something doesn't work!
Here's an W-E example, which I suppose would only occur... ~1% of the time?
"""
RANDOM ANGLE TRANSECT LINES CENTERED ON SAMPLING POINTS
"""
#-----------------------------------------------------------------------------#
# IMPORTS
import math
import os
import random
import sys
from typing import Literal, Union, Tuple
import arcpy
#-----------------------------------------------------------------------------#
# INPUTS
# Full path to the geodatabase where the sampling points are, and where outputs will go.
workspace = r"your_full_path\geodatabase_name.gdb"
# Name of the point feature class in the above geodatabase that houses the
# sampling points.
in_pnt_fc = "_SAMPLE_POINT"
# Distance between sampling points (along x and y, not diagonals).
# Units match coordinate system.
point_spacing = 400 # feet, in my case
box_width = 80
# NORTH_BOX: squares will be oriented to grid-north
# NORTH_DIAMOND: diamonds will be oriented to grid-north
# ANGLE_BOX: square sides oriented to the tranect lines
# ANGLE_DIAMOND: diamonds angled to the transect lines
box_rotation = 'ANGLE_BOX' # OPTIONS: (NORTH_BOX, NORTH_DIAMOND, ANGLE_BOX, ANGLE_DIAMOND)
# Name of the output feature class to be created.
fc_line = "TRANSECT_LINE"
fc_endpoints = "TRANSECT_LINE_ENDPOINTS"
fc_box = f"TRANSECT_{box_rotation}"
# -----------------------------------------------------------------------------#
# SETTINGS
# This just overwrites previous outputs to help with testing.
arcpy.env.overwriteOutput = True
#-----------------------------------------------------------------------------#
# FUNCTIONS
def rotate_point(coordinate: tuple[int, int],
pivot_point: tuple[int, int],
angle: Union[int, float],
angle_format: Literal['DEGREES', 'RADIANS']='DEGREES'
) -> tuple[float, float]:
"""Rotates coordinate values around another point.
:coordinate: (x, y) coordinate pair to be rotated
:pivot_point: (x, y) coordinate pair around which point will be rotated
:angle: angle to rotate the coordinate
:angle_format: angle format, RADIANS or DEGREES
|(x_prime, y_prime)| rotated point coordinates"""
if angle_format == 'DEGREES':
angle = math.radians(angle)
x, y = coordinate
X, Y = pivot_point
x_trans = x - X
y_trans = y - Y
sin = math.sin(angle)
cos = math.cos(angle)
# Counter-clockwise rotation:
x_transprime = cos * x_trans - sin * y_trans
y_transprime = sin * x_trans + cos * y_trans
x_prime = x_transprime + X
y_prime = y_transprime + Y
# print(f"XY TRANS: {x_trans:,.0f}, {y_trans:,.0f}\n"
# f"XY TRANSPRIME: {x_transprime:,.0f}, {y_transprime:,.0f}\n"
# f"XY PRIME: {x_prime:,.0f}, {y_prime:,.0f}\n")
return (x_prime, y_prime)
def create_box(method: Literal['NORTH_BOX', 'NORTH_DIAMOND', 'ANGLE_BOX', 'ANGLE_DIAMOND'],
center_pnt: arcpy.Point,
width: Union[int, float],
angle: Union[int, float],
spat_ref: int) -> arcpy.Polygon:
"""Create a square box around a given point. Optionally rotate it.
:center_pnt: centroid of the square box
:width: width of side of box
:angle: angle the box should be rotated
|box_geom| geometry object representing the box"""
angle_dict = {'NORTH_BOX': 45,
'NORTH_DIAMOND': 0,
'ANGLE_BOX': angle + 45,
'ANGLE_DIAMOND': angle}
# Distance to move point based on known hypotenuse (1/2 width of box).
offset = math.sqrt(2 * ((width / 2) ** 2))
x, y = center_pnt.X, center_pnt.Y
Nx, Ny = x, y + offset
Ex, Ey = x + offset, y
Sx, Sy = x, y - offset
Wx, Wy = x - offset, y
box_corners = []
for coordinate in [(Nx, Ny), (Ex, Ey), (Sx, Sy), (Wx, Wy), (Nx, Ny)]:
box_corners.append(rotate_point(coordinate=coordinate,
pivot_point=(x, y),
angle=angle_dict[method],
angle_format='DEGREES'))
pnt_array = arcpy.Array([arcpy.Point(*corner) for corner in box_corners])
box_geom = arcpy.Polygon(pnt_array, spatial_reference=spat_ref)
return box_geom
def get_angle_quadrant(a):
"""As previously defined, angle resctricted [1-360]."""
# Check for 90 degree cardinal directions.
if a == 90:
return "N"
elif a == 180:
return "W"
elif a == 270:
return "S"
elif a == 360:
return "E"
# Determine non-90 quadrants.
if 0 < a < 90:
return "NE"
elif 90 < a < 180:
return "NW"
elif 180 < a < 270:
return "SW"
elif 270 < a < 360:
return "SE"
else:
raise ValueError("Angle value must be numeric and between 1-360, inclusive.")
#-----------------------------------------------------------------------------#
# MAIN
in_points = os.path.join(workspace, in_pnt_fc)
transect_lines = []
sample_box_dict = {}
pnt_desc = arcpy.da.Describe(in_points)
pnt_oid_fld = pnt_desc['OIDFieldName']
pnt_spat_ref = pnt_desc['spatialReference'].factoryCode
with arcpy.da.SearchCursor(in_points, [pnt_oid_fld, 'SHAPE@XY']) as scurs:
for oid, (x, y) in scurs:
# Random angle in degrees, 1-360; convert to radians.
angle_deg = random.randint(1, 360)
angle_rad = math.radians(angle_deg)
print(f"POINT {oid}, ANGLE: {angle_deg}")
# Get x/y offset of random point on the imaginary circle around each point.
start_x = (point_spacing/2) * math.cos(angle_rad)
start_y = (point_spacing/2) * math.sin(angle_rad)
# Using the sample point coordinate as the starting point,
# calculate real-world x and y of both start and end of transect line.
start_point = arcpy.Point(x + start_x, y + start_y)
end_point = arcpy.Point(x - start_x, y - start_y)
# -S/-E suffixes differentiate start/end boxes. Not currently used, but
# documented just in case (also, unique keys are required).
sample_box_dict[f'{oid}-S'] = create_box(box_rotation, start_point,
box_width, angle_deg, pnt_spat_ref)
sample_box_dict[f'{oid}-E'] = create_box(box_rotation, end_point,
box_width, angle_deg, pnt_spat_ref)
# Create a Polyline Geometry Object. Append dict entry to list.
new_line = arcpy.Polyline(arcpy.Array([start_point, end_point]))
transect_lines.append({'OID': oid, 'GEOM': new_line, 'ANGLE': angle_deg})
# Create three feature classes.
for fc, geom in [(fc_line, 'POLYLINE'), (fc_endpoints, 'POINT'), (fc_box, 'POLYGON')]:
arcpy.management.CreateFeatureclass(out_path=workspace, out_name=fc,
geometry_type=geom, spatial_reference=pnt_spat_ref)
print(f"CREATED FC {fc}")
# FC paths.
ln_fc = os.path.join(workspace, fc_line)
pt_fc = os.path.join(workspace, fc_endpoints)
box_fc = os.path.join(workspace, fc_box)
# Add and ID field to the points, the lines, and box FCs.
for fc in (ln_fc, pt_fc, box_fc):
arcpy.management.AddField(fc, 'POINT_OID', 'SHORT')
# Add angle field to the line FC.
arcpy.management.AddField(ln_fc, 'LINE_ANGLE', 'SHORT')
# Add the point type and the point angle to the point FC.
arcpy.management.AddField(pt_fc, 'POINT_TYPE', 'TEXT', field_length=5)
arcpy.management.AddField(pt_fc, 'POINT_QUAD', 'TEXT', field_length=5)
print("ALL FIELDS ADDED")
# Write transect line features----------------------------------#
with arcpy.da.InsertCursor(ln_fc, ['SHAPE@', 'POINT_OID', 'LINE_ANGLE']) as icurs:
for feature_dict in transect_lines:
icurs.insertRow([feature_dict['GEOM'],
feature_dict['OID'],
feature_dict['ANGLE']])
print("TRANSECTS WRITTEN")
# Write endpoint features----------------------------------#
invert_quad = {"NE": "SW", "SE": "NW", "SW": "NE", "NW": "SE",
"N": "S", "S": "N", "E": "W", "W": "E", }
with arcpy.da.InsertCursor(pt_fc, ['SHAPE@', 'POINT_OID', 'POINT_TYPE', 'POINT_QUAD']) as icurs:
for feature_dict in transect_lines:
start_pnt = feature_dict['GEOM'].firstPoint
end_pnt = feature_dict['GEOM'].lastPoint
angle = feature_dict['ANGLE']
quadrant_1 = get_angle_quadrant(angle)
quadrant_2 = invert_quad[quadrant_1]
icurs.insertRow([start_pnt, feature_dict['OID'], "START", quadrant_1])
icurs.insertRow([end_pnt, feature_dict['OID'], "END", quadrant_2])
print("ENDPOINTS WRITTEN")
# Create Sampling Zones----------------------------------#
# Transect Line Oriented
with arcpy.da.InsertCursor(box_fc, ['SHAPE@', 'POINT_OID']) as icurs:
for oid, geom in sample_box_dict.items():
icurs.insertRow([geom, oid.split('-')[0]])
print("BOXES WRITTEN")
Hi Vince, I just tested it and it works great! Thanks again for your help!!