Haversine/Pythagoras Distance (again)

Please start any new threads on our new site at https://forums.sqlteam.com. We've got lots of great SQL Server experts to answer whatever question you can come up with.

All Forums

SQL Server 2005 Forums

Transact-SQL (2005)

Haversine/Pythagoras Distance (again)

Author

Topic

Kristen
Test

22859 Posts

Posted - 2007-10-29 : 16:55:30

I've been looking at the "cheapest" CPU cost of calculating distance between two points.

Generally I am asked to answer questions along the lines of:

1) Order by distance
2) Show distance (in miles / kilometres)
3) Find locations within X miles / kilometres

My inclination is to store the Lat/Long on each "location" record - in either degrees or radians. I think that radians saves one function call in the subsequent calculations.

And from that to use Pythagoras to provide a "rough" calculation, as it needs to work at my Latitudes; but based on my sample data the accuracy seems to be way off - although its quite possible that I've mucked up the calculations.

My test data is based around 1 mile, 10 mile and 50 mile radius. (These were taken from a map, so hopefully I have measured the Lat/Long correctly, but it is worth checking that my data is sensible)


DECLARE	@TestData TABLE
(
	ID	int identity(1,1),
	LatD	float, 
	LongD	float, 
	Place	varchar(20),
	LatR	float, 
	LongR	float,
	HaversineTemp		float,
	HaversineDistance	float,
	PythagorasDistance	float
)

INSERT INTO @TestData(LatD, LongD, Place)
-- Lat, Long, Place
SELECT [LatD] = 52.40858, [LongD] = -1.50361, [Place]='Coventry' UNION ALL
SELECT 52.42315, -1.50361, '1 mile North' UNION ALL
SELECT 52.39426, -1.50361, '1 mile South' UNION ALL
SELECT 52.40858, -1.47996, '1 mile East' UNION ALL
SELECT 52.40858, -1.52732, '1 mile West' UNION ALL

SELECT 52.42315, -1.47996, '1 mile NorthEast' UNION ALL
SELECT 52.39426, -1.47996, '1 mile SouthEast' UNION ALL
SELECT 52.42315, -1.52732, '1 mile NorthWest' UNION ALL
SELECT 52.39426, -1.52732, '1 mile SouthWest' UNION ALL

SELECT 52.55094, -1.50361, '10 mile North' UNION ALL
SELECT 52.26536, -1.50361, '10 mile South' UNION ALL
SELECT 52.40858, -1.26690, '10 mile East' UNION ALL
SELECT 52.40858, -1.74169, '10 mile West' UNION ALL

SELECT 52.55094, -1.26690, '10 mile NorthEast' UNION ALL
SELECT 52.26536, -1.26690, '10 mile SouthEast' UNION ALL
SELECT 52.55094, -1.74169, '10 mile NorthWest' UNION ALL
SELECT 52.26536, -1.74169, '10 mile SouthWest' UNION ALL

SELECT 53.13351, -1.50361, '50 mile North' UNION ALL
SELECT 51.68883, -1.50361, '50 mile South' UNION ALL
SELECT 52.40858, -0.32553, '50 mile East' UNION ALL
SELECT 52.40858, -2.69110, '50 mile West' UNION ALL

SELECT 53.13351, -0.32553, '50 mile NorthEast' UNION ALL
SELECT 51.68883, -0.32553, '50 mile SouthEast' UNION ALL
SELECT 53.13351, -2.69110, '50 mile NorthWest' UNION ALL
SELECT 51.68883, -2.69110, '50 mile SouthWest' UNION ALL

SELECT 58.64068, -3.07649, 'John o Groats' UNION ALL
SELECT 59.35703, -3.07649, 'JoG 50 mile North' UNION ALL
SELECT 57.92433, -3.07649, 'JoG 50 mile South' UNION ALL
SELECT 58.64068, -1.69155, 'JoG 50 mile East' UNION ALL
SELECT 58.64068, -4.46144, 'JoG 50 mile West' UNION ALL

SELECT 49.95374, -5.21061, 'Lizard' UNION ALL
SELECT 50.67168, -5.21061, 'Liz 50 mile North' UNION ALL
SELECT 49.23580, -5.21061, 'Liz 50 mile South' UNION ALL
SELECT 49.95374, -4.08781, 'Liz 50 mile East' UNION ALL
SELECT 49.95374, -6.33341, 'Liz 50 mile West' 

UPDATE	U
SET	LatR = radians(LatD),
	LongR = radians(LongD)
FROM	@TestData AS U

SELECT	*
FROM	@TestData
ORDER BY ID

DECLARE	@MyLatR float,
	@MyLongR float

-- Get "Home" location
SELECT	@MyLatR = [LatR],
	@MyLongR = [LongR]
FROM	@TestData
WHERE	Place = 'Coventry'

-- Distance
UPDATE U
SET	HaversineTemp = sqrt(square(sin((LatR - @MyLatR)/2.0E)) + 
	(cos(@MyLatR) * cos(LatR) * square(sin((LongR - @MyLongR)/2.0E))) )
FROM	@TestData AS U

UPDATE U
SET	HaversineDistance = 6371.0E * ( 2.0E * asin(case when 1.0E < HaversineTemp then 1.0E else HaversineTemp end ))
FROM	@TestData AS U


UPDATE U
SET	PythagorasDistance = 6371.0E * SQRT(POWER(LatR - @MyLatR, 2) + Power(LongR - @MyLongR, 2))
FROM	@TestData AS U

SELECT	*
FROM	@TestData
ORDER BY ID

EDIT: Fixed some typos in the data
EDIT: Data added for John O'Groats (North) and Lizard (South)
Kristen

Kristen
Test

22859 Posts

Posted - 2007-10-30 : 10:25:55

<bump>

tm
Posting Yak Master

160 Posts

Posted - 2007-10-30 : 15:08:53

quote:
Originally posted by Kristen

<bump>

What does <bump> mean?

I thought I would have a look at what your doing and it's giving me a headache.

I checked out .. http://en.wikipedia.org/wiki/Great-circle_distance
It's all greek to me

Kristen
Test

22859 Posts

Posted - 2007-10-30 : 16:50:09

"What does <bump> mean?"

Just giving the thread a fresh date in the hope that folk could come up with some advice.

I can't understand why the variation is so great between the two methods because I thought it was only fractions of a percent, so I reckon I have done something fundamentally wrong - can't spot what though

Kristen

Michael Valentine Jones
Yak DBA Kernel (pronounced Colonel)

7020 Posts

Posted - 2007-10-30 : 17:59:50

It looks to me like your Haversine distance calculation is correct. It produces the same result as the function here:
http://www.sqlteam.com/forums/topic.asp?TOPIC_ID=81360

That means the problem is in the PythagorasDistance, but I guess you already figured that out.

You have use a distance correction factor for the latitude difference. One degree (or radian) is less distance East/West at a higher latitude than it is at the Equator.

When I ran the following update on your test data right after the insert, the results were very close. The update sets the latitude very close to the Equator, so there is less error.

update @TestData set LatD = LatD - 52.00E

For your test data, this correction works well:

UPDATE U
SET	PythagorasDistance = 
	(6371.0E*0.6101451639012192E) * 
	SQRT(POWER(LatR - @MyLatR, 2) + Power(LongR - @MyLongR, 2))
FROM	@TestData AS U

I got the correction factor this way,


select
	[Latitude Correction Factor] =
	-- Distance of .0001 degree at 52.4 Degrees latitude
	dbo.F_GREAT_CIRCLE_DISTANCE( 52.4, 0, 52.4, 0.0001 )/
	-- Distance of .0001 degree at equator latitude
	dbo.F_GREAT_CIRCLE_DISTANCE( 0, 0, 0, 0.0001)

CODO ERGO SUM

Kristen
Test

22859 Posts

Posted - 2007-10-30 : 18:18:40

You are a star MVJ, many thanks for that. I would have been scratching my head all the way to the North Pole!

Kristen

spirit1
Cybernetic Yak Master

11752 Posts

Posted - 2007-10-30 : 18:38:02

things will be much easier with geo coding datatype in Katmai... finnaly...

_______________________________________________
Causing trouble since 1980
blog: http://weblogs.sqlteam.com/mladenp
SSMS Add-in that does a few things: www.ssmstoolspack.com

Kristen
Test

22859 Posts

Posted - 2007-10-30 : 18:40:24

Hmmm ... so distances Due North/South short significant discrepancy between the two methods.

After that experiement I'm not sure that Pythagoras is going to be any use to me ... on a 50 mile journey due North/South its going to be 20KM/12Miles out - that's quite a bit ...

Unless I've misinterpreted the figures?

I fixed some typos in the data above (which I had fixed locally before making this discovery)

Kristen

Michael Valentine Jones
Yak DBA Kernel (pronounced Colonel)

7020 Posts

Posted - 2007-10-30 : 21:34:09

If you just want to get the distance, you might as well use the Haversine calculation. I doubt that it costs that much more to calculate, and you don't have to worry much about accuracy.

What is the nature of your application? Are you doing searches to find locations within a certain distance of a location? Or just calculating the distance between two points?

CODO ERGO SUM

SwePeso
Patron Saint of Lost Yaks

30421 Posts

Posted - 2007-10-31 : 15:55:09

What if you for pythagoras do not convert to radians beforehand?
Just apply theorem to coordinates and apply factor multiplication afterwards?

E 12°55'05.25"
N 56°04'39.16"

Kristen
Test

22859 Posts

Posted - 2007-11-01 : 09:03:53

"What is the nature of your application? Are you doing searches to find locations within a certain distance of a location? Or just calculating the distance between two points?"

1) Order by distance
2) Show distance (in miles / kilometres)
3) Find locations within X miles / kilometres

Generally its "order by distance" and then "show actual distance". This would be on data that is paged, so I was looking for the cheapest CPU cost route for the Order By (at the expense of some accuracy), and then doing an accurate calculation to show the Distance on a page of data.

This is a very high volume query, selecting from a large number of rows, and producing larg-ish resultsets (i.e. several pages). The user is going to only look at the first couple of pages, and "distance" is going to be the most commonly used sort order for paging.

I have Lat/Long for the User and each of the Results. I can store Lat/Long as Radians, Degrees, or anything else that would be helpful to this query

I've changed MVJ's suggestion to:


UPDATE U
SET	PythagorasDistance = 
	(6371.0E *0.6101451639012192E) * 
	SQRT(POWER(LatR - @MyLatR, 2) + Power((LongR - @MyLongR) * 0.6101451639012192E, 2))
FROM	@TestData AS U

and that reduces the East/West error (for 50 mile radius) from 50KM to 0.5KM.

Is this a valid adjustment?

Should I "tune" the size of the adjustment according to the latitude value of @MyLatR to make it more generic? (that's obviously a once-per-query Great Circle Distance calculation, rather than a once-per-row

Thanks,

Kristen

Kristen
Test

22859 Posts

Posted - 2007-11-01 : 09:58:00

I've added some data for the North and South of the UK.

Using the formula above for calculating the adjustment, and actually minimising the adjust by trial-and-error I get:


-- Location              Calculated            Trial-and-Error
-- ====================  ==================    ===============
-- John o'Groats (North) 0.5204034795563881	0.5639
-- Coventry = (Middle)   0.6101451639012192
-- Lizard = (South)      0.64338371678284456	0.6265

so less adjustment further North / South than Great Circle would suggest. Might just be chance that it worked OK for Coventry for starter?

Kristen

Subscribe to SQLTeam.com

SQLTeam.com Articles via RSS

SQLTeam.com Weblog via RSS

- Advertisement -

Resources