# Create 3D lines at a specific distance and bearing from a 2D polygon using Shapely and Geopandas

I have a shapefile of the runways(polygons) of different airports: I want to build a line that represents the 3D trajectory followed by an airplane 9 km after take off and 18 km before landing and 1000 meters above ground.

My questions:

1. How to calculate the coordinates of the starting points?
2. How to calculate the coordinates of the target 3D points?
3. How to use shapely to create the 3D Multilines

I'm trying to solve this as follows:

1. Creating a line that goes through the polygon exactly in the middle, obtain the coordinates of the vertices of this line and supply them as the coordinates of the known point, making a lineal regression it will achieved it.

2. Calculate the coordinates of a point 9 kilometers away from the runway and keeping the orientation of the runway, that is, it will be a 3D line projected from the runway, the same to represent the landing.

3. Use the equation of the direct problem in geodesy, to calculate the coordinates of the points 9km and 18km away from the runway

The output I would like to achieve, is like that:  • Which polygon are you talking about, I can only see lines on your images... Feb 14 at 19:06

After several attempts, I have achieved the expected result in the following way, however if someone knows a more pythonic way it will be too well received:

``````python
# Compute the mid line endpoints of a polygon runway
# the parameters that this function receives come from my *runways.shp* shapefile:
# the geometry (polygon) of a shapefile read with geopandas (runways)
# the orientation of the runway head measured from north to south
# the orientation of the runway tail measured from north to south
def stringify_runway(polygon_runway, orientation, or_to):

# Extract geographic coordinates of polygon vertex
(x, y) = polygon_runway.exterior.coords.xy

# Check whether the best geographical precision can be achieved with a x-> or y-> linear interpolation
x_first = (orientation + 45) % 180 > 90
if x_first:
# Use linear regression to find the runway's center equation
model = np.polyfit(x[0:-1], y[0:-1], 1)
predict = np.poly1d(model)
# Thanks to the linear equation, compute the coordinates of the runway's endppoints
x_lin_reg = [min(x), max(x)]
y_lin_reg = predict(x_lin_reg)
else:
# Use linear regression to find the runway's center equation
model = np.polyfit(y[0:-1], x[0:-1], 1)
predict = np.poly1d(model)
# Thanks to the linear equation, compute the coordinates of the runway's endppoints
y_lin_reg = [min(y), max(y)]
x_lin_reg = predict(y_lin_reg)

# Creates the runway segment
r = LineString([[x_lin_reg, y_lin_reg, 0.0], [x_lin_reg[-1], y_lin_reg[-1], 0.0]])

# Creates the TakeOff trajectorie
# knowing that the planes use approximately 60% of the runway to take off
v1 = r.interpolate(0.4, normalized=True)
runway = LineString([[v1.x, v1.y, 0.0], [x_lin_reg[-1], y_lin_reg[-1], 0.0]])

# Direct problem equation in Geodesy, for calculated the coord of the proyected points
# lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
# lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))

R = 6378.1  # Radius of the Earth km
d = 9  # Distance in km

lat2 = np.arcsin(np.sin(lat1)*np.cos(d/R) + np.cos(lat1)*np.sin(d/R)*np.cos(brng))
lon2 = lon1 + np.arctan2(np.sin(brng)*np.sin(d/R)*np.cos(lat1), np.cos(d/R)-np.sin(lat1)*np.sin(lat2))

t = LineString([[v1.x, v1.y, 0.0], [lon2, lat2, 1000.0]])

#   creates the APP segment

d_app = 18  # Distance in km

lat2_app = np.arcsin(np.sin(lat1_app)*np.cos(d_app/R) + np.cos(lat1_app)*np.sin(d_app/R)*np.cos(brng_app))
lon2_app = lon1_app + np.arctan2(np.sin(brng_app)*np.sin(d_app/R)*np.cos(lat1_app), np.cos(d_app/R)-np.sin(lat1_app)*np.sin(lat2_app))

a = LineString([[x_lin_reg[-1], y_lin_reg[-1], 0.0], [lon2_app, lat2_app, 1000.0]])

#     agrouping segments
lines = [a, runway, t]
string_runway = MultiLineString(lines)

return string_runway
``````

calling up the function:

``````python