view release on metacpan or  search on metacpan

lib/Geo/HelmertTransform.pm  view on Meta::CPAN

=cut

use constant M_PI => 3.141592654;

=item rad_to_deg RADIANS

Convert RADIANS to degrees.

=cut
sub rad_to_deg ($) {
    return 180. * $_[0] / M_PI;
}

=item deg_to_rad DEGREES

Convert DEGREES to radians.

=cut
sub deg_to_rad ($) {
    return M_PI * $_[0] / 180.;
}

=item geo_to_xyz DATUM LAT LON H

Return the Cartesian (X, Y, Z) coordinates for the geographical coordinates
(LAT, LON, H) in the given DATUM.

=cut
sub geo_to_xyz ($$$$) {
    my ($datum, $lat, $lon, $h) = @_;
    $lat = deg_to_rad($lat);
    $lon = deg_to_rad($lon);
    
    my $v = $datum->a() / sqrt(1 - $datum->e2() * sin($lat) ** 2);
    return (
            ($v + $h) * cos($lat) * cos($lon),
            ($v + $h) * cos($lat) * sin($lon),
            ((1 - $datum->e2()) * $v + $h) * sin($lat)
        );
}

=item xyz_to_geo DATUM X Y Z

Return the geographical (LAT, LON, H) coordinates for the Cartesian coordinates
(X, Y, Z) in the given DATUM. This is an iterative procedure.

=cut
sub xyz_to_geo ($$$$) {
    my ($datum, $x, $y, $z) = @_;
    my ($lat, $lat2, $lon, $h, $v, $p);
    $lon = atan2($y, $x);
    
    $p = sqrt($x**2 + $y**2);
    $lat2 = atan2($z, $p);

    my $niter = 0;
    do {
        $lat = $lat2;
        $v = $datum->a() / sqrt(1 - $datum->e2() * sin($lat) ** 2);
        $lat2 = atan2(($z + $datum->e2() * $v * sin($lat)), $p);
        die "exceeded 10000 iterations without converging in Geo::HelmertTransform::xyz_to_geo"
            if (++$niter > 10000);
    } while (abs($lat2 - $lat) > 2e-6); # about 1/10000 mile

    $h = $p / cos($lat) - $v;

    return (rad_to_deg($lat), rad_to_deg($lon), $h);
}

=item convert_datum D1 D2 LAT LON H

Given geographical coordinates (LAT, LON, H) in datum D1, return the
corresponding coordinates in datum D2. This assumes that the transformations
are small, and always converts via WGS84.

=cut
sub convert_datum ($$$$$) {
    my ($d1, $d2, $lat, $lon, $h) = @_;
    my ($x1, $y1, $z1) = geo_to_xyz($d1, $lat, $lon, $h);
    my ($x, $y, $z) = ($x1, $y1, $z1);
    if (!$d1->is_wgs84()) {
        # Transform into WGS84.
        $x = $d1->tx()
                + (1 + $d1->s()) * $x1
                - $d1->rz()      * $y1
                + $d1->ry()      * $z1;
        $y = $d1->ty()
                + $d1->rz()      * $x1
                + (1 + $d1->s()) * $y1
                - $d1->rx()      * $z1;
        $z = $d1->tz()
                - $d1->ry()      * $x1
                + $d1->rx()      * $y1
                + (1 + $d1->s()) * $z1;
    }

    my ($x2, $y2, $z2) = ($x, $y, $z);
    if (!$d2->is_wgs84()) {
        $x2 = -$d2->tx()
                + (1 - $d2->s()) * $x
                + $d2->rz()      * $y
                - $d2->ry()      * $z;
        $y2 = -$d2->ty()
                - $d2->rz()      * $x
                + (1 - $d2->s()) * $y
                + $d2->rx()      * $z;
        $z2 = -$d2->tz()
                + $d2->ry()      * $x
                - $d2->rx()      * $y
                + (1 - $d2->s()) * $z;
    }

    return xyz_to_geo($d2, $x2, $y2, $z2);
}

=item datum NAME

Return the datum of the given NAME. Currently implemented are:

=over 4