2

I have a shapefile of the runways(polygons) of different airports:

runways digitalized

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: aerial view of the proyected runways horizontal view of the proyected runways

1
  • Which polygon are you talking about, I can only see lines on your images... Feb 14 at 19:06
0

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[0], y_lin_reg[0], 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
    brng = np.deg2rad(orientation)  # Bearing is degrees converted to radians.
    d = 9  # Distance in km

    lat1 = np.deg2rad(v1.y)  # Current lat point converted to radians
    lon1 = np.deg2rad(v1.x)  # Current long point converted to radians

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

    lat2 = np.rad2deg(lat2)
    lon2 = np.rad2deg(lon2)

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

#   creates the APP segment

    brng_app = np.deg2rad(or_to)
    d_app = 18  # Distance in km

    lat1_app = np.deg2rad(y_lin_reg[-1])  # Current lat point converted to radians
    lon1_app = np.deg2rad(x_lin_reg[-1])  # Current long point converted to radians

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

    lat2_app = np.rad2deg(lat2_app)
    lon2_app = np.rad2deg(lon2_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
runways = geopandas.read_file(runways.shp, crs={'init': 'epsg:3832'})

# reprojecting the dataset
runways = runways.to_crs("EPSG:4326")

# Compute the mid line endpoints of each polygon runway
runways['geometry'] = runways.apply(lambda row: stringify_runway(row.geometry, float(row.HETRUE_HDG), float(row.LETRUE_HDG)), axis=1)

# Write the new runway geometries 
runways.to_file(string_runways.shp)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.