1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
/**
* \file OSGB.hpp
* \brief Header for GeographicLib::OSGB class
*
* Copyright (c) Charles Karney (2010-2024) <karney@alum.mit.edu> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_OSGB_HPP)
#define GEOGRAPHICLIB_OSGB_HPP 1
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/TransverseMercator.hpp>
namespace GeographicLib {
/**
* \brief Ordnance Survey grid system for Great Britain
*
* The class implements the coordinate system used by the Ordnance Survey for
* maps of Great Britain and conversions to the grid reference system.
*
* See
* - <a href="https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf">
* A guide to coordinate systems in Great Britain</a>
* - <a href="https://www.ordnancesurvey.co.uk/documents/resources/guide-to-nationalgrid.pdf">
* Using the National Grid</a>
*
* \warning the latitudes and longitudes for the Ordnance Survey grid
* system do not use the WGS84 datum. Do not use the values returned by this
* class in the UTMUPS, MGRS, or Geoid classes without first converting the
* datum (and vice versa).
*
* Example of use:
* \include example-OSGB.cpp
**********************************************************************/
class GEOGRAPHICLIB_EXPORT OSGB {
private:
typedef Math::real real;
static const char* const letters_;
static const char* const digits_;
static const TransverseMercator& OSGBTM();
static constexpr int base_ = 10;
static constexpr int tile_ = 100000;
static constexpr int tilelevel_ = 5;
static constexpr int tilegrid_ = 5;
static constexpr int tileoffx_ = 2 * tilegrid_;
static constexpr int tileoffy_ = 1 * tilegrid_;
static constexpr int minx_ = - tileoffx_ * tile_;
static constexpr int miny_ = - tileoffy_ * tile_;
static constexpr int maxx_ = (tilegrid_*tilegrid_ - tileoffx_) * tile_;
static constexpr int maxy_ = (tilegrid_*tilegrid_ - tileoffy_) * tile_;
// Maximum precision is um
static constexpr int maxprec_ = 5 + 6;
static real computenorthoffset();
static void CheckCoords(real x, real y);
OSGB() = delete; // Disable constructor
public:
/**
* Forward projection, from geographic to OSGB coordinates.
*
* @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees).
* @param[out] x easting of point (meters).
* @param[out] y northing of point (meters).
* @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
*
* \e lat should be in the range [−90°, 90°].
**********************************************************************/
static void Forward(real lat, real lon,
real& x, real& y, real& gamma, real& k) {
OSGBTM().Forward(OriginLongitude(), lat, lon, x, y, gamma, k);
x += FalseEasting();
y += computenorthoffset();
}
/**
* Reverse projection, from OSGB coordinates to geographic.
*
* @param[in] x easting of point (meters).
* @param[in] y northing of point (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
*
* The value of \e lon returned is in the range [−180°,
* 180°].
**********************************************************************/
static void Reverse(real x, real y,
real& lat, real& lon, real& gamma, real& k) {
x -= FalseEasting();
y -= computenorthoffset();
OSGBTM().Reverse(OriginLongitude(), x, y, lat, lon, gamma, k);
}
/**
* OSGB::Forward without returning the convergence and scale.
**********************************************************************/
static void Forward(real lat, real lon, real& x, real& y) {
real gamma, k;
Forward(lat, lon, x, y, gamma, k);
}
/**
* OSGB::Reverse without returning the convergence and scale.
**********************************************************************/
static void Reverse(real x, real y, real& lat, real& lon) {
real gamma, k;
Reverse(x, y, lat, lon, gamma, k);
}
/**
* Convert OSGB coordinates to a grid reference.
*
* @param[in] x easting of point (meters).
* @param[in] y northing of point (meters).
* @param[in] prec precision relative to 100 km.
* @param[out] gridref National Grid reference.
* @exception GeographicErr if \e prec, \e x, or \e y is outside its
* allowed range.
* @exception std::bad_alloc if the memory for \e gridref can't be
* allocatied.
*
* \e prec specifies the precision of the grid reference string as follows:
* - prec = 0 (min), 100km
* - prec = 1, 10km
* - prec = 2, 1km
* - prec = 3, 100m
* - prec = 4, 10m
* - prec = 5, 1m
* - prec = 6, 0.1m
* - prec = 11 (max), 1μm
*
* The easting must be in the range [−1000 km, 1500 km) and the
* northing must be in the range [−500 km, 2000 km). These bounds
* are consistent with rules for the letter designations for the grid
* system.
*
* If \e x or \e y is NaN, the returned grid reference is "INVALID".
**********************************************************************/
static void GridReference(real x, real y, int prec, std::string& gridref);
/**
* Convert OSGB grid reference to coordinates.
*
* @param[in] gridref National Grid reference.
* @param[out] x easting of point (meters).
* @param[out] y northing of point (meters).
* @param[out] prec precision relative to 100 km.
* @param[in] centerp if true (default), return center of the grid square,
* else return SW (lower left) corner.
* @exception GeographicErr if \e gridref is illegal.
*
* The grid reference must be of the form: two letters (not including I)
* followed by an even number of digits (up to 22).
*
* If the first 2 characters of \e gridref are "IN", then \e x and \e y are
* set to NaN and \e prec is set to −2.
**********************************************************************/
static void GridReference(const std::string& gridref,
real& x, real& y, int& prec,
bool centerp = true);
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the Airy 1830 ellipsoid (meters).
*
* This is 20923713 ft converted to meters using the rule 1 ft =
* 10<sup>9.48401603−10</sup> m. The Airy 1830 value is returned
* because the OSGB projection is based on this ellipsoid. The conversion
* factor from feet to meters is the one used for the 1936 retriangulation
* of Britain; see Section A.1 (footnote 10 on p. 44) of <i>A guide to
* coordinate systems in Great Britain</i>, v3.6 (2020).
**********************************************************************/
static Math::real EquatorialRadius() {
// result is about 6377563.3960320664406 m
using std::pow;
return pow(real(10), real(48401603 - 100000000) / 100000000)
* real(20923713);
}
/**
* @return \e f the inverse flattening of the Airy 1830 ellipsoid.
*
* For the Airy 1830 ellipsoid, \e a = 20923713 ft and \e b = 20853810 ft;
* thus the flattening = (20923713 − 20853810)/20923713 =
* 7767/2324857 = 1/299.32496459... (The Airy 1830 value is returned
* because the OSGB projection is based on this ellipsoid.)
**********************************************************************/
static Math::real Flattening()
{ return real(20923713 - 20853810) / real(20923713); }
/**
* @return \e k0 central scale for the OSGB projection (0.9996012717...).
*
* C. J. Mugnier, Grids & Datums, PE&RS, Oct. 2003, states that
* this is defined as 10<sup>9.9998268−10</sup>.
**********************************************************************/
static Math::real CentralScale() {
using std::pow;
return pow(real(10), real(9998268 - 10000000) / 10000000);
}
/**
* @return latitude of the origin for the OSGB projection (49 degrees).
**********************************************************************/
static Math::real OriginLatitude() { return real(49); }
/**
* @return longitude of the origin for the OSGB projection (−2
* degrees).
**********************************************************************/
static Math::real OriginLongitude() { return real(-2); }
/**
* @return false northing the OSGB projection (−100000 meters).
**********************************************************************/
static Math::real FalseNorthing() { return real(-100000); }
/**
* @return false easting the OSGB projection (400000 meters).
**********************************************************************/
static Math::real FalseEasting() { return real(400000); }
///@}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_OSGB_HPP