electronic Terrain and Obstacle Data, Salzburg, Austria

Table of Contents

• Introduction to electronic Terrain and Obstacle Data
• Description of Data
• Flight Safety Zone Modeling
• Approach
• Python Implementation
  • Data Preparation
  • Model Creation
  • Point Modeling
  • Polygon Construction
  • Obstacle Detection
• Result
• References

Introduction to electronic Terrain and Obstacle Data

electronic terrain and obstacle data is used to maintain safety margins and enhance ground proximity warning for the safety of flights. Obstacle data represents individual points, lines or contours indicating horizontal and vertical extent of a structure.

Description of Data

The study area is Salzburg Airport. The obstacle data is in the form of vectors that shows the shape of the obstacles in the area and may also stores the heights of the obstacles. Data required are the location of the airport's runway, Digital Terrain Model of the Area, and the Digital Surface Model of the Area. First, a geometric model is created based on the location of the airport's runway. from the plane that intersects with the runway, we created two planes on the direction of the arrival and departure of the aircraft to and from the runway. The arrival/departure plane has 2% inclination to and from the runway.

Flight Safety Zone Modeling

Figure 1: Obstacle Detection Area.

Approach

The obstacles data can be obtained by using intersection of a raster which modeled the safety zone and with the DSM. The safety zone are modeled as vector data which then rasterized. after the intersection with DSM has been done, the raster file will be vectorized and can be used as obstacle data.

Python Implementation

Data Preparation

The first part of the modeling is the preparation of the points that are being used in the safezone planes. starting with the two points that corresponds to the runway, we use height from the points or use height information from the DTM. for this purpose we use opalsImport to get the points. If we decide to use DTM as points' heights then we use Module AddInfo to add height information from grid file (DTM) into the ODM z component

    imp = Import.Import()
    imp.inFile = axis_file
    imp.outFile = axisODM
    imp.run()
    imp = None
 
    #
    # get height from DTM
    #
    if not point_height:
        logger.verbose("Loading point height from DTM")
        addin = AddInfo.AddInfo()
        addin.inFile = axisODM
        addin.gridFile = dtm_file
        addin.attribute = "_z=r[0]"
        addin.run()
        addin = None
    
    inDM = DM.Datamanager.load(axisODM, True, False)
    lf = DM.AddInfoLayoutFactory()
    lf.addColumn(DM.ColumnType.double_, "_z")
    layout = lf.getLayout()
    
    apts = []
    for p in inDM.points(layout):
        if point_height:
            apts.append(np.array([p.x, p.y, p.z]))
        else:
            apts.append(np.array([p.x, p.y, p.info().get(0)]))
            if abs(p.info().get(0) - p.z):
                print (
                    "point height : {}\ndtm height : {}".format(
                        p.z, p.info().get(0)))
            ##
            ## height_difference_threshold = 5
            ##
            # if height difference (p.info().get(0)-p.z) > height_difference_threshold
            ##
    
    inCRS = inDM.getCRS()
    logger.verbose("CRS\t"+inCRS)
    inDM = None
    
    pts = np.vstack(apts)

Model Creation

We use two points from axis file as starting points to refer to both end of the airport's runway

    A = pts[0, :]
    B = pts[1, :]
 
    #
    # airport reference point used for building circular plane zone
    # it is assumed that A and B has same height
    #
    mid = 0.5 * (A + B)
    height_ref = mid[2]

and calculate direction vectors to compute safety zone vertices

    #
    # calculate direction vector (AB)/|AB|
    #
    strip_length_3d = np.linalg.norm(B - A)
    strip_length_2d = np.linalg.norm(B[:2] - A[:2])
    
    dir = (B[:2] - A[:2]) / strip_length_2d
    
    dir2d = np.array([dir[0], dir[1], 0])
    dir3d = (B - A) / strip_length_3d  # 3d
    
    dir_perp = np.array([dir[1], -dir[0]])
    dir_perp_3d = np.array([dir[1], -dir[0], 0])

here are the constant parameters used for constructing safety zone vertices

    parallel_offset = 60.
    perp_offset = 150
    side_ramp_offset = 850.
 
    far_ramp_parallel_length = 10000.
    far_ramp_perp_length = 1500.
 
    far_ramp_height = 120.
    side_ramp_height = 10.2
 
    parallel_inclination = 1.2 / 100.
    perp_inclination = 1.2 / 100.
 
    side_ramp_height_difference = side_ramp_offset * perp_inclination
    far_ramp_parallel_height_difference = far_ramp_parallel_length * parallel_inclination
 
    max_far_ramp_height_difference = 100.
    max_side_height_difference = 45.
    side_radius_flat = 5100.
    side_radius_ramp = 5100.

Point Modeling

we start by constructing extension of the starting points outward 30 meters each

    #
    # create line from two point AB. call this AB Line.
    # extend the line AB by 60m
    #
    A_ = A - dir3d * parallel_offset
    B_ = B + dir3d * parallel_offset

then we compute points perpendicular to the runway

    #
    # from each endpoint of this line, create a line perpendicular of length 300m
    # (150m from the line's endpoint)
    #
    A_L = A_ - dir_perp_3d * perp_offset
    A_R = A_ + dir_perp_3d * perp_offset
    B_L = B_ - dir_perp_3d * perp_offset
    B_R = B_ + dir_perp_3d * perp_offset

then we compute points in the tilted plane on each side of the runway

    #
    # from pair of points in the left and right side create a tilted plane with 1.2% slope inclination
    # in the direction perpendicular to the line AB
    #
    A_LL = A_L - dir_perp_3d * side_ramp_offset + \
        np.array([0, 0, side_ramp_height_difference])
    A_RR = A_R + dir_perp_3d * side_ramp_offset + \
        np.array([0, 0, side_ramp_height_difference])
    B_LL = B_L - dir_perp_3d * side_ramp_offset + \
        np.array([0, 0, side_ramp_height_difference])
    B_RR = B_R + dir_perp_3d * side_ramp_offset + \
        np.array([0, 0, side_ramp_height_difference])

we continue with computing tilted plane in the direction of aircraft departure and arrival

    #
    # from pair of points in the front and back side create a tilted plane with 1.2% slope
    # inclination and 15% width expansion.
    #
    #
    A_L1 = A_L - dir2d * far_ramp_parallel_length - dir_perp_3d * \
        far_ramp_perp_length + np.array([0, 0, far_ramp_parallel_height_difference])
    A_R1 = A_R - dir2d * far_ramp_parallel_length + dir_perp_3d * \
        far_ramp_perp_length + np.array([0, 0, far_ramp_parallel_height_difference])
    B_L1 = B_L + dir2d * far_ramp_parallel_length - dir_perp_3d * \
        far_ramp_perp_length + np.array([0, 0, far_ramp_parallel_height_difference])
    B_R1 = B_R + dir2d * far_ramp_parallel_length + dir_perp_3d * \
        far_ramp_perp_length + np.array([0, 0, far_ramp_parallel_height_difference])

to fill in the hole gap between tilted plane on the side and tilted plane in the departure/arrival, we compute the intersection points

    #
    # compute intersection of side point to the front/rear tilted plane
    #
    A_L0, al0ratio = get_intersect(A_L, A_L1, A_LL)
    A_R0, ar0ratio = get_intersect(A_R, A_R1, A_RR)
    B_L0, bl0ratio = get_intersect(B_L, B_L1, B_LL)
    B_R0, br0ratio = get_intersect(B_R, B_R1, B_RR)

Polygon Construction

    #
    # construct 9 planes
    #
    dm = DM.Datamanager.create(inFile, False)
    dm.setCRS(inCRS)
    f = DM.PolygonFactory()
 
    if draw_runway_plane:
        # runway plane
        # A_L A_R B_R B_L
        f.addPoint(*A_L)
        f.addPoint(*A_R)
        f.addPoint(*B_R)
        f.addPoint(*B_L)
        f.closePart()
        p1 = f.getPolygon()
        dm.addPolygon(p1)
 
    # runway leftside ramp
    # A_L B_L B_LL A_LL
    f.addPoint(*B_LL)
    f.addPoint(*A_LL)
    f.addPoint(*A_L)
    f.addPoint(*B_L)
    f.closePart()
    p2 = f.getPolygon()
    dm.addPolygon(p2)
 
    # runway rightside ramp
    # A_R A_RR B_RR B_R
    f = DM.PolygonFactory()
    f.addPoint(*A_R)
    f.addPoint(*A_RR)
    f.addPoint(*B_RR)
    f.addPoint(*B_R)
    f.closePart()
    p3 = f.getPolygon()
    dm.addPolygon(p3)
 
    # runway rearside ramp
    # A_L A_L1 A_R1 A_R
    f = DM.PolygonFactory()
    f.addPoint(*A_L)
    f.addPoint(*A_L1)
    f.addPoint(*A_R1)
    f.addPoint(*A_R)
    f.closePart()
    p4 = f.getPolygon(DM.Orientation.ccw)
    dm.addPolygon(p4)
 
    # runway frontside ramp
    # B_R B_R1 B_L1 B_L
    f.addPoint(*B_R)
    f.addPoint(*B_R1)
    f.addPoint(*B_L1)
    f.addPoint(*B_L)
    f.closePart()
    p5 = f.getPolygon()
    dm.addPolygon(p5)
 
    # rear-left triangle ramp
    # A_L A_LL A_L0
    f.addPoint(*A_L)
    f.addPoint(*A_LL)
    f.addPoint(*A_L0)
    f.closePart()
    p6 = f.getPolygon()
    dm.addPolygon(p6)
 
    # rear-right triangle ramp
    # A_R A_R0 A_RR
    f.addPoint(*A_R)
    f.addPoint(*A_R0)
    f.addPoint(*A_RR)
    f.closePart()
    p7 = f.getPolygon()
    dm.addPolygon(p7)
 
    # front-left triangle ramp
    # B_L B_L0 B_LL
    f.addPoint(*B_L)
    f.addPoint(*B_L0)
    f.addPoint(*B_LL)
    f.closePart()
    p8 = f.getPolygon()
    dm.addPolygon(p8)
 
    # front-right triangle ramp
    # B_R B_RR B_R0
    f.addPoint(*B_R)
    f.addPoint(*B_RR)
    f.addPoint(*B_R0)
    f.closePart()
    p9 = f.getPolygon()
    dm.addPolygon(p9)

Obstacle Detection

Obstacles are detected by intersecting rasterized safety zone model with DSM and DTM

    #
    # intersect with DSM and create difference mask
    #
    alg.reset()
    alg.inFile = [step2File, dsm_file]
    alg.outFile = intersectFile
    alg.formula = "return r[1]-r[0] if r[0] is not None and r[1]-r[0]>1 else None"
    alg.limit = "dataset:0"
    alg.run()
    
    #
    # intersect with DTM and create difference mask
    #
    alg.reset()
    alg.inFile = [step2File, dtm_file]
    alg.outFile = intersectDtmFile
    alg.formula = "return r[1]-r[0] if r[0] is not None and r[1]-r[0]>1 else None"
    alg.limit = "dataset:0"
    alg.run()

Since opalsContouring only supports 8-bit raster file, The intersected file needs to be converted. This can be done using opalsShade which can also used as Visualization.

    #
    # build 8-bit version of the mask
    #
    shd = Shade.Shade()
    shd.inFile = intersectFile
    shd.outFile = intersectShade
    shd.run()
 
    from opals import Contouring
    cnt = Contouring.Contouring()
    cnt.inFile = intersectShade
    cnt.outFile = obstacleFile
    cnt.run()

Result

Figure 2: Output of obstacle detection on DTM.

    #
    # create outline model of the obstacle field
    #
    shd.reset()
    shd.inFile = step2File
    shd.sunPosition = (350, 95)
    shd.run()
 
    cnt.reset()
    cnt.inFile = step2Shade
    cnt.outFile = regionFile
    cnt.run()
 
    opals_step_log()  # 5
 
    # visualize temporary model
    shd.reset()
    shd.inFile = step1File
    shd.sunPosition = (350, 95)
    shd.run()

References

• TO COME!