# Source code for libpysal.cg.sphere

"""
sphere: Tools for working with spherical geometry.

Author(s):
Charles R Schmidt schmidtc@gmail.com
Luc Anselin luc.anselin@asu.edu
Xun Li xun.li@asu.edu

"""

__author__ = "Charles R Schmidt <schmidtc@gmail.com>, Luc Anselin <luc.anselin@asu.edu, Xun Li <xun.li@asu.edu"

import math
import numpy
import scipy.spatial
import scipy.constants
from scipy.spatial.distance import euclidean
from math import pi, cos, sin

__all__ = [
'brute_knn', 'fast_knn', 'fast_threshold', 'linear2arcdist', 'toLngLat',
'toXYZ', 'lonlat', 'harcdist', 'geointerpolate', 'geogrid'
]

scipy.constants.kilo) / scipy.constants.mile

"""
Arc distance between two points on a sphere.

Parameters
----------
pt0 : point
assumed to be in form (lng,lat)
pt1 : point
assumed to be in form (lng,lat)

Source: http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html

Returns
-------
The arc distance between pt0 and pt1 using supplied radius

Examples
--------
>>> pt0 = (0,0)
>>> pt1 = (180,0)
True
"""

"""
Convert an arc distance (spherical earth) to a linear distance (R3) in the unit sphere.

Examples
--------
>>> pt0 = (0,0)
>>> pt1 = (180,0)
True
2.0
"""
c = 2 * math.pi * radius
d = (2 - (2 * math.cos(math.radians((arc_dist * 360.0) / c))))**(0.5)
return d

"""
Convert a linear distance in the unit sphere (R3) to an arc distance based on supplied radius

Examples
--------
>>> pt0 = (0,0)
>>> pt1 = (180,0)
True
"""
if linear_dist == float('inf'):
return float('inf')
elif linear_dist > 2.0:
raise ValueError(
"linear_dist, must not exceed the diameter of the unit sphere, 2.0"
)
c = 2 * math.pi * radius
a2 = linear_dist**2
theta = math.degrees(math.acos((2 - a2) / (2.)))
d = (theta * c) / 360.0
return d

[docs]def toXYZ(pt):
"""
Convert a point's latitude and longitude to x,y,z

Parameters
----------
pt0 : point
assumed to be in form (lng,lat)
pt1 : point
assumed to be in form (lng,lat)

Returns
-------
x, y, z
"""
phi, theta = phi + pi, theta + (pi / 2)
x = 1 * sin(theta) * cos(phi)
y = 1 * sin(theta) * sin(phi)
z = 1 * cos(theta)
return x, y, z

[docs]def toLngLat(xyz):
"""
Convert x,y,z to latitude and longitude

"""
x, y, z = xyz
if z == -1 or z == 1:
phi = 0
else:
phi = math.atan2(y, x)
if phi > 0:
phi = phi - math.pi
elif phi < 0:
phi = phi + math.pi
theta = math.acos(z) - (math.pi / 2)
return phi, theta

[docs]def brute_knn(pts, k, mode='arc'):
"""
valid modes are ['arc','xrz']
"""
n = len(pts)
full = numpy.zeros((n, n))
for i in range(n):
for j in range(i + 1, n):
if mode == 'arc':
lng0, lat0 = pts[i]
lng1, lat1 = pts[j]
elif mode == 'xyz':
dist = euclidean(pts[i], pts[j])
full[i, j] = dist
full[j, i] = dist
w = {}
for i in range(n):
w[i] = full[i].argsort()[1:k + 1].tolist()
return w

[docs]def fast_knn(pts, k, return_dist=False):
"""
Computes k nearest neighbors on a sphere.

Parameters
----------
pts :  list of x,y pairs
k   :  int
Number of points to query
return_dist : bool
Return distances in the 'wd' container object

Returns
-------
wn  :  list
list of neighbors
wd  : list
list of neighbor distances (optional)

"""
pts = numpy.array(pts)
kd = scipy.spatial.KDTree(pts)
d, w = kd.query(pts, k + 1)
w = w[:, 1:]
wn = {}
for i in range(len(pts)):
wn[i] = w[i].tolist()
if return_dist:
d = d[:, 1:]
wd = {}
for i in range(len(pts)):
wd[i] = [
for x in d[i].tolist()
]
return wn, wd
return wn

"""
Find all neighbors on a sphere within a threshold distance

Parameters
----------

pointslist : list of lat-lon tuples (Note, has to be a list, even for one point)

dist: float
threshold distance

Returns
-------

dict:
key is id of point, value is a list of ids for other points within dist of key point

"""
kd = scipy.spatial.KDTree(pts)
r = kd.query_ball_tree(kd, d)
wd = {}
for i in range(len(pts)):
l = r[i]
l.remove(i)
wd[i] = l
return wd

########### new functions

[docs]def lonlat(pointslist):
"""
Converts point order from lat-lon tuples to lon-lat (x,y) tuples

Parameters
----------

pointslist : list of lat-lon tuples (Note, has to be a list, even for one point)

Returns
-------

newpts      : list with tuples of points in lon-lat order

Examples
--------
>>> points = [(41.981417, -87.893517), (41.980396, -87.776787), (41.980906, -87.696450)]
>>> newpoints = lonlat(points)
>>> newpoints
[(-87.893517, 41.981417), (-87.776787, 41.980396), (-87.69645, 41.980906)]

"""
newpts = [(i[1], i[0]) for i in pointslist]
return newpts

def haversine(x):
"""
Computes the haversine formula

Parameters
----------

Returns
-------
: square of sine of half the radian (the haversine formula)

Examples
--------
>>> haversine(math.pi)     # is 180 in radians, hence sin of 90 = 1
1.0

"""

x = math.sin(x / 2)
return x * x

# Lambda functions

d2r = lambda x: x * math.pi / 180.0

r2d = lambda x: x * 180.0 / math.pi

"""
Radian angle between two points on a sphere in lon-lat (x,y)

Parameters
----------
p0    : first point as a lon,lat tuple
p1    : second point as a lon,lat tuple

Returns
-------

Examples
--------
>>> p0 = (-87.893517, 41.981417)
>>> p1 = (-87.519295, 41.657498)
0.007460167953189258

Notes
-----
Uses haversine formula, function haversine and degree to radian
conversion lambda function d2r

"""
x0, y0 = d2r(p0[0]), d2r(p0[1])
x1, y1 = d2r(p1[0]), d2r(p1[1])
d = 2.0 * math.asin(
math.sqrt(
haversine(y1 - y0) +
math.cos(y0) * math.cos(y1) * haversine(x1 - x0)))
return d

"""
Alternative arc distance function, uses haversine formula

Parameters
----------
p0       : first point as a tuple in decimal degrees
p1       : second point as a tuple in decimal degrees
lonx     : boolean to assess the order of the coordinates,
for lon,lat (default) = True, for lat,lon = False
radius   : radius of the earth at the equator as a sphere

Returns
-------
d        : distance in units specified, km, miles or radians (for None)

Examples
--------
>>> p0 = (-87.893517, 41.981417)
>>> p1 = (-87.519295, 41.657498)
>>> harcdist(p0,p1)
47.52873002976876
0.007460167953189258

Notes
-----

"""
if not (lonx):
p = lonlat([p0, p1])
p0 = p[0]
p1 = p[1]

return d

[docs]def geointerpolate(p0, p1, t, lonx=True):
"""
Finds a point on a sphere along the great circle distance between two points
on a sphere
also known as a way point in great circle navigation

Parameters
----------
p0       : first point as a tuple in decimal degrees
p1       : second point as a tuple in decimal degrees
t        : proportion along great circle distance between p0 and p1
e.g., t=0.5 would find the mid-point
lonx     : boolean to assess the order of the coordinates,
for lon,lat (default) = True, for lat,lon = False

Returns
-------
x,y      : tuple in decimal degrees of lon-lat (default) or lat-lon,
depending on setting of lonx; in other words, the same
order is used as for the input

Examples
--------
>>> p0 = (-87.893517, 41.981417)
>>> p1 = (-87.519295, 41.657498)
>>> geointerpolate(p0,p1,0.1)          # using lon-lat
(-87.85592403438788, 41.949079912574796)
>>> p3 = (41.981417, -87.893517)
>>> p4 = (41.657498, -87.519295)
>>> geointerpolate(p3,p4,0.1,lonx=False)   # using lat-lon
(41.949079912574796, -87.85592403438788)

"""

if not (lonx):
p = lonlat([p0, p1])
p0 = p[0]
p1 = p[1]

k = 1.0 / math.sin(d)
t = t * d
A = math.sin(d - t) * k
B = math.sin(t) * k

x0, y0 = d2r(p0[0]), d2r(p0[1])
x1, y1 = d2r(p1[0]), d2r(p1[1])

x = A * math.cos(y0) * math.cos(x0) + B * math.cos(y1) * math.cos(x1)
y = A * math.cos(y0) * math.sin(x0) + B * math.cos(y1) * math.sin(x1)
z = A * math.sin(y0) + B * math.sin(y1)

newpx = r2d(math.atan2(y, x))
newpy = r2d(math.atan2(z, math.sqrt(x * x + y * y)))
if not (lonx):
return newpy, newpx
return newpx, newpy

[docs]def geogrid(pup, pdown, k, lonx=True):
"""
Computes a k+1 by k+1 set of grid points for a bounding box in lat-lon
uses geointerpolate

Parameters
----------
pup     : tuple with lat-lon or lon-lat for upper left corner of bounding box
pdown   : tuple with lat-lon or lon-lat for lower right corner of bounding box
k       : number of grid cells (grid points will be one more)
lonx    : boolean to assess the order of the coordinates,
for lon,lat (default) = True, for lat,lon = False

Returns
-------
grid    : list of tuples with lat-lon or lon-lat for grid points, row by row,
starting with the top row and moving to the bottom; coordinate tuples
are returned in same order as input

Examples
--------
>>> pup = (42.023768,-87.946389)    # Arlington Heights IL
>>> pdown = (41.644415,-87.524102)  # Hammond, IN
>>> geogrid(pup,pdown,3,lonx=False)
[(42.023768, -87.946389), (42.02393997819538, -87.80562679358316), (42.02393997819538, -87.66486420641684), (42.023768, -87.524102), (41.897317, -87.94638900000001), (41.8974888973743, -87.80562679296166), (41.8974888973743, -87.66486420703835), (41.897317, -87.524102), (41.770866000000005, -87.94638900000001), (41.77103781320412, -87.80562679234043), (41.77103781320412, -87.66486420765956), (41.770866000000005, -87.524102), (41.644415, -87.946389), (41.64458672568646, -87.80562679171955), (41.64458672568646, -87.66486420828045), (41.644415, -87.524102)]

"""
if lonx:
corners = [pup, pdown]
else:
corners = lonlat([pup, pdown])
tpoints = [float(i) / k for i in range(k)[1:]]
leftcorners = [corners[0], (corners[0][0], corners[1][1])]
rightcorners = [(corners[1][0], corners[0][1]), corners[1]]
leftside = [leftcorners[0]]
rightside = [rightcorners[0]]
for t in tpoints:
newpl = geointerpolate(leftcorners[0], leftcorners[1], t)
leftside.append(newpl)
newpr = geointerpolate(rightcorners[0], rightcorners[1], t)
rightside.append(newpr)
leftside.append(leftcorners[1])
rightside.append(rightcorners[1])

grid = []
for i in range(len(leftside)):
grid.append(leftside[i])
for t in tpoints:
newp = geointerpolate(leftside[i], rightside[i], t)
grid.append(newp)
grid.append(rightside[i])
if not (lonx):
grid = lonlat(grid)
return grid