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
// examples/00-transformations.rs
// Using Rust Geodesy to transform geodata.
// Run with:
// cargo run --example 00-transformations
use geodesy::prelude::*;
// Use Anyhow for convenient error handling
fn main() -> anyhow::Result<()> {
// The context provider is the entry point to all transformation functionality:
let mut ctx = Minimal::default();
// The concept of a "context data structure" will be well known to
// PROJ users, where the context plays a somewhat free-flowing role,
// and only becomes truly visible in multithreaded cases.
// In Rust Geodesy, the context plays a much more visible role, as
// most transformation functionality is implemented directly as
// methods of the context data structure.
// We need some coordinates to test the code. The convenience methods
// `gis` and `geo` produce 2D coordinate tuples and automatically handle
// conversion of the angular parts to radians, and geographical coordinates
// in latitude/longitude order, to the GIS convention of longitude/latitude.
// Here using the GIS convention - longitude before latitude
let cph = Coor2D::gis(12., 55.); // Copenhagen
let osl = Coor2D::gis(10., 60.); // Oslo
// And here using the geodesy/navigation convention - latitude before longitude
let sth = Coor2D::geo(59., 18.); // Stockholm
let hel = Coor2D::geo(60., 25.); // Helsinki
// `gis` and `geo` have a sibling `raw` which generates coordinate tuples
// from raw numbers, in case your point coordinates are already given in
// radians.
let cph_raw = Coor2D::raw(12_f64.to_radians(), 55_f64.to_radians());
// But since a coordinate tuple is really just an array of double
// precision numbers, you may also generate it directly using plain
// Rust syntax. Note that Coor2D, like f64, provides the to_radians
// method, but it operates in place, so we need two steps in this case.
let cph_direct = Coor2D([12., 55.]).to_radians();
// The three versions of Copenhagen coordinates should be identical.
assert_eq!(cph, cph_raw);
assert_eq!(cph, cph_direct);
// The Rust Geodesy interface is based on transformation of *arrays* of
// coordinate tuples, rather than single points. So let's make an array:
let mut data = [osl, cph, sth, hel];
// Since all operations are carried out in place, the array needs to
// be mutable, hence `let mut`.
// Note that on the inside, Rust Geodesy operates in 4 dimensions, but
// externally you may use whatever fits the problem at hand. Here, we
// use the two-dimensional Coor2D, since we're doing a plain projection
// from 2D geographical coordinates (latitude, longitude) to 2D projected
// coordinates (northing, easting)
// Let's obtain a handle to a transformation element ("an operation"),
// turning geographical coordinates into UTM zone 32 coordinates.
// Since this may go wrong (e.g. due to syntax errors in the operator
// definition), use the Rust `?`-operator to handle errors.
let utm32 = ctx.op("utm zone=32")?;
// Now, let's use the utm32-operator to transform some coordinate data.
// The data are transformed in place, so we pass them by mutable reference.
// The `Fwd` constant indicates that the operation runs in the forward
// direction
ctx.apply(utm32, Fwd, &mut data)?;
println!("utm32:");
println!(" {:#?}", data);
for coord in data {
println!(" {:?}", coord);
}
// Specifying `Inv` rather than `Fwd` takes us on the inverse trip back
// to geographic coordinates
ctx.apply(utm32, Inv, &mut data)?;
println!("Roundtrip to geo:");
// Note the use of `to_geo`, which transforms lon/lat in radians
// to lat/lon in degrees. It is defined for Coor2D as well as for
// arrays, vectors and slices of Coor2D
for coord in data {
println!(" {:?}", coord.to_geo());
}
// To geo again, but using slices - in two different ways
println!("Sliced to_geo:");
let mut data = [osl, cph, sth, hel];
let slice = &mut data[..2];
for coord in slice {
println!(" {:?}", coord.to_geo());
}
for coord in (&mut data[2..]).iter_mut() {
println!(" {:?}", coord.to_geo());
}
// ctx.apply(...) with slices (a step towards parallel operation)
println!("Sliced utm32:");
let mut data = [osl, cph, sth, hel];
let slice = &mut data[..];
let (mut first, mut last) = slice.split_at_mut(2);
ctx.apply(utm32, Fwd, &mut first)?;
ctx.apply(utm32, Fwd, &mut last)?;
for coord in data {
println!(" {:?}", coord);
}
// Now a slightly more complex case: Transforming the coordinates,
// which we consider given in WGS84, back to the older ED50 datum.
// The EPSG:1134 method handles that through a 3 parameter Helmert
// transformation. But since the Helmert transformation works on
// cartesian coordinates, rather than geographic, we need to add
// pre- and post-processing steps, taking us from geographical
// coordinates to cartesian, and back. Hence, we need a pipeline
// of 3 steps:
let pipeline = "cart ellps=intl | helmert x=-87 y=-96 z=-120 | cart inv=true ellps=GRS80";
// Here we're doing a 3D operation, so we need a 3D coordinate representation
let cph = Coor3D::geo(55., 12., 0.); // Copenhagen
let osl = Coor3D::geo(60., 10., 0.); // Oslo
let sth = Coor3D::geo(59., 18., 0.); // Stockholm
let hel = Coor3D::geo(60., 25., 0.); // Helsinki
let mut data = [osl, cph, sth, hel];
let ed50_wgs84 = ctx.op(pipeline)?;
// Since the forward transformation goes *from* ed50 to wgs84, we use
// the inverse method to take us the other way, back in time to ED50
ctx.apply(ed50_wgs84, Inv, &mut data)?;
println!("ed50:");
for coord in data {
println!(" {:?}", coord.to_geo());
}
// Finally an example of handling bad syntax:
println!("Bad syntax example:");
let op = ctx.op("aargh zone=23");
if op.is_err() {
println!("Deliberate error - {:?}", op);
}
// ------------------------------------------------------------------------------
// Don't do this!
// ------------------------------------------------------------------------------
// As PROJ-conoisseurs know, the `Context` type in PROJ plays a somewhat more
// passive role than the `Context` type in RG. But actually, the RG Context API
// is built on top of a much more PROJ-like interface, where the operators take
// center stage. In general, the use of this interface is *not* recommended, but
// for the rare cases where it is preferable, we demonstrate its use below, by
// repeating the exercises above, while swapping the roles of the `Op`/`OpHandle`
// and the `Context`.
println!("------------------------------------");
println!(" Doing the same in the wrong way... ");
println!("------------------------------------");
use geodesy::authoring::Op;
// Create an `Op`, turning geographical coordinates into UTM zone 32 coordinates
let utm32 = Op::new("utm zone=32", &ctx)?;
// Now, let's use the utm32-operator to transform some data
utm32.apply(&ctx, &mut data, Fwd);
// println!("utm32:");
// for coord in data {
// println!(" {:?}", coord);
// }
// Take the inverse road back to geographic coordinates
utm32.apply(&ctx, &mut data, Inv);
// println!("Roundtrip to geo:");
// println!(" {:#?}", data);
// for coord in data {
// println!(" {:?}", coord.to_geo());
// }
// EPSG:1134
let ed50_wgs84 = Op::new(pipeline, &ctx)?;
let mut data = [osl, cph, sth, hel];
ed50_wgs84.apply(&ctx, &mut data, Inv);
// println!("ed50:");
// println!(" {:?}", data.to_geo());
// Handling bad syntax:
println!("Bad syntax example:");
let op = Op::new("aargh zone=23", &ctx);
//let op = ctx.define_operation("aargh zone: 23");
if op.is_err() {
println!("Deliberate error - {:?}", op);
}
Ok(())
}