-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
241 lines (189 loc) · 8.4 KB
/
README.Rmd
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
238
239
240
241
---
title: "sgo - Simple Geographical Operations (with OSGB36)"
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
fig.width = 5,
out.width = "50%"
)
```
<!-- badges: start -->
[](https://github.com/clozanoruiz/sgo/actions?workflow=R-CMD-check)
[](https://app.codecov.io/gh/clozanoruiz/sgo?branch=main)
[](https://cran.r-project.org/package=sgo)
[](https://github.com/clozanoruiz/sgo/releases)
[](https://opensource.org/licenses/BSD-2-Clause)
<!-- badges: end -->
The goal of `sgo` is to provide a set of methods to transform ETRS89 coordinates to the equivalent
OSGB36 National Grid coordinates, and vice-versa, using the [Ordnance Survey's National Grid Transformation Model OSTN15](https://www.ordnancesurvey.co.uk/blog/2016/09/ostn15-new-geoid-britain/). Additionally, it also has the ability to calculate distances and areas from groups of points defined in any of the coordinate systems supported.
The Coordinate Systems supported by `sgo` are:
* ETRS89: EPSGs 4258, 4937, 4936 and 3035
* WGS84: EPSGs 4326, 4979, 4978 and 3857
* OSGB36: EPSGs 27700, 7405 and 4277
Please note that this package assumes that the Coordinate Reference Systems
(CRS) ETRS89 and WGS84 are the same within the UK, but this shouldn't be a
problem for most civilian use of GPS satellites. If a high-precision
transformation between WGS84 and ETRS89 is required then it is recommended
to use a different package to do the conversion.
According to the Transformations and OSGM15 User Guide, p. 8:
*"...ETRS89 is a precise version of the better known WGS84 system
optimised for use in Europe; however, for most purposes it can be considered
equivalent to WGS84."* and *"For all navigation, mapping, GIS, and engineering
applications within the tectonically stable parts of Europe (including UK and
Ireland), the term ETRS89 should be taken as synonymous with WGS84."*.
## Installation
You can install the released version of `sgo` from [CRAN](https://CRAN.R-project.org) with:
``` r
install.packages("sgo")
```
And the development version from [GitHub](https://github.com/) with:
``` r
# install.packages("remotes")
remotes::install_github("clozanoruiz/sgo")
```
## Usage
```{r, message = FALSE}
library(sgo)
## sgo_points creation (from lists, matrices or dataframes)
sr <- sgo_points(list(-3.9369, 56.1165), epsg=4326)
sr
ln <- c(-4.22472, -2.09908)
lt <- c(57.47777, 57.14965)
n <- c("Inverness", "Aberdeen")
df <- data.frame(n, ln, lt, stringsAsFactors = FALSE)
locs <- sgo_points(df, coords=c("ln", "lt"), epsg=4326)
library(maps)
map('world', regions=('uk'), xlim=c(-9, 0), ylim=c(54.5, 60.9))
points(x=sr$x, y=sr$y, pch=0, col="green")
points(locs, pch=0, col="red")
text(locs, labels=locs$n, pos=1, cex=0.9)
text(sr, labels="Stirling", pos=1, cex=0.9)
## Convert coordinates to OS National Grid reference
bv <- sgo_points(list(x=247455, y=706338, name="Ben Venue"),
coords=c("x", "y"), epsg=27700)
sgo_bng_ngr(bv)$ngr
# notice the truncating:
sgo_bng_ngr(bv, digits=8)$ngr
## Convert from lon/lat to BNG
lon <- c(-4.25181,-3.18827)
lat <- c(55.86424, 55.95325)
pts <- sgo_points(list(longitude=lon, latitude=lat), epsg=4326)
sgo_lonlat_bng(pts)
## sgo_transform is a wrapper for all the conversions available.
# to BNG:
sgo_transform(locs, to=27700)
# to OSGB36 (historical):
sgo_transform(locs, to=4277)
# to Pseudo-Mercator:
sgo_transform(locs, to=3857)
## Convert OS National Grid references to ETRS89
munros <- data.frame(
ngr=c("NN 1652471183", "NN 98713 98880",
"NN9521499879", "NN96214 97180"),
name=c("Beinn Nibheis", "Beinn Macduibh",
"Am Bràigh Riabhach", "Càrn an t-Sabhail"))
sgo_transform(sgo_ngr_bng(munros, col="ngr"), to=4258)
## Calculate distances
# Distance in metres from one point to 2 other points
p1 <- sgo_points(list(-3.9369, 56.1165), epsg=4326)
lon <- c(-4.25181,-3.18827)
lat <- c(55.86424, 55.95325)
pts <- sgo_points(list(longitude=lon, latitude=lat), epsg=4326)
p1.to.pts <- sgo_distance(p1, pts, by.element = TRUE)
p1.to.pts
# Perimeter in metres of a polygon defined as a series of ordered points:
lon <- c(-6.43698696, -6.43166843, -6.42706831, -6.42102546,
-6.42248238, -6.42639092, -6.42998435, -6.43321409)
lat <- c(58.21740316, 58.21930597, 58.22014035, 58.22034112,
58.21849188, 58.21853606, 58.21824033, 58.21748949)
pol <- sgo_points(list(lon, lat), epsg=4326)
# Let's create a copy of the polygon with its coordinates shifted one position
# so that we can calculate easily the distance between vertices
coords <- sgo_coordinates(pol)
pol.shift.one <- sgo_points(rbind(coords[-1, ], coords[1, ]), epsg=pol$epsg)
perimeter <- sum(sgo_distance(pol, pol.shift.one, by.element=TRUE))
perimeter
## Area of an ordered set of points
A <- sgo_area(pol)
sprintf("%1.2f", A)
# Interpolate new vertices (here every 10 metres) if more accuracy is needed
A <- sgo_area(pol, interpolate=10)
sprintf("%1.2f", A)
```
## Benchmark
Very simple benchmark in a Intel(R) Core(TM) i7-6700HQ CPU @2.60GHz, and 16 GB of RAM.
See the code of this README file to see the R statements used.
```{r, eval=FALSE, echo=FALSE}
xmin <- -3.5
xmax <- 2.5
ymin <- 53.5
ymax <- 56.5
each_n <- function(n, reps) {
calc_coords <- function(i) {
p <- cbind(round(runif(n, xmin, xmax), 8), round(runif(n, ymin, ymax), 8))
#points(p)
#map('world', regions=('uk'), xlim=c(xmin, xmax), ylim=c(ymin, ymax))
#rect(xmin, ymin, xmax, ymax)
#2D lonlat to bng
lonlat <- sgo_points(p, epsg=4258)
start_time <- Sys.time()
bng <- sgo_transform(lonlat, to=27700)
t1 <- Sys.time()- start_time
#2D bng to lonlat
start_time <- Sys.time()
bng2ll <- sgo_transform(bng, to=4258)
t2 <- Sys.time()- start_time
#3D lonlat to bng
p.3d <- cbind(p, rep(0, n))
lonlat.3d <- sgo_points(p.3d, epsg=4937)
start_time <- Sys.time()
bng.3d <- sgo_transform(lonlat.3d, to=7405, OD=FALSE)
t3 <- Sys.time() - start_time
#3D bng to lonlat
start_time <- Sys.time()
bng.3d2ll <- sgo_transform(bng.3d, to=4979, OD=FALSE)
t4 <- Sys.time() - start_time
#3D lonlat to bng (OD=TRUE)
p.3d <- cbind(p, rep(0, n))
lonlat.3d <- sgo_points(p.3d, epsg=4937)
start_time <- Sys.time()
bng.3d <- sgo_transform(lonlat.3d, to=7405, OD=TRUE)
t5 <- Sys.time() - start_time
#3D bng to lonlat (OD=TRUE)
start_time <- Sys.time()
bng.3d2ll <- sgo_transform(bng.3d, to=4979, OD=TRUE)
t6 <- Sys.time() - start_time
out <- c(t6, t4, t2, t5, t3, t1)
}
apply(sapply(1:reps, calc_coords), 1, mean)
}
n <- c(1000, 10000, 100000, 1000000, 5000000, 10000000)
reps <- 2
result <- sapply(n, each_n, reps)
```
```{r, echo=FALSE, fig.width = 7, fig.height = 3.5, out.width="100%"}
##### PLOT RESULTS #####
n <- c(1000, 10000, 100000, 1000000, 5000000, 10000000)
result <- matrix(NA,nrow=6, ncol=6)
series <- c("BNG to LL (3D+OD)", "BNG to LL (3D)", "BNG to LL (2D)",
"LL to BNG (3D+OD)", "LL to BNG (3D)", "LL to BNG (2D)")
result[1,] <- c(0.006632924, 0.04785049, 0.3512925, 3.412248, 17.286100, 35.01287)
result[2,] <- c(0.005520940, 0.04217410, 0.2820065, 2.825136, 14.126171, 28.38339)
result[3,] <- c(0.008078456, 0.03790104, 0.2718840, 2.567024, 12.890800, 26.27742)
result[4,] <- c(0.004070640, 0.01709044, 0.1517755, 1.505867, 7.350618, 14.95144)
result[5,] <- c(0.003953457, 0.01806092, 0.1187755, 1.205132, 6.302003, 12.52924)
result[6,] <- c(0.004461050, 0.01568353, 0.1435580, 1.070004, 5.667857, 11.35794)
color <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e", "#e6ab02")
par(mar=c(5.1, 4.1, 4.1, 11.5), xpd=TRUE)
matplot(n, t(result), type="b", pch = 20, lty=1, col=color, xaxt = 'n', ylim=c(0,45), xlim=c(1, 10000000),
main="OSTN15 transformations", xlab=" pairs of coordinates", ylab="seconds")
myTicks = axTicks(1)
axis(1, at = myTicks, labels = formatC(myTicks, format = 'd', big.mark=","))
legend("topright", inset=c(-0.52,0), legend = series,
lty = c("6a"), pch=20, col = color, box.lty=0, cex=0.8)
```