Creating point grid within a polygon, but aligned with polygon boundary

Question

I'm attempting to find an easier way to create a grid of 15m x 15m points within a specified area, but following the path of the polygon boundary. I have been creating points along a line, by distance, then using the array function and creating the appropriate rows and columns (and rotating them as necessary). However, this really only works well when I'm doing it within a square-shaped polygon. I think if I did it within a more complex shape, it wouldn't work (the points would get too close together/too far away at the turns). Something like the fishnet tool seems right, but the examples I see are always at 90°.

Also, if possible, I would like to give the points a specific label. Essentially starting at A1, B1, C1, etc. for the columns and continuing numerically for each row A1, A2, A3, etc. I currently have a function that does it for a single column, but it's quite time consuming to do it for each one.

I would appreciate any assistance on these matters. Thanks!

Matthew_Muehlhauser · Answer

Ok, I've come up with a slightly better version, although it currently assumes that it is a square / rectangle.

import arcpy
import math

def DegreestoRadians(deg):
    rad = deg*(math.pi/180.0)
    return rad

def getDistance(pnt1, pnt2):
    return math.sqrt((pnt2.X - pnt1.X)**2 + (pnt2.Y - pnt1.Y)**2)

def getAngleBetweenPoints(pnt1, pnt2):
    angle = math.degrees(math.atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X))
    if angle < 0:
        angle += 360
    return angle

def movePoint(pnt, angle, distance, spatial_reference):
    if isinstance(pnt, arcpy.PointGeometry):
        pnt = pnt.getPart(0)
    newX = pnt.X + math.cos(DegreestoRadians(float(angle))) * float(distance)
    newY = pnt.Y + math.sin(DegreestoRadians(float(angle))) * float(distance)
    new_point = arcpy.PointGeometry(arcpy.Point(newX, newY), spatial_reference)
    return new_point

project= arcpy.mp.ArcGISProject("CURRENT")
mv = project.activeMap
lyr = mv.listLayers("Test")[0]
with arcpy.da.SearchCursor(lyr, ["SHAPE@"]) as cursor:
    for row in cursor:
        poly = row[0]
sr = poly.spatialReference

# GET ANGLE
min_y = None
index = None

arr = poly.getPart(0)
n = len(arr)-1
arr.remove(n)

for i, pnt in enumerate(arr):
    if not min_y or pnt.Y < min_y:
        min_y = pnt.Y
        index = i

pnt1 = arr[index]
next = index+1 if index+1 < len(arr) else 0
prev = index-1

dist1 = getDistance(pnt1, arr[prev])
dist2 = getDistance(pnt1, arr[next])

closest_pnt = arr[prev] if dist1 < dist2 else arr[next]
furthest_pnt = arr[prev] if dist1 > dist2 else arr[next]

short_angle = getAngleBetweenPoints(pnt1, closest_pnt)
long_angle = getAngleBetweenPoints(pnt1, furthest_pnt)

short_pnts_n = int((getDistance(pnt1, closest_pnt) - 10) / 15)
long_pnts_n = int((getDistance(pnt1, furthest_pnt) - 10) / 15)

point_row = []
point_array = []

startPoint = movePoint(movePoint(pnt1, short_angle, 10, sr), long_angle, 10, sr)
point_row.append(startPoint)
point_array.append({'geom': startPoint, 'lbl': "A1"})

for i in range(short_pnts_n):
    new_point = movePoint(startPoint, short_angle, 15, sr)
    point_row.append(new_point)
    point_array.append({'geom': new_point, 'lbl': f'A{i+2}'})
    startPoint = new_point

alpha = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
i=0
while i < long_pnts_n:
    i += 1
    for j, pnt in enumerate(point_row):
        new_point = movePoint(pnt, long_angle, 15, sr)
        point_row[j] = new_point
        letter_1 = '' if i < 26 else alpha[math.floor(i/26)-1]
        letter_2 = alpha[i % 26]
        lbl = f'{letter_1}{letter_2}{j+1}'
        point_array.append({'geom': new_point, 'lbl': lbl})

with arcpy.da.InsertCursor('Test_Point', ["SHAPE@", "Label"]) as cursor:
    for ftr in point_array:
        new_ftr = [ftr["geom"], ftr["lbl"]]
        cursor.insertRow(new_ftr)

And this is the result:

Matthew_Muehlhauser · Answer

So, this would still require a little work to make easily repeatable, but I think the majority is there. Here is a potential solution in Python.

project= arcpy.mp.ArcGISProject("CURRENT")
mv = project.activeMap
lyr = mv.listLayers("Test")[0]
with arcpy.da.SearchCursor(lyr, ["SHAPE@"]) as cursor:
    for row in cursor:
        poly = row[0]

ext = poly.extent

point_array = []
alpha = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
i = 0
k = -1

x = ext.XMin
y = ext.yMin

while x < ext.XMax:
    j = 1
    while y < ext.YMax:
        new_point = arcpy.PointGeometry(arcpy.Point(x, y), poly.spatialReference)
        if k >= 0:
            label = f'{alpha[k]}{alpha[i]}{j}'
        else:
            label = f'{alpha[i]}{j}'
        if not poly.disjoint(new_point):
            point_array.append({'pnt': new_point, 'lbl': label})
        y += 15
        j += 1
    y = ext.YMin
    x += 15
    i += 1
    if i == len(alpha):
        i = 0
        k += 1

with arcpy.da.InsertCursor('Test_Point', ["SHAPE@", "Label"]) as cursor:
    for ftr in point_array:
        new_ftr = [ftr["pnt"], ftr["lbl"]]
        cursor.insertRow(new_ftr)

I do assume that the polygon that I'm working with is the last feature in the feature class, so actually getting the polygon shape attribute could be re-worked. I also assume that the spatial reference of this feature works in meters and not decimal degrees. And then the output feature was already created and only had a label field. This is what the output looks like though:

One thing I noticed though was that the first row and column were cut off by disjoint. So if you wanted to make sure that those were kept, you could add an initial displacement to the XMin / YMin values (eg. x = ext.XMin + 10).