# Ticket #89: HelmertTransform.pm

File HelmertTransform.pm, 8.7 KB (added by , 15 years ago) |
---|

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 | |

19 | package Geo::HelmertTransform; |

20 | |

21 | use strict; |

22 | |

23 | =head1 NAME |

24 | |

25 | Geo::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 | |

42 | Perform transformations between geographical coordinates in different datums. |

43 | |

44 | It is usual to describe geographical points in terms of their polar coordinates |

45 | (latitude, longitude and altitude) referenced to a "datum ellipsoid", which is |

46 | used to approximate the Earth's geoid. The latitude, longitude and altitude of |

47 | a given physical point vary depending on which datum ellipsoid is in use. |

48 | Unfortunately, a number of ellipsoids are in everyday use, and so it is often |

49 | necessary to transform geographical coordinates between different datum |

50 | ellipsoids. |

51 | |

52 | Two different datum ellipsoids may differ in the locations of their centers, or |

53 | in their shape; and there may be an angle between their equatorial planes or |

54 | the meridians relative to which longitude is measured. The Helmert Transform, |

55 | which this module implements, is a linear transformation of coordinates between |

56 | pairs of datum ellipsoids in the limit of small angles of deviation between |

57 | them. |

58 | |

59 | =head1 CONVENTIONS |

60 | |

61 | Latitude is expressed in degrees, positive-north; longitude in degrees, |

62 | positive-east. Heights (ellipsoid) and cartesian coordinates are in meters. |

63 | |

64 | =head1 FUNCTIONS |

65 | |

66 | =over 4 |

67 | |

68 | =cut |

69 | |

70 | use constant M_PI => 3.141592654; |

71 | |

72 | =item rad_to_deg RADIANS |

73 | |

74 | Convert RADIANS to degrees. |

75 | |

76 | =cut |

77 | sub rad_to_deg ($) { |

78 | return 180. * $_[0] / M_PI; |

79 | } |

80 | |

81 | =item deg_to_rad DEGREES |

82 | |

83 | Convert DEGREES to radians. |

84 | |

85 | =cut |

86 | sub deg_to_rad ($) { |

87 | return M_PI * $_[0] / 180.; |

88 | } |

89 | |

90 | =item geo_to_xyz DATUM LAT LON H |

91 | |

92 | Return the Cartesian (X, Y, Z) coordinates for the geographical coordinates |

93 | (LAT, LON, H) in the given DATUM. |

94 | |

95 | =cut |

96 | sub 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 | |

111 | Return 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 |

115 | sub 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 | |

139 | Given geographical coordinates (LAT, LON, H) in datum D1, return the |

140 | corresponding coordinates in datum D2. This assumes that the transformations |

141 | are small, and always converts via WGS84. |

142 | |

143 | =cut |

144 | sub 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 | |

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

186 | |

187 | =over 4 |

188 | |

189 | =item Airy1830 |

190 | |

191 | The 1830 Airy ellipsoid to which the British Ordnance Survey's National Grid is |

192 | referenced. |

193 | |

194 | =item Airy1830Modified |

195 | |

196 | The modified 1830 Airy ellipsoid to which the Irish Grid (as used by Ordnance |

197 | Survey Ireland and Ordnance Survey Northern Ireland); also known as the Ireland |

198 | 1975 datum. |

199 | |

200 | =item WGS84 |

201 | |

202 | The global datum used for GPS. |

203 | |

204 | =back |

205 | |

206 | =cut |

207 | sub 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"). |

217 | package Geo::HelmertTransform::Datum; |

218 | |

219 | use 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 | # -------------- --------------- --------- --------- --------- --------- -------- -------- ------- |

226 | my %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 | |

235 | sub 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 | |

271 | foreach (qw(a b tx ty tz s rx ry rz is_wgs84)) { |

272 | eval <<EOF; |

273 | sub $_ (\$) { |

274 | return \$_[0]->{$_}; |

275 | } |

276 | EOF |

277 | } |

278 | |

279 | sub 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 | |

293 | I<A guide to coordinate systems in Great Britain>, |

294 | http://www.gps.gov.uk/guidecontents.asp |

295 | |

296 | I<Making maps compatible with GPS>, |

297 | http://www.osni.gov.uk/downloads/Making%20maps%20GPS%20compatible.pdf |

298 | |

299 | =head1 AUTHOR AND COPYRIGHT |

300 | |

301 | Written by Chris Lightfoot, chris@mysociety.org |

302 | |

303 | Copyright (c) UK Citizens Online Democracy. |

304 | |

305 | This module is free software; you can redistribute it and/or modify it under |

306 | the 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 | |

314 | 1; |