Implement Web Mercator (#67)

Author: c42f

README.md

@@ -419,6 +419,18 @@ maximal efficiency when performing multiple `transform`s. The transformation can
 be inverted with `inv` to perform the reverse transformation with respect to the
 same origin.
 
+#### Web Mercator support
+
+We support the Web Mercator / Pseudo Mercator projection with the
+`WebMercatorfromLLA` and `LLAfromWebMercator` transformations for
+interoperability with many web mapping systems. The scaling of the
+northing and easting is defined to be meters at the Equator, the same as how
+proj handles this (see https://proj.org/operations/projections/webmerc.html ).
+
+If you need to deal with web mapping tile coordinate systems (zoom levels and
+pixel coordinates, etc) these could be added by composing another
+transformation on top of the web mercator projection defined in this package.
+
 #### Composed transformations
 
 Many other methods are defined as convenience constructors for composed

src/Geodesy.jl

@@ -42,6 +42,7 @@ export
     UTMfromLLA, LLAfromUTM, UTMfromECEF, ECEFfromUTM, ENUfromUTM, UTMfromENU,
     UTMZfromLLA, LLAfromUTMZ, UTMZfromECEF, ECEFfromUTMZ, ENUfromUTMZ, UTMZfromENU,
     UTMZfromUTM, UTMfromUTMZ,
+    WebMercatorfromLLA, LLAfromWebMercator,
 
     # Datum transformations
     datum_shift_ECEF,
@@ -55,6 +56,7 @@ include("datums.jl")
 include("points.jl")
 include("transverse_mercator.jl")
 include("polar_stereographic.jl")
+include("web_mercator.jl")
 include("transformations.jl")
 include("conversion.jl")
 include("distances.jl")

src/web_mercator.jl

@@ -0,0 +1,78 @@
+"""
+    WebMercatorfromLLA(datum=wgs84)
+
+Convert from LLA to Web Mercator / Pseudo Mercator, following the convention of
+proj, which uses a scaling factor of the semi-major axis of the ellipsoid
+https://proj.org/operations/projections/webmerc.html.
+
+!!! warning
+    For web mapping applications this projection is ubiquitous due to its
+    simplicity of implementation, but this simplicity gives rise to poor
+    mathematical properties: it doesn't preserve angles away from the equator.
+    Other projections should be preferred if possible, and especially when
+    measurement is important. See
+    https://earth-info.nga.mil/GandG/wgs84/web_mercator/(U)%20NGA_SIG_0011_1.0.0_WEBMERC.pdf
+    for an extended discussion.
+"""
+struct WebMercatorfromLLA <: Transformation
+    el::Ellipsoid
+end
+
+WebMercatorfromLLA(d::Datum=wgs84) = WebMercatorfromLLA(ellipsoid(d))
+
+"""
+    LLAfromWebMercator
+
+Inverse of WebMercatorfromLLA — see the docs for that transformation.
+"""
+struct LLAfromWebMercator <: Transformation
+    el::Ellipsoid
+end
+
+LLAfromWebMercator(d::Datum=wgs84) = LLAfromWebMercator(ellipsoid(d))
+
+
+function (trans::WebMercatorfromLLA)(lla::LLA)
+    xy = trans(LatLon(lla.lat, lla.lon))
+    SA[xy[1], xy[2], lla.alt]
+end
+
+function (trans::WebMercatorfromLLA)(ll::LatLon)
+    lat_lim = 85.06 # according to https://epsg.io/3857
+    if abs(ll.lat) > lat_lim
+        throw(ArgumentError("Exceeded Web Mercator latitude bounds"))
+    end
+    x, y = web_mercator_forward(ll.lat, ll.lon, trans.el.a)
+    SA[x, y]
+end
+
+function (trans::LLAfromWebMercator)(point::AbstractVector)
+    x = point[1]
+    y = point[2]
+    lat, lon = web_mercator_reverse(x, y, trans.el.a)
+    if length(point) == 2
+        LatLon(lat, lon)
+    elseif length(point) == 3
+        LLA(lat, lon, point[3])
+    else
+        throw(ArgumentError("Expected input vector of length 2 or 3"))
+    end
+end
+
+Base.inv(trans::WebMercatorfromLLA) = LLAfromWebMercator(trans.el)
+Base.inv(trans::LLAfromWebMercator) = WebMercatorfromLLA(trans.el)
+
+# Web / Psuedo Mercator. Following the scaling convention of proj, the scaling
+# will be set to the ellipsoid semi-major axis.
+#   https://proj.org/operations/projections/webmerc.html
+function web_mercator_forward(lat, lon, scaling)
+    x = scaling * deg2rad(lon)
+    y = scaling * log(tand((90 + lat)/2))
+    (x,y)
+end
+
+function web_mercator_reverse(x, y, scaling)
+    lon = rad2deg(x / scaling)
+    lat = 2 * atand(exp(y / scaling)) - 90
+    (lat,lon)
+end

test/transformations.jl

@@ -73,6 +73,29 @@
     @test inv(enu_ecef) == ecef_enu
     @test inv(ecef_enu) == enu_ecef
 
+    # Web Mercator
+    points = [LLA(0, 0, 0),
+              LLA(49, 2, 0),
+			  LLA(-19, -128, 4),
+			  LLA(36, -77, 6)]
+    # Reference generated using Proj4
+    # [Proj4.transform(Proj4.Projection("+proj=longlat +datum=WGS84"),
+    #                  Proj4.Projection("+proj=webmerc +datum=WGS84"),
+    #                  [point.lon, point.lat, point.alt]) for point in points]
+    points_wm = [[0.0, 0.0, 0.0],
+                 [222638.98158654713, 6.274861394006577e6, 0.0],
+                 [-1.4248894821539015e7, -2.1549359150858945e6, 4.0],
+                 [-8.571600791082064e6, 4.300621372044271e6, 6.0]]
+    wm_lla = WebMercatorfromLLA()
+    lla_wm = LLAfromWebMercator()
+    @testset "WebMercator" for (lla,wm) in zip(points, points_wm)
+        @test wm_lla(lla) ≈ wm
+        @test lla_wm(wm) ≈ lla
+        ll = LatLon(lla.lat, lla.lon)
+        wm_2d = wm[1:2]
+        @test wm_lla(ll) ≈ wm_2d
+        @test lla_wm(wm_2d) ≈ ll
+	end
 
     # Composed transformation functions
     @test LLAfromENU(ecef, wgs84) == LLAfromECEF(wgs84) ∘ ECEFfromENU(ecef, wgs84)