Question

New to Pysal and just getting to grips with loading, reading etc. I have some test point data, in Decimal Degrees (WGS84) and want to use the get_point_dist function. After running it, I get the correct response but would like the distance in meters/KM etc rather the DD.

Any ideas how to do this, or is the answer, "project the input data"?

Was it helpful?

Solution

The get_points_dist method in PySAL operates only in euclidean space, so you would need to project your input data to use it.

However you could use the arcdist method.

Note: Point tuples should be provided in (Long, Lat) order.

For example:

>>> import pysal
>>> pt0 = (-74.0, 40.7)
>>> pt1 = (-77.0, 38.8)
>>> pysal.cg.arcdist(pt0, pt1)

arcdist does have some limitations, it assumes a perfect sphere and will not take into account the geoid (WGS84). It defaults to KiloMeters, to change this you must specify the radius of the sphere (earth) in whatever units you need. Constants for KM and Miles are provides, pysal.cg.sphere.RADIUS_EARTH_KM pysal.cg.sphere.RADIUS_EARTH_MILES

Example:

>>> pysal.cg.arcdist(pt0, pt1, pysal.cg.sphere.RADIUS_EARTH_KM)

The KDTree in PySAL will also work with spherical coordinates and will provide accurate nearest neighbor results, even at the anit-meridian (dateline/180 lng/-180 lng) and the polls.

>>> pts = [(-179.0,0.0), (170,0.0), (179.0,0.0)]
>>> kd = pysal.cg.KDTree(pts, distance_metric='Arc', radius = pysal.cg.sphere.RADIUS_EARTH_KM)
>>> d, i = kd.query((-180.0,0.0), k=3)
>>> d
array([  111.19492664,   111.19492664,  1111.94926645])
>>> i
array([0, 2, 1])
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top