← Back to branch summary

~ubuntu-branches/ubuntu/trusty/qiime/trusty

« back to all changes in this revision

Viewing changes to qiime/distance_matrix_from_mapping.py

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2013-06-17 18:28:26 UTC
  • mfrom: (9.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20130617182826-376az5ad080a0sfe
Tags: 1.7.0+dfsg-1
Upload preparations done for BioLinux to Debian

Show diffs side-by-side

added added

removed removed

Lines of Context:
qiime/distance_matrix_from_mapping.py
5
5
 
6
6
__author__ = "Antonio Gonzalez Pena"
7
7
__copyright__ = "Copyright 2011, The QIIME Project"
8
 
__credits__ = ["Antonio Gonzalez Pena"]
 
8
__credits__ = ["Antonio Gonzalez Pena", "Andrew J. King"]
9
9
__license__ = "GPL"
10
 
__version__ = "1.5.0"
 
10
__version__ = "1.7.0"
11
11
__maintainer__ = "Antonio Gonzalez Pena"
12
12
__email__ = "antgonza@gmail.com"
13
13
__status__ = "Release"
27
27
"""
28
28
 
29
29
from qiime.util import FunctionWithParams
30
 
from qiime.parse import parse_mapping_file_to_dict
31
 
from qiime.format import format_distance_matrix
32
 
from numpy import array, reshape
33
 
from sys import exit
 
30
from numpy import array, reshape, nan, radians, zeros
 
31
from math import atan, tan, sin, cos, pi, sqrt, atan2, acos, asin
34
32
 
35
 
def distance_matrix(input_path, column):
 
33
def compute_distance_matrix_from_metadata(column_data):
36
34
    """ calculates distance matrix on a single column of a mapping file
37
35
    
38
36
    inputs:
39
 
     input_path (file handler)
40
 
     column (str)
41
 
    """
42
 
    data, comments = parse_mapping_file_to_dict(input_path)
43
 
    column_data = []
44
 
    column_headers = []
45
 
    for i in data:
46
 
        if column not in data[i]:
47
 
            stderr.write("\n\nNo column: '%s' in the mapping file. Existing columns are: %s\n\n" % (column,data[i].keys()))
48
 
            exit(1)
49
 
        try:
50
 
            column_data.append(float(data[i][column]))
51
 
        except ValueError:
52
 
            stderr.write("\n\nall the values in the column '%s' must be numeric but '%s' has '%s'\n\n"\
53
 
                % (column,i,data[i][column]))
54
 
            exit(1)
55
 
            
56
 
        column_headers.append(i)
57
 
    
 
37
     column_data (list of values)
 
38
    """     
58
39
    data_row = array(column_data)
59
40
    data_col = reshape(data_row, (1, len(data_row)))
60
41
    dist_mtx = abs(data_row-data_col.T)
61
42
    
62
 
    return format_distance_matrix(column_headers, dist_mtx)
 
 
b'\\ No newline at end of file'
 
43
    return dist_mtx
 
44
    
 
45
    
 
46
def dist_vincenty(lat1, lon1, lat2, lon2, iterations=20):
 
47
    """Returns distance in meters between two lat long points
 
48
       
 
49
       Vincenty's formula is accurate to within 0.5mm, or 0.000015" (!),
 
50
       on the ellipsoid being used. Calculations based on a spherical model,
 
51
       such as the (much simpler) Haversine, are accurate to around 0.3%
 
52
       (which is still good enough for most purposes, of course).
 
53
       from: http://www.movable-type.co.uk/scripts/latlong-vincenty.html
 
54
 
 
55
       Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
 
56
       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975    */
 
57
       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf  
 
58
              
 
59
       This code was modified from geopy and movable-type:
 
60
       http://code.google.com/p/geopy/source/browse/trunk/geopy/distance.py?r=105
 
61
       http://www.movable-type.co.uk/scripts/latlong-vincenty.html
 
62
    """
 
63
    if lat1<-90 or lat1>90 or lat2<-90 or lat2>90 or lon1<-180 or lon1>180 or lon2<-180 or lon2>180:
 
64
        raise ValueError, "Latitude values shoulds range from (-90,90) and longitude from (-180,180) but one of the input values is out of bounds. Latitude_1: %f, Logitude_1: %f, Latitude_2: %f, Logitude_2: %f" % (lat1, lon1, lat2, lon2)
 
65
    
 
66
    major, minor, f = 6378137, 6356752.314245, 1/298.257223563
 
67
    
 
68
    lat1, lng1, lat2, lng2 = radians(lat1), radians(lon1), radians(lat2), radians(lon2)
 
69
    delta_lng = lng2 - lng1
 
70
    reduced_lat1, reduced_lat2 = atan((1 - f) * tan(lat1)), atan((1 - f) * tan(lat2))
 
71
 
 
72
    sin_reduced1, cos_reduced1 = sin(reduced_lat1), cos(reduced_lat1)
 
73
    sin_reduced2, cos_reduced2 = sin(reduced_lat2), cos(reduced_lat2)
 
74
 
 
75
    lambda_lng = delta_lng
 
76
    lambda_prime = 2 * pi
 
77
    while abs(lambda_lng - lambda_prime) > 10e-12 and iterations > 0:
 
78
        sin_lambda_lng, cos_lambda_lng = sin(lambda_lng), cos(lambda_lng)
 
79
 
 
80
        sin_sigma = sqrt(
 
81
            (cos_reduced2 * sin_lambda_lng) ** 2 +
 
82
            (cos_reduced1 * sin_reduced2 -
 
83
            sin_reduced1 * cos_reduced2 * cos_lambda_lng) ** 2
 
84
        )
 
85
        if sin_sigma == 0:
 
86
            return 0 # Coincident points
 
87
 
 
88
        cos_sigma = (
 
89
            sin_reduced1 * sin_reduced2 +
 
90
            cos_reduced1 * cos_reduced2 * cos_lambda_lng
 
91
        )
 
92
        sigma = atan2(sin_sigma, cos_sigma)
 
93
 
 
94
        sin_alpha = ( cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma )
 
95
        cos_sq_alpha = 1 - sin_alpha ** 2
 
96
 
 
97
        if cos_sq_alpha != 0:
 
98
            cos2_sigma_m = cos_sigma - 2 * ( sin_reduced1 * sin_reduced2 / cos_sq_alpha )
 
99
        else:
 
100
            cos2_sigma_m = 0.0 # Equatorial line
 
101
 
 
102
        C = f / 16. * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))
 
103
 
 
104
        lambda_prime = lambda_lng
 
105
        lambda_lng = (
 
106
            delta_lng + (1 - C) * f * sin_alpha * (
 
107
               sigma + C * sin_sigma * ( 
 
108
                   cos2_sigma_m + C * cos_sigma * ( -1 + 2 * cos2_sigma_m ** 2 )
 
109
               )
 
110
            )
 
111
        )
 
112
        iterations -= 1
 
113
 
 
114
    if iterations == 0:
 
115
        raise ValueError("Vincenty formula failed to converge!")
 
116
 
 
117
    u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2
 
118
    A = 1 + u_sq / 16384. * ( 4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq)) )
 
119
    B = u_sq / 1024. * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))
 
120
    delta_sigma = B * sin_sigma * ( 
 
121
        cos2_sigma_m + B / 4. * ( cos_sigma * ( -1 + 2 * cos2_sigma_m ** 2 ) -
 
122
            B / 6. * cos2_sigma_m * ( -3 + 4 * sin_sigma ** 2 ) *
 
123
            ( -3 + 4 * cos2_sigma_m ** 2 ) )
 
124
        )
 
125
    s = minor * A * (sigma - delta_sigma)
 
126
    
 
127
    return round(s,3) # round to 1mm precision
 
128
    
 
129
 
 
130
def calculate_dist_vincenty(latitudes, longitudes):
 
131
    """Returns the distance matrix from calculating dist_Vicenty
 
132
    
 
133
       latitudes, longitudes: list of values, have to be the same size
 
134
    """
 
135
    assert len(latitudes) == len(longitudes), "latitudes and longitudes must be lists of exactly the same size"
 
136
    
 
137
    size = len(latitudes)
 
138
    dtx_mtx = zeros([size, size])
 
139
        
 
140
    for i in range(size):
 
141
        for j in range(i,size):
 
142
            dtx_mtx[i,j] = dtx_mtx[j,i] = dist_vincenty(latitudes[i], longitudes[i], latitudes[j], longitudes[j])
 
143
        
 
144
    return dtx_mtx