calculate sphere distance

Topics: Developer Forum, Project Management Forum, User Forum
Sep 6, 2007 at 5:27 AM
test env:
srid = 50000, srtext = PROJCS"Mercator Spheric", GEOGCS["WGS84based_GCS", DATUM["WGS84based_Datum",SPHEROID["WGS84based_Sphere", 6378137, 0], TOWGS84[0, 0, 0, 0, 0, 0, 0]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]], AXIS["E", EAST], AXIS["N", NORTH]], PROJECTION["Mercator_1SP"], PARAMETER["False_Easting", 0], PARAMETER["False_Northing", 0], PARAMETER["Central_Meridian", 0], PARAMETER["Latitude_of_origin", 0], UNIT["metre", 1, AUTHORITY["EPSG", "9001"]], AXIS["East", EAST], AXIS["North", NORTH]

test step:
1. select ST.distance(ST.GeomFromText('POINT(-85.04344635735094 36.8243139564142)',50000),ST.GeomFromText('POINT(-85.078125 36.73888412439431)',50000))
it returns 0.0922001326188406
2. i calculate these two point using google, it is 10KM
3. i calculate these two points using postgis using distance function, it returns 0.0922001326188406, but when i use distance_sphere function, it returns 9.988KM.

Question: how can i get sphere distance in sql spatial extension????

Sep 6, 2007 at 7:28 AM
ST.distance returns the euclidian distance between the two points. That is:
distance^2 = (x2-x1)^2+(y2-y1)^2
This doesn't work very well with geographic coordinates, since they are in a spheric coordinate system.

Try to Google for "distance longitude latitude" and you will get plenty of ways for calculating the distance. Some are more accurate, others are very simple. You will need the radius of the earth, which depends on the datum you are using, in this case 6378.137 km. Others also use the flattening of the earth. The datum you specified here is a spheric datum that will give you inaccurate distances especially for points close to the poles.
Sep 6, 2007 at 8:22 AM
because i am not just calculating the distance between points, but also between point and polygon. What you said is fine for points i think. I am wondering if spatial extension support sphere distance? Is it related with the projection which i use? what is the unit of result i got?
Jun 10, 2009 at 10:09 PM

I'm new to the MsSqlSpatial tool and running into this same issue... so, does anyone have a good solution?

Thanks

Developer
Jun 10, 2009 at 11:06 PM

Hi akovar, MsSqlSpatial doesn't contain a spherical distance function you will need to create your own. see http://en.wikipedia.org/wiki/Great-circle_distance hth jd

Jun 11, 2009 at 9:50 PM

Thanks for the quick response johndiss!  A Great-circle distance is fine if I were only looking at points, but how about distances between higher dimension geometries?

I've been playing around with translating the Lat/Long coordinate into spatial cartesian coordinates (XYZ), finding the straight line distance between the two, and then correcting for arc.  This seems to be rather accurate for the locations I've tried (the largest error I've seen has been about 6 meters), though for short distances this is a much larger percentage error.  Also, for short distances, the correction for arc isn't as vital, thus the closest distance between two non-point geometries should be pretty accurate as well.  The equations I've used are as follows (assuming lat/lng has already been translated to radians):

X = cos(lng)*sin(90 - lat)

Y = sin(lng)*sin(90 - lat)

Z = cos(90 - lat)

dist = 2*R*asin(ST.Distance(xyzPoint1, xyzPoint2) / 2)

 

Again, the errors I'm seeing are about 3 meters, which (I believe) is well within the error caused by assuming the earth is a sphere.  Am I missing anything major with this process?

Thanks, again

Developer
Jun 12, 2009 at 8:42 AM

Hi again Akovar, I am no expert on such things but you may like to investigate the ST.Transform function so that you can get the points transformed onto a local grid - this may improve your accuracy further so something like:

DECLARE @targetSrid int
SET @targetSrid  = xxxxx

DECLARE @lonLatGeom1 VarBinary(8000)
DECLARE @lonLatGeom2 VarBinary(8000)

SET @lonLatGeom1.....

SELECT ST.Distance(
ST.Transform(@lonlatGeom1, @targetSrid), 
ST.Transform(@lonLatGeom2, @targetSrid)
) as Distance
(this will require that the input geometries have an appropriate SRID, if they dont you can use ST.SetSRID())
hth jd
Jun 12, 2009 at 8:12 PM

It seems I made a rather big error with my sphere logic... MsSqlSpatial doesn't handle 3D geometry!  For some reason the first tests I ran kept a constant latitude, so it didn't really matter that the Z value of the point was dropped, but with any significant changes in latitude, the process fails horribly!

I think the Transform function is probably the best bet, but now the question becomes how to select the correct targetSRID.  And then correctly handling conditions near the borders between different SRID change overs... if that makes any sense...

Thanks again for the help jd, now I just have to make it work!

Developer
Jun 15, 2009 at 9:44 AM

Hi akovar, you may find http://www.spatialreference.org/ useful. Wrt 'changeover' regions I am not sure how you would need to deal with them - MsSqlSpatial will insist that both geometries are in the same reference system.