Hi @MeleKoneya,
The real bulk of of the lines takes place at line 77. When the Geodesic are created. Everything else I preparing to set the start and end locations of the flight path.
lines 77-109 is the mathematical expression to create the circle. Then print the start and end values
(line 78) num_points=100 it will create 100 vertices along the line. You can raise or lower this to smooth it out but I find 100 to be antiquate.
lines 138-146 set the height of the Geodesic line by changes the divisional value
lines 148-164 creates an array to hold the new z-values
lines 166-175 updates the z-vales of the line work to push the geodesic line work to 3D
I hope this helps. If what you are working in is publicly available I would love to see what it looks like. Reply with your site or an image if you'd rather.
Let me know if you need anything else.
Cheer,
Heath
Here is the example of my working site and script. Just for fun https://gis.appleton.org/grinch.html
import arcpy, random, webbrowser, math, datetime
import numpy as np
import math
import matplotlib.pyplot as plt
from IPython.display import display
from xml.sax.saxutils import escape
from arcpy import GetPortalInfo, GetActivePortalURL
#--- AGOL Login Info ---
agol_org = '[PORTAL AND AGOL URL]'
agol_user = '[USERNAME]' # See line 172 and make that match
agol_pass = '[PASSWORD]'
gis = arcpy.SignInToPortal(agol_org, agol_user, agol_pass)
print(arcpy.GetPortalInfo(portal_URL=arcpy.GetActivePortalURL()))
# Checks to see if the user is signed into Portal
##def portal_sign_in_check() -> None:
## portal_info: dict = GetPortalInfo(portal_URL = GetActivePortalURL())
## if portal_info['role'] == '':
## print('User is not signed-in to portal ' + str(portal_info))
## else:
## print('User is signed-in to portal ' + str(portal_info))
##
##print(portal_sign_in_check())
# get Grinch current location
grinchCurXY = "https://services1.arcgis.com/KowZ0VFM04w7hjgW/arcgis/rest/services/Grinch_Location/FeatureServer/0"
for loc in arcpy.da.SearchCursor(grinchCurXY, ["SHAPE@XY", "CITY_NAME", "POP", "CNTRY_NAME"]):
grXY = loc[0]
grCity = loc[1]
grPop = loc[2]
grCnt = loc[3]
print ("PointXY: {0}, City: {1}, Population: {2}, Country: {3}".format (loc[0], loc[1], loc[2], loc[3]))
# Update grinch location as historical
griHist = "https://services1.arcgis.com/KowZ0VFM04w7hjgW/arcgis/rest/services/grinchHistoricLocs/FeatureServer/0"
values = [(grXY, grCity, grPop, grCnt)]
griCur = arcpy.da.InsertCursor(griHist, ["SHAPE@XY", "CITY_NAME", "POP", "CNTRY_NAME"])
for loc in values:
griCur.insertRow(loc)
del griCur
# Updates historical points with xy flight data
arcpy.management.CalculateField(
in_table=griHist,
field="START_X",
expression="-123.105556",
expression_type="SQL",
code_block="",
field_type="TEXT",
enforce_domains="NO_ENFORCE_DOMAINS"
)
arcpy.management.CalculateField(
in_table=griHist,
field="START_Y",
expression="49.714444",
expression_type="SQL",
code_block="",
field_type="TEXT",
enforce_domains="NO_ENFORCE_DOMAINS"
)
arcpy.management.CalculateGeometryAttributes(
in_features=griHist,
geometry_property="POINT_X POINT_X;POINT_Y POINT_Y",
length_unit="",
area_unit="",
coordinate_system='GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]',
coordinate_format="SAME_AS_INPUT"
)
# Create a function to generate points along the geodesic
def great_circle_points(start, end, num_points=100):
lat1, lon1 = np.radians(start[1]), np.radians(start[0])
lat2, lon2 = np.radians(end[1]), np.radians(end[0])
d = 2 * np.arcsin(np.sqrt(np.sin((lat2 - lat1) / 2) ** 2 +
np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2))
lats = np.linspace(lat1, lat2, num_points)
lons = np.linspace(lon1, lon2, num_points)
arc_points = []
for i in range(num_points + 1):
a = math.sin((1 - i / num_points) * d) / math.sin(d)
b = math.sin((i / num_points) * d) / math.sin(d)
x = a * np.cos(lat1) * np.cos(lon1) + b * np.cos(lat2) * np.cos(lon2)
y = a * np.cos(lat1) * np.sin(lon1) + b * np.cos(lat2) * np.sin(lon2)
z = a * np.sin(lat1) + b * np.sin(lat2)
lat = np.arctan2(z, np.sqrt(x * x + y * y))
lon = np.arctan2(y, x)
arc_points.append((np.degrees(lon), np.degrees(lat)))
return arc_points
#set geodesic startup
for loc in arcpy.da.SearchCursor(griHist, ["SHAPE@XY", "CITY_NAME", "POP", "CNTRY_NAME", "START_X", "START_Y", "POINT_X", "POINT_Y"]):
start_point = (loc[4], loc[5])
end_point = (loc[6], loc[7])
print ("Start point: {0}, End point: {1}".format (start_point, end_point))
# calculate geodesic
# Generate geodesic points
arc_points = great_circle_points(start_point, end_point)
# Convert to 3D with Z values
line_length = arcpy.PointGeometry(arcpy.Point(*start_point)).distanceTo(arcpy.PointGeometry(arcpy.Point(*end_point)))
z_values = line_length * np.sin(np.linspace(0, math.pi, len(arc_points)))
arc_points_3d = [(lon, lat, z) for (lon, lat), z in zip(arc_points, z_values)]
# Use arcpy to create a polyline and save it to a feature class
arcpy.env.workspace = "[WORK SPACE ENVIRONMENT]"
arcpy.env.overwriteOutput = True
sr = arcpy.SpatialReference(4326) # WGS84
# Create the polyline geometry
polyline = arcpy.Polyline(arcpy.Array([arcpy.Point(*coords) for coords in arc_points_3d]), sr)
# Save the polyline to a feature class
arcpy.management.CreateFeatureclass(arcpy.env.workspace, "geodesicArc3D", "POLYLINE", spatial_reference=sr, has_z="ENABLED")
with arcpy.da.InsertCursor("geodesicArc3D", ["SHAPE@"]) as cursor:
cursor.insertRow([polyline])
#TEMPFlights = "holding"
TEMPFlights = "geodesicArc3D"
# Calculate the geodesic distance (length of the line)
with arcpy.da.SearchCursor(TEMPFlights, ["SHAPE@"]) as cursor:
for row in cursor:
line = row[0]
start_point = line.firstPoint
end_point = line.lastPoint
line_length = line.getLength("GEODESIC")/5
del cursor
# Function to adjust Z-values in a semi-circle arc
def adjust_z_values_semi_circle(line, max_z):
num_points = len(line.getPart(0))
angles = np.linspace(0, np.pi, num_points)
z_values = max_z * np.sin(angles)
# Print the Z-values
print("Z-values for each vertex:")
for i, z in enumerate(z_values):
print(f"Vertex {i + 1}: Z = {z}")
# Create new vertices with updated Z-values
new_vertices = []
for i, point in enumerate(line.getPart(0)):
new_vertices.append(arcpy.Point(point.X, point.Y, z_values[i]))
return arcpy.Array(new_vertices), z_values
# Extract all vertices and update their Z-values
with arcpy.da.UpdateCursor(TEMPFlights, ["SHAPE@"]) as cursor:
for row in cursor:
line = row[0]
new_vertices, z_values = adjust_z_values_semi_circle(line, line_length)
new_line = arcpy.Polyline(arcpy.Array([new_vertices]), line.spatialReference, has_z=True)
row[0] = new_line
cursor.updateRow(row)
del cursor
print("Z-values have been updated successfully.")
# Plot the Z-values
'''
plt.plot(z_values)
plt.title("Z-values Along the Arc (Semi-Circle Shape)")
plt.xlabel("Vertex Index")
plt.ylabel("Z Value (meters)")
plt.grid(True)
plt.show()
'''
# Grinch Flight lines
grinLine = "https://services1.arcgis.com/KowZ0VFM04w7hjgW/arcgis/rest/services/flightpath/FeatureServer/0"
arcpy.management.Append(
inputs="geodesicArc3D",
target=grinLine,
schema_type="NO_TEST",
field_mapping='POINT_X "POINT_X" true true false 0 Double 0 0,First,#;POINT_Y "POINT_Y" true true false 0 Double 0 0,First,#;START_X "START_X" true true false 0 Double 0 0,First,#;START_Y "START_Y" true true false 0 Double 0 0,First,#',
subtype="",
expression="",
match_fields=None,
update_geometry="NOT_UPDATE_GEOMETRY"
)
print ("done")
