11

I am working with an ESRI stack, storing my layers in a sql-spatial-enabled-SDE- geodatabase (Geometry type, Web Mercator-3857).

I am building a web mapping application, so by default, the tiles are also in web mercator, 3857.

Via stored procs, i use STDistance to query distances from a user's location (Coordinates also in web mercator) to the various layers.

Problem is, due to the distortion of web mercator, my distance calcs are increasingly off, the further from the equator they are made.

I've thought of storing my layers in sql-spatial-geography (rather than geometry) type, but:

  • I imagine my distance queries will take much longer (distance calcs on spherical surface)
  • i will need to reimport a lot of data
  • arcgis services won't be as fast as they will need to project on the fly

If I go to Google maps, and do a distance calc, the distance returned is much more accurate, even in Nortern/southern regions, so I assume Google must be correcting for distortion caused by the web mercator projection.

My question then: Is there a simple factor value that can be applied to distance calcs done in web mercator projection to get the "correct" distance?

5 Answers 5

15

For short distances, you could multiply the calculated distance with cos(lat), since the scale of the mercator projection is proportional to the secant of the latitude (secant is 1/cos). Also see http://en.wikipedia.org/wiki/Mercator_projection#Mathematics_of_the_projection


Addendum thanks to @jeremiah-england: while the above correction would be correct when it comes to true Mercator projections, the Web Mercator (EPSG:3857) is not a Mercator. EPSG calls it a "Pseudo-Mercator." The problem is that it uses WGS84, an eliptical model, and projects it using spherical mercator calculations (which Google used because they are faster). If you scale your distances by 1/cos(phi) with the web Mercator, you'll be about 0.6% off at the equator. See Noel Zinn's presentation on Web Mercator for more details.

According to the above-mentioned presentation, the following method could be used to calculate more accurate distances from web mercator coordinates. Given dx — horizontal coordinate difference (W-E direction), and dy — vertical coordinate difference (S-N direction):

e = 0.081819191
adjustedX = dx * cos(lat) / sqrt(1 - e^2 * sin(lat)^2)
adjustedY = dy * cos(lat) * (1 - e^2) / pow(1 - e^2 * sin(lat)^2, 3/2)
adjustedDistance = hypot(adjustedX, adjustedY)

The ratio between this adjustment and cos(lat) is larger in the S-N direction, ranging from 0.9933 at the equator to 1.0034 at the poles. The W-E direction ratio starts with 1 at the equator and grows to 1.0034 at the poles.

Note that this correction still works reasonably well only for short distances, where planar geometry of the Earth's surface can be assumed.

9
  • 1
    And for longer distances, you can use the mean latitude of the two endpoints: cos((l1 + l2)/2) which will give you the rhumb-line/constant course distance, rather than a great-circle distance. Commented Sep 12, 2011 at 13:08
  • that only changes the spheroid, but does not gives the same precision like a local projection.
    – falcacibar
    Commented Sep 12, 2011 at 14:23
  • 1
    @falcibar - i don't see how choosing the mean latitude changes the spheroid
    – mkadunc
    Commented Sep 12, 2011 at 14:36
  • @MerseyViking: thanks, forgot to mention that the best latitude to use for the calculation would be the mean of the two compared points.
    – mkadunc
    Commented Sep 12, 2011 at 14:38
  • 1
    +1 for the answer. @Mersey: Why does the mean latitude correction work? After all, distortions in the Mercator can become arbitrarily large and departures of the loxodromes from the geodesics can also be large. It seems that there are some potentially enormous errors that could be made for which a simple correction would be unreliable.
    – whuber
    Commented Sep 12, 2011 at 15:31
6

I'd consider your second option of storing your data in the GEOGRAPHY format again if you are looking for accurate results on global data.

There is nothing stopping you from having two spatial fields in a table - one in Mercator as a GEOMETRY type and one in WGS84 as a GEOGRAPHY type (at least not in SQL Server, I am not sure about ArcSDE).

You should be able to create a simple geoprocessing script that will populate both fields with your original data. If there is editing and frequent updates going on then this may not be an option.

Once you have the two fields you can carry on using your Mercator for fast display, and convert user entered points to Lat/Lon for distances.

This has two major advantages:

  • you can also get more accurate areas and lengths of features based on user queries
  • you don't need to worry about forgetting to add custom code for each query which deals with measurements

The query speeds will be more complex, but it may be unnoticeable to a user. You may want to test with one feature class before deciding on a solution. You'll also have to convert the user entered points from Mercator (which should be trivial, and could be done in the browser).

Even with the geography type though there is still an error margin:

The error tolerance for the geography methods can be as large as 1.0e-7 * extents

From MSDN

2
  • Unfortunately SDE will only allow for one spatial column: help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… - I haven't tried to create a view and register it with SDE, but that might be a possible approach. Commented Sep 12, 2011 at 17:00
  • @Allan - good to know. Another positive for using SDE then! A table with just the 4326 geometry and a join back to the original table + spatial view should achieve the same purpose Commented Sep 15, 2011 at 14:23
3

An alternative suggestion from ESRI is to use a "Geometry service", which involves sending the geometries to an ArcGIS Server and having the result returned. If you don't expect any bottleneck issues by using a web service, this can be a very effective approach. I have used the same approach in a Silverlight application.

Here is an original blog entry from ESRI: http://blogs.esri.com/Dev/blogs/arcgisserver/archive/2010/03/05/Measuring-distances-and-areas-when-your-map-uses-the-Mercator-projection.aspx

In the blog there is a simplified JavaScript example, and there is also a full example application here: http://serverapps.esri.com/javascript_examples/compare_measurements.htm

0

As google earth does, you could transform the geometry to local WGS84 zone projection, or a enhanced WGS84 projection like SIRGAS in South America.

You could look in http://www.spatialreference.org for the coordinates of almost all "zonified" projection, and you could make a table, etc, to transform them depends on the zone.

We imagine a table zones

|   minx     |    miny    |    maxx    |   maxy     |  srid  |
+------------+------------+------------+------------+--------+
|234567.34314|234567.34334|234567.34334|234567.34334|  1234  |

all minx, miny, maxx, maxy values stored in Spherical Mercator (Web Mercator). so i will use postgis as example.

SELECT ST_Distance(
         ST_Transform( -- transform/reproject
            ST_SetSRID(geom_line, 3857) -- geom_line with forced srid assignation to web mercator
            , ( -- here we get the first SRID from the spatial position of geom_line
                 SELECT  srid
                 FROM    zones
                 WHERE   ST_Contains(
                           ST_MakeBox2d(ST_Point(minx,miny),ST_Point(maxx,maxy))
                           , geom_line
                         )
                 LIMIT 1
            ))

Hope this be useful, have a nice day.

3
  • I can't be certain, but based on Blomster's use of "STDistance" in his question it sounds like he isn't using PostGIS. If Blomster is using SQL Server there is no ST_Transform functionality. Commented Sep 12, 2011 at 14:21
  • I give the process and the best way that i could solve it, he could deal with the rest, anway he could use a software reprojection, but is really a pity that a SQL Server doesn't deal with projections.
    – falcacibar
    Commented Sep 12, 2011 at 14:42
  • maybe CLR enabled and the .net library nettopologysuite could helps?
    – falcacibar
    Commented Sep 12, 2011 at 18:33
0

OpenLayers have a utility method to calculate the distance between two EPSG:4326 (WGS84) points on an ellipsoid. Since OpenLayers is open source, you can see how it's implemented here: https://github.com/openlayers/openlayers/blob/release-2.12/lib/OpenLayers/Util.js#L750

For those using OpenLayers and the Measure control, just enable geodesic in the options litteral to use this calculation.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.