Code covered by the BSD License  

Highlights from
Geodesic projections for an ellipsoid

image thumbnail

Geodesic projections for an ellipsoid



07 Dec 2012 (Updated )

Four map projections based on geodesics

function geodproj
%GEODPROJ  Geodesic projections for an ellipsoid
%   This package implements four projections based on geodesics:
%     * the azimuthal equidistant projection (eqdazim)
%     * the Cassini-Soldner projection (cassini)
%     * the transverse Mercator projection (tranmerc)
%     * the ellipsoidal gnomonic projection (gnomonic)
%   The package implements the forward projection (from geographic to
%   projected coordinates) and inverse projection (from projected to
%   geographic coordinates) with abbreviated function names (listed in
%   parentheses in the list above) suffixed by _FWD and _INV.  In addition,
%   implementations of the UTM projection are provided with UTM_FWD and
%   UTM_INV.  For each function, metric properties of the projection are
%   also returned.
%   This package requires the availability of MATLAB File Exchange package
%   "Geodesics on an ellipsoid of revolution" for performing the necessary
%   geodesic computations:
%   The azimuthal equidistant inverse projection is a geodesic projection
%   defined by
%     X = S.*SIND(AZI0), Y = S.*COSD(AZI0)
%   Thus all distances and azimuths relative to the center point are
%   correct.
%   The Cassini-Solder projection is a rectangular geodesic projection
%   whose inverse is specified by
%   Cassini-Soldner is a transverse cylindrical projection where the
%   meridian LON0 maps to a straight line with constant scale and distances
%   perpendicular to the central meridian are true.  Cassini-Soldner was
%   widely used for large scale maps until about 1930 (when it was
%   supplanted by various conformal projections, in particular, by the
%   transverse Mercator projection).
%   The transverse Mercator projection is also a cylindrical projection
%   which maps the central meridian to a straight line at constant scale.
%   However the behavior either side of the central meridian is determined
%   by the condition of conformality.  The implementation used in this
%   package is based on the series method described in
%     C. F. F. Karney, Transverse Mercator with an accuracy of a few
%     nanometers, J. Geodesy 85(8), 475-485 (Aug. 2011);
%     Addenda:
%   The ellipsoidal gnomonic projection is defined by
%     [~,AZI0,~,~,m,M] = GEODDISTANCE(LAT0,LON0,LAT,LON)
%     RHO = m./M, X = RHO.*SIND(AZI0), Y = RHO.*COSD(AZI0)
%   Obviously this is an azimuthal projection.  It also enjoys approximately
%   the property of the spherical gnomonic projection, that geodesics map
%   to straight lines.  The projection is derived in Section 8 of
%     C. F. F. Karney, Algorithms for geodesics,
%     J. Geodesy 87, 43-55 (2013);
%     Addenda:
%   The parameters of the ellipsoid are specified by the optional ELLIPSOID
%   argument to the routines.  This is a two-element vector of the form
%   [a,e], where a is the equatorial radius, e is the eccentricity e =
%   sqrt(a^2-b^2)/a, and b is the polar semi-axis.  Typically, a and b are
%   measured in meters and the linear and area quantities returned by the
%   routines are then in meters and meters^2.  However, other units can be
%   employed.  If ELLIPSOID is omitted, then the WGS84 ellipsoid (more
%   precisely, the value returned by DEFAULTELLIPSOID) is assumed [6378137,
%   0.0818191908426215] corresponding to a = 6378137 meters and a
%   flattening f = (a-b)/a = 1/298.257223563.  The flattening and
%   eccentricity are related by
%       e = sqrt(f * (2 - f))
%       f = e^2 / (1 + sqrt(1 - e^2))
%   (The functions ECC2FLAT and FLAT2ECC implement these conversions.)  For
%   a sphere, set e = 0; for a prolate ellipsoid (b > a), specify e as a
%   pure imaginary number.
%   All angles (latitude, longitude, azimuth) are measured in degrees with
%   latitudes increasing northwards, longitudes increasing eastwards, and
%   azimuths measured clockwise from north.  For a point at a pole, the
%   azimuth is defined by keeping the longitude fixed, writing lat =
%   +/-(90-eps), and taking the limit eps -> 0+.
%   Restrictions on the inputs:
%     * All latitudes must lie in [-90, 90].
%     * All longitudes and azimuths must lie in [-540, 540).  On output,
%       these quantities lie in [-180, 180).
%     * The equatorial radius, a, must be positive.
%     * The eccentricity, e, should be satisfy abs(e) < 0.2 in order to
%       retain full accuracy (this corresponds to flattenings satisfying
%       abs(f) <= 1/50, approximately).  This condition holds for most
%       applications in geodesy.
%   These routines are fully vectorized so that their speed is competitive
%   with the compiled C++ code.  These MATLAB versions are transcriptions
%   of several C++ classes provided by GeographicLib which is available at
%   The routines duplicate some of the functionality of the EQDAZIM,
%   CASSINISTD, TRANMERC, and GNOMONIC projections in the the MATLAB
%   mapping toolbox.  The major improvements offered by this package are
%     * The azimuthal equidistant and gnomonic projections are defined
%       for the ellipsoid (instead of just for the sphere).
%     * The Cassini-Soldner projection is essentially exact (instead of
%       being defined in terms of an approximate series).
%     * The transverse Mercator projection uses a much more accurate
%       series which greatly extends its domain of applicability.
%   The primary importance of the two azimuthal projections is that they
%   offer a convenient way of solving various geometrical problems on the
%   ellipsoid.  In particular the azimuthal equidistant projection allows
%   problems associated with determining maritime boundaries to be solved
%   easily.  Similarly the gnomonic projection allows the intersection of
%   two geodesics to be determined quickly.

% Copyright (c) Charles Karney (2012) <>.
% This file was distributed with GeographicLib 1.29.

  help geodproj

Contact us