Standard Molodensky Datum Transform Source Code

Below is a listing of a Java implementation of the standard Molodensky datum transformation method.

Peter H. Dana has an excellent discussion of geodetic datums at http://www.colorado.edu/geography/gcraft/notes/datum/datum_f.html. A diagram showing the Molodensky formulas is located at http://www.colorado.edu/geography/gcraft/notes/datum/gif/molodens.gif, and an extensive table of datums and ellipsoid parameters are at http://www.colorado.edu/geography/gcraft/notes/datum/edlist.html.

For information on transformations that yield results with higher accuracy than Molodensky, see the document High-Accuracy Datum Transformations.

Source Listing

In the listing below, the class GeodeticPosition has three members, lon, lat, and h. They are double-precision values indicating the longitude and latitude in radians, and height in meters above the ellipsoid.

The source code in the listing below may be copied and reused without restriction, but it is offered AS-IS with NO WARRANTY.

    /*
    * transform
    *
    * Parameters:
    *     from:     The geodetic position to be translated.
    *     from_a:   The semi-major axis of the "from" ellipsoid.
    *     from_f:   Flattening of the "from" ellipsoid.
    *     from_esq: Eccentricity-squared of the "from" ellipsoid.
    *     da:       Change in semi-major axis length (meters); "to" minus "from"    *     df:       Change in flattening; "to" minus "from"
    *     dx:       Change in x between "from" and "to" datum.
    *     dy:       Change in y between "from" and "to" datum.
    *     dz:       Change in z between "from" and "to" datum.
    */
    public GeodeticPosition transform (GeodeticPosition from,
                                       double from_a, double from_f,

                                       double from_esq, double da, double df,
                                       double dx, double dy, double dz)
    {
        double slat = Math.sin (from.lat);
        double clat = Math.cos (from.lat);
        double slon = Math.sin (from.lon);
        double clon = Math.cos (from.lon);
        double ssqlat = slat * slat;
        double adb = 1.0 / (1.0 - from_f);  // "a divided by b"
        double dlat, dlon, dh;

        double rn = from_a / Math.sqrt (1.0 - from_esq * ssqlat);
        double rm = from_a * (1. - from_esq) / Math.pow ((1.0 - from_esq * ssqlat), 1.5);

        dlat = (((((-dx * slat * clon - dy * slat * slon) + dz * clat)
                    + (da * ((rn * from_esq * slat * clat) / from_a)))
                + (df * (rm * adb + rn / adb) * slat * clat)))
            / (rm + from.h);

        dlon = (-dx * slon + dy * clon) / ((rn + from.h) * clat);

        dh = (dx * clat * clon) + (dy * clat * slon) + (dz * slat)
             - (da * (from_a / rn)) + ((df * rn * ssqlat) / adb);

        return new GeodeticPosition (from.lon + dlon, from.lat + dlat, from.h + dh);
    }

Chuck Taylor -- (Copyright) -- (Contact)