source: trunk/npemap.org.uk/scripts/exporter/Geo/HelmertTransform.pm @ 148

Last change on this file since 148 was 148, checked in by Nick Burch, 15 years ago

Make a start on an exporter

File size: 8.7 KB
Line 
1#!/usr/bin/perl
2#
3# Geo/HelmertTransform.pm:
4# Perform "Helmert" (linear) transformations between coordinates referenced to
5# different datums.
6#
7# Reference:
8#   http://www.gps.gov.uk/additionalInfo/images/A_guide_to_coord.pdf
9#
10# Copyright (c) 2005 UK Citizens Online Democracy.  This module is free
11# software; you can redistribute it and/or modify it under the same terms as
12# Perl itself.
13#
14# Email: chris@mysociety.org; WWW: http://www.mysociety.org/
15#
16# $Id: HelmertTransform.pm,v 1.6 2006/06/21 17:10:24 francis Exp $
17#
18
19package Geo::HelmertTransform;
20
21use strict;
22
23=head1 NAME
24
25Geo::HelmertTransform
26
27=head1 SYNOPSIS
28
29    use Geo::HelmertTransform;
30
31    my ($lat, $lon, $height) = ...; # from OS map
32    my $airy1830 = Geo::HelmertTransform::datum('Airy1830');
33    my $wgs84 = Geo::HelmertTransform::datum('WGS84');
34
35    ($lat, $lon, $height)
36        = Geo::HelmertTransform::convert_datum($airy1830, $wgs84,
37                                                $lat, $lon, $h);
38   
39
40=head1 DESCRIPTION
41
42Perform transformations between geographical coordinates in different datums.
43
44It is usual to describe geographical points in terms of their polar coordinates
45(latitude, longitude and altitude) referenced to a "datum ellipsoid", which is
46used to approximate the Earth's geoid. The latitude, longitude and altitude of
47a given physical point vary depending on which datum ellipsoid is in use.
48Unfortunately, a number of ellipsoids are in everyday use, and so it is often
49necessary to transform geographical coordinates between different datum
50ellipsoids.
51
52Two different datum ellipsoids may differ in the locations of their centers, or
53in their shape; and there may be an angle between their equatorial planes or
54the meridians relative to which longitude is measured. The Helmert Transform,
55which this module implements, is a linear transformation of coordinates between
56pairs of datum ellipsoids in the limit of small angles of deviation between
57them.
58
59=head1 CONVENTIONS
60
61Latitude is expressed in degrees, positive-north; longitude in degrees,
62positive-east. Heights (ellipsoid) and cartesian coordinates are in meters.
63
64=head1 FUNCTIONS
65
66=over 4
67
68=cut
69
70use constant M_PI => 3.141592654;
71
72=item rad_to_deg RADIANS
73
74Convert RADIANS to degrees.
75
76=cut
77sub rad_to_deg ($) {
78    return 180. * $_[0] / M_PI;
79}
80
81=item deg_to_rad DEGREES
82
83Convert DEGREES to radians.
84
85=cut
86sub deg_to_rad ($) {
87    return M_PI * $_[0] / 180.;
88}
89
90=item geo_to_xyz DATUM LAT LON H
91
92Return the Cartesian (X, Y, Z) coordinates for the geographical coordinates
93(LAT, LON, H) in the given DATUM.
94
95=cut
96sub geo_to_xyz ($$$$) {
97    my ($datum, $lat, $lon, $h) = @_;
98    $lat = deg_to_rad($lat);
99    $lon = deg_to_rad($lon);
100   
101    my $v = $datum->a() / sqrt(1 - $datum->e2() * sin($lat) ** 2);
102    return (
103            ($v + $h) * cos($lat) * cos($lon),
104            ($v + $h) * cos($lat) * sin($lon),
105            ((1 - $datum->e2()) * $v + $h) * sin($lat)
106        );
107}
108
109=item xyz_to_geo DATUM X Y Z
110
111Return the geographical (LAT, LON, H) coordinates for the Cartesian coordinates
112(X, Y, Z) in the given DATUM. This is an iterative procedure.
113
114=cut
115sub xyz_to_geo ($$$$) {
116    my ($datum, $x, $y, $z) = @_;
117    my ($lat, $lat2, $lon, $h, $v, $p);
118    $lon = atan2($y, $x);
119   
120    $p = sqrt($x**2 + $y**2);
121    $lat2 = atan2($z, $p);
122
123    my $niter = 0;
124    do {
125        $lat = $lat2;
126        $v = $datum->a() / sqrt(1 - $datum->e2() * sin($lat) ** 2);
127        $lat2 = atan2(($z + $datum->e2() * $v * sin($lat)), $p);
128        die "exceeded 10000 iterations without converging in Geo::HelmertTransform::xyz_to_geo"
129            if (++$niter > 10000);
130    } while (abs($lat2 - $lat) > 2e-6); # about 1/10000 mile
131
132    $h = $p / cos($lat) - $v;
133
134    return (rad_to_deg($lat), rad_to_deg($lon), $h);
135}
136
137=item convert_datum D1 D2 LAT LON H
138
139Given geographical coordinates (LAT, LON, H) in datum D1, return the
140corresponding coordinates in datum D2. This assumes that the transformations
141are small, and always converts via WGS84.
142
143=cut
144sub convert_datum ($$$$$) {
145    my ($d1, $d2, $lat, $lon, $h) = @_;
146    my ($x1, $y1, $z1) = geo_to_xyz($d1, $lat, $lon, $h);
147    my ($x, $y, $z) = ($x1, $y1, $z1);
148    if (!$d1->is_wgs84()) {
149        # Transform into WGS84.
150        $x = $d1->tx()
151                + (1 + $d1->s()) * $x1
152                - $d1->rz()      * $y1
153                + $d1->ry()      * $z1;
154        $y = $d1->ty()
155                + $d1->rz()      * $x1
156                + (1 + $d1->s()) * $y1
157                - $d1->rx()      * $z1;
158        $z = $d1->tz()
159                - $d1->ry()      * $x1
160                + $d1->rx()      * $y1
161                + (1 + $d1->s()) * $z1;
162    }
163
164    my ($x2, $y2, $z2) = ($x, $y, $z);
165    if (!$d2->is_wgs84()) {
166        $x2 = -$d2->tx()
167                + (1 - $d2->s()) * $x
168                + $d2->rz()      * $y
169                - $d2->ry()      * $z;
170        $y2 = -$d2->ty()
171                - $d2->rz()      * $x
172                + (1 - $d2->s()) * $y
173                + $d2->rx()      * $z;
174        $z2 = -$d2->tz()
175                + $d2->ry()      * $x
176                - $d2->rx()      * $y
177                + (1 - $d2->s()) * $z;
178    }
179
180    return xyz_to_geo($d2, $x2, $y2, $z2);
181}
182
183=item datum NAME
184
185Return the datum of the given NAME. Currently implemented are:
186
187=over 4
188
189=item Airy1830
190
191The 1830 Airy ellipsoid to which the British Ordnance Survey's National Grid is
192referenced.
193
194=item Airy1830Modified
195
196The modified 1830 Airy ellipsoid to which the Irish Grid (as used by Ordnance
197Survey Ireland and Ordnance Survey Northern Ireland); also known as the Ireland
1981975 datum.
199
200=item WGS84
201
202The global datum used for GPS.
203
204=back
205
206=cut
207sub datum ($) {
208    return new Geo::HelmertTransform::Datum(Name => $_[0]);
209}
210
211=back
212
213=cut
214
215# Datum class for internal use (alternative spelling: "I can't be bothered to
216# document it now").
217package Geo::HelmertTransform::Datum;
218
219use fields qw(name a b e2 tx ty tz s rx ry rz is_wgs84);
220
221# Fields are: semi-major and -minor axes; and the x-, y-, and z-displacements,
222# scale change, and rotations to transform from this datum into WGS84.
223#
224#                             a (m)          b               tx        ty        tz        s (ppm)   rx (sec) ry       rz
225#                             -------------- --------------- --------- --------- --------- --------- -------- -------- -------
226my %known_datums = (
227            # from OS article above
228        Airy1830          => [6_377_563.396, 6_356_256.910,  +446.448, -125.157, +542.060, -20.4894, +0.1502, +0.2470, +0.8421],
229            # from http://www.osni.gov.uk/downloads/Making%20maps%20GPS%20compatible.pdf
230        Airy1830Modified  => [6_377_340.189, 6_356_034.447,  +482.530, -130.596, +564.557,  +8.150,  +1.042,  +0.214,  +0.631],
231#        International1924 => [6_378_388.000, 6_356_911.946,  ??? ],
232        WGS84             => [6_378_137.000, 6_356_752.3141,   0.000,    0.000,    0.000,   0.0000,  0.0000,  0.0000,  0.0000]
233    );
234
235sub new ($%) {
236    my ($class, %p) = @_;
237    if (exists($p{Name})) {
238        die "datum \"$p{Name}\" not known"
239            if (!exists($known_datums{$p{Name}}));
240
241        # Take a copy of the datum data
242        my $reald = $known_datums{$p{Name}};
243        my $d;
244        foreach my $df (@$reald) {
245            push @$d,$df;
246        }
247
248        my $s = fields::new($class);
249        foreach (qw(a b tx ty tz)) {
250            $s->{$_} = shift(@$d);
251        }
252        $s->{s} = shift(@$d) / 1_000_000;               # ppm
253        foreach (qw(rx ry rz)) {
254            $s->{$_} = Geo::HelmertTransform::deg_to_rad(shift(@$d) / 3600.);  # seconds
255        }
256        $s->{is_wgs84} = ($p{Name} eq 'WGS84');
257        return $s;
258    } elsif (!exists($p{a}) || !exists($p{b})) {
259        die "must specify semi-major axis a and semi-minor axis b";
260    } else {
261        my $s = fields::new($class);
262        foreach (qw(a b tx ty tz s rx ry rz)) {
263            $s->{$_} = 0;
264            $s->{$_} = $p{$_} if (exists($p{$_}));
265        }
266        $s->{is_wgs84} = 0;
267        return $s;
268    }
269}
270
271foreach (qw(a b tx ty tz s rx ry rz is_wgs84)) {
272    eval <<EOF;
273sub $_ (\$) {
274    return \$_[0]->{$_};
275}
276EOF
277}
278
279sub e2 ($) {
280    my $s = shift;
281    if (!exists($_[0]->{e2})) {
282        if($s->a() == 0) {
283            $s->{e2} = 1;
284        } else {
285            $s->{e2} = 1 - ($s->b() / $s->a()) ** 2;
286        }
287    }
288    return $s->{e2}
289}
290
291=head1 SEE ALSO
292
293I<A guide to coordinate systems in Great Britain>,
294http://www.gps.gov.uk/guidecontents.asp
295
296I<Making maps compatible with GPS>,
297http://www.osni.gov.uk/downloads/Making%20maps%20GPS%20compatible.pdf
298
299=head1 AUTHOR AND COPYRIGHT
300
301Written by Chris Lightfoot, chris@mysociety.org
302
303Copyright (c) UK Citizens Online Democracy.
304
305This module is free software; you can redistribute it and/or modify it under
306the same terms as Perl itself.
307
308=head1 VERSION
309
310$Id: HelmertTransform.pm,v 1.6 2006/06/21 17:10:24 francis Exp $
311
312=cut
313
3141;
Note: See TracBrowser for help on using the repository browser.