Select to view content in your preferred language

Creating randomly angled lines

1415
9
Jump to solution
11-11-2022 09:58 AM
Sam44Moles
Emerging Contributor

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.

0 Kudos
1 Solution

Accepted Solutions
VinceE
by
Frequent Contributor

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.

VinceE_1-1669060765014.png

 

#-----------------------------------------------------------------------------#
# 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")

 

View solution in original post

0 Kudos
9 Replies
VinceE
by
Frequent Contributor

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.

0 Kudos
Sam44Moles
Emerging Contributor

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.

canvas.png

0 Kudos
VinceE
by
Frequent Contributor

I think your image makes sense, but let me know based on this output.

VinceE_1-1668744871065.png

Here's another random run:

VinceE_2-1668744944693.png

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:

  • You can run this in an IDE with Pro's built-in Python environment, or simply paste the whole thing in the Python window if that is easier for you.
  • You need to update the five values in the "inputs" section so they are tailored to your situation.
  • I would use projected coordinate systems for all of this.
  • If two adjacent north-south points also had two perfectly north-south transect lines, you could have overlapping endpoints. There are checks for this, but if it was me, I would just re-run the script to try again. Currently the full length of a line is the same as the length between any two adjacent (along x-y, not diagonal) grid points. If you have other specifications, feel free to add that.
  • This script also creates the origin points for these lines (see small points). This was for testing (the "NOT NEEDED" comments in the code below show those parts). It does not create the endpoints.
  • To create your endpoints at each end of the lines (since it sounds like you need those), I would just use the built-in Feature Vertices to Points

 

 

#-----------------------------------------------------------------------------#
# 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])

 

0 Kudos
Sam44Moles
Emerging Contributor

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?

0 Kudos
VinceE
by
Frequent Contributor

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.

VinceE_1-1669060765014.png

 

#-----------------------------------------------------------------------------#
# 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")

 

0 Kudos
Sam44Moles
Emerging Contributor

This is fantastic, thank you again for your help, problem solved!

Sam44Moles
Emerging Contributor

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 🙂

0 Kudos
VinceE
by
Frequent Contributor

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 (SW). 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!
VinceE_0-1701115619227.png

Here's an W-E example, which I suppose would only occur... ~1% of the time?

VinceE_0-1701116510429.png

 

 

"""
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")

 

 

0 Kudos
Sam44Moles
Emerging Contributor

Hi Vince, I just tested it and it works great! Thanks again for your help!!