greatcircledistance_calc_from_osm_file.py

import math
from xml.dom import minidom

def coord(way_id):
    abs_file_path = "/Users/zia/Documents/Test/osm_road_length_data/itu_maslak.osm";
    xmldoc = minidom.parse(abs_file_path)
    way = xmldoc.getElementsByTagName("way")

    def node_coord(node_id):
        node = xmldoc.getElementsByTagName("node")
        latlon = list();
        for n in node:
            if int(n.getAttribute("id")) == node_id:
                latlon.append(float(n.getAttribute("lat")))
                latlon.append(float(n.getAttribute("lon")))
                return latlon

    coord_list = list();
    for w in way:
        if int(w.getAttribute("id")) == way_id:
            node_id = w.getElementsByTagName("nd")
            for nd in node_id:
                lat = node_coord(int(nd.getAttribute("ref")))[0]
                lon = node_coord(int(nd.getAttribute("ref")))[1]
                latlon = list()
                latlon.append(lat)
                latlon.append(lon)
                coord_list.append(latlon)
    return coord_list

def gcd(st, end):
    lat1 = st[0]
    long1 = st[1]
    lat2 = end[0]
    long2 = end[1]
    rLat1 = math.radians(lat1)
    rLong1 = math.radians(long1)
    rLat2 = math.radians(lat2)
    rLong2 = math.radians(long2)
    dLat = rLat2 - rLat1
    dLong = rLong2 - rLong1
    a = math.sin(dLat/2)**2 + math.cos(rLat1) * math.cos(rLat2) \
    * math.sin(dLong/2)**2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))

    theta = math.atan((6356752.3**2)*(math.tan((rLat1+rLat2)/2))/(6378137**2))
    rad_new = (((math.cos(theta))**2/(6378137**2))+((math.sin(theta))**2/(6356752.3**2)))**-0.5

    ret = list()
    ret.append(c)
    ret.append(rad_new)
    return ret

def distance(way_id):
    coord_list = list(coord(way_id));
    final_dist_wgs84 = 0.0;
    final_dist_localdatum = 0.0;
    rad_new = 0.0;
    for count in range(0,len(coord_list)-1):
        final_dist_wgs84 += (gcd(coord_list[count], coord_list[count+1])[0])*6371000;
        final_dist_localdatum += (gcd(coord_list[count], coord_list[count+1])[0])*(gcd(coord_list[count], coord_list[count+1])[1]);
        rad_new += gcd(coord_list[count], coord_list[count+1])[1]
    rad_new /= (len(coord_list)-1)
    return (6371000,final_dist_wgs84,rad_new,final_dist_localdatum)



print distance(20);