Channel Tunnel/1989-Traverse.py

From OpenStreetMap Wiki
Jump to: navigation, search

See Channel Tunnel#1989 Traverse for context, and Channel Tunnel/1989-Traverse.gpx for generated output.

#!/usr/bin/env python
# Paul Sladen 2013-05-12
import math
# Korittke, N.; "Influence of horizontal Refraction on the traverse Measurements
# in Tunnels with Small Diameters"
# 
# Institute for Deposits and Surveying, DMT, Bochum, W-Germany
# http://www.slac.stanford.edu/econf/C9009106/papers/023.PDF
earth_radius = 6378137
bt = {'Great Britain 1989-03 Zig-Zag-Traverse': {
        # WGS84, corner of top of ramp, from Bing imagery/OSM mapping
        'initial': ('A2/1', 51.1064385,1.2782458),
        'data': (
        # Each ring-number is 1.5-metres in length
        ("A2/2",49.4450,0.372), 
        ("MTS2",56.5071,0.499), # Base of Adit A2
        ("MTS1",122.8067,0.655),  # Base of Adit A1
        ("17",119.81560,1.042),
        ("203",121.4786,1.302),
        ("345",120.8653,1.516),
        ("384",128.1011,1.575),
        ("423",119.5116,1.634),
        ("485",127.8503,1.728),
        ("548",123.1891,1.822),
        ("644",128.8531,1.967),
        ("691",124.4496,2.038),
        ("770",131.4425,2.158),
        ("852",128.4224,2.280),
        ("952",133.2877,2.431),
        ("1113",130.6273,2.673),
        ("1280",132.7052,2.924),
        ("1389",130.0523,3.088),
        ("1680",132.2658,3.526),
        ("1775",129.8136,3.669),
        ("1880",133.3176,3.827),
        ("2070",130.7208,4.113),
        ("2175",133.3090,4.271),
        ("2469",131.0453,4.713),
        ("2569",133.3729,4.863),
        ("2695",130.2317,5.052),
        ("2815",133.0986,5.232),
        ("3066",130.9558,5.609),
        ("3178",133.2237,5.777),
        ("3289",130.1029,5.944),
        ("3410",131.6857,6.126),
        ("3315",131.6927,6.284),
        ("3360",131.6781,6.487)
        )},
      'Great Britain 1989-12 Centre-Line-Traverse': {
        # WGS84, centre of top of ramp via Bing imagery/OSM mapping
        'initial': ('A2T', 51.106446,1.278199),
        'data': (
        # Each ring-number is 1.5-metres in length
        ("A2M",50.9911,0.496),
        ("ENT",120.5663,0.762),
        ("T5",120.7411,1.031),
        ("171",120.4814,1.255),
        ("296",121.4060,1.443),
        ("436",123.5396,1.654),
        ("568",125.7719,1.855),
        ("709",127.6997,2.066),
        ("859",130.1423,2.292),
        ("1018",131.6562,2.531),
        ("1294",131.6672,2.946),
        ("1548",131.6636,3.328),
        ("1827",131.6631,3.748),
        ("1982",131.6358,3.981),
        ("2095",131.6327,4.151),
        ("2234",131.6988,4.360),
        ("2369",131.5948,4.563),
        ("2498",131.6456,4.756),
        ("2622",131.6204,4.943),
        ("2868",131.6293,5.312),
        ("3121",131.6624,5.693),
        ("3270",131.6767,5.916),
        ("3410",131.6813,6.127),
        ("3682",131.6758,6.536),
        ("3850",131.6875,6.787),
        ("4036",131.6758,7.068),
        ("4297",131.6532,7.461),
        ("4441",129.5053,7.677),
        ("4584",126.1020,7.893),
        ("4728",127.2179,8.109),
        ("4868",130.6272,8.319),
        ("5008",133.9648,8.530),
        ("5139",137.1472,8.727),
        ("5279",136.4754,8.938),
        ("5419",133.1182,9.148),
        ("5546",131.6829,9.339),
        ("5722",131.6894,9.412),
        ("5857",131.6793,9.615),
        ("5994",132.5355,9.821),
        ("6136",134.6789,10.035),
        ("6268",136.8510,10.234),
        ("6412",139.0447,10.447),
        ("6558",141.3381,10.670),
        ("6705",143.7011,10.890),
        ("6846",145.9892,11.103),
        ("6985",147.6905,11.312),
        ("7187",147.8379,11.615),
        ("7392",147.8528,11.924),
        ("7574",147.8278,12.198),
        ("7725",147.8456,12.425),
        ("7903",147.8596,12.692),
        ("8050",147.8313,12.913),
        ("8222",147.8339,13.172),
        ("8396",147.8538,13.433),
        ("8544",147.8352,13.656))
        }
}

output = """<?xml version="1.0" encoding="UTF-8"?>
<gpx version="1.0">
<!--
     Paul Sladen, May 2013, for OpenStreetMap reference purposes

     Projected (using a simple spheroid model) into WGS84 from two
     estimated starting points at the top of Shakespeare Cliff Adit A2
     and then integrating over the two sets of Bearing+Distance pairs
     surveyed in 1989, and published in:

     Korittke, Norbert (1989)
     "Influence of horizontal Refraction on the traverse Measurements
     in Tunnels with Small Diameters"
     Deutsche Montan Technologie/Institute for Deposits and Surveying, Bochum, W-Germany
     http://www.slac.stanford.edu/econf/C9009106/papers/023.PDF p.13,14

     Tunnel chainages for the waypoint markers have been auto-generated from the
     stated marine service tunnel bore ring-number using an estimate of:

     chainage 19.841 + 0.0015062 * (ring_count - 1)
-->\n"""

number = 1
for name, survey in bt.items():
    bearing_traverse = survey['data']
    initial_latlon = survey['initial'][1:3]
    initial_name = survey['initial'][0]
    year_month = name[name.index('1989'):name.index('1989')+7]
    lat1, lon1 = math.radians(initial_latlon[0]), math.radians(initial_latlon[1])
    R = earth_radius
    path = [(initial_latlon[0], initial_latlon[1], initial_name)]
    output += """<wpt lat="%f" lon="%f"><name>%s</name></wpt>\n""" % path[-1]
    already_traversed = 0
    for ring, bearing, traverse in bearing_traverse:
        # surveyors use gradians (gon) 
        brng = math.radians(bearing*360/400)
        d = 1000 * (traverse - already_traversed)
        # based on http://www.movable-type.co.uk/scripts/latlong.html#destPoint
        # the original fails to mention that you need to work in radians everywhere
        # which is obvious, I guess
        lat2 = math.asin( math.sin(lat1)*math.cos(d/R) + math.cos(lat1)*math.sin(d/R)*math.cos(brng) )
        lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1), math.cos(d/R)-math.sin(lat1)*math.sin(lat2));
        try:
            chainage = str((int(ring)-1)*0.0015062+19.841)[:6].replace('.','+')
            ring = ring + ' Ch' + chainage
        except:
            pass
        path.append((math.degrees(lat2),math.degrees(lon2),ring))
        already_traversed = traverse
        lat1, lon1 = lat2, lon2
        output += """<wpt lat="%f" lon="%f"><name>%s</name></wpt>\n""" % path[-1]
    output += """<time>%s-31T00:00:00Z</time>
<trk><name>%s</name><number>%d</number><trkseg>\n""" % (year_month, name, number)
    for p in path:
        output += """<trkpt lat="%f" lon="%f"></trkpt>\n""" % p[0:2]
    output += """</trkseg></trk>\n"""
    number += 1
output += """</gpx>\n"""
print output