4 ms·

I get a bit over 200 metres from the following R code. (Edit: I forgot how to format code for HM, sorry.) library(oce) lat <- c(57+18.280/60, 56+54.15

I get a bit over 200 metres from the following R code. (Edit: I forgot how to format code for HM, sorry.) library(oce) lat <- c(57+18.280/60, 56+54.152/60) lon <- -c(3+55.948/60, 4+56.664/60) m <- lm(lat ~ lon) # Find points along a great-circle route gc <- oce::geodGc(lon, lat, 0.001) # 'scale' is distance associated with 1 degree latitude shift scale <- geodDist(mean(lon), mean(lat), mean(lon), mean(lat)+1) # 'e' is latitude shift from great-circle route to 'straight' line route e <- scale\*sapply(seq_along(gc$longitude), function(i) gc$latitude[i] - predict(m, data.frame(lon=gc$longitude[i]))) plot(gc$longitude, e, type="l", xlab="Longitude [degE]", ylab="Distance error [km]", col=2) max(e)

codeflo 4y ago
I’m not familiar with R, but the lm(…) line looks like it’s linearly interpolating the latitude and longitude. However, that’s not what straight lines in a Mercator projection are. You might be confusing it with the cylindrical projection. The Mercator projection stretches things away from the equator vertically to compensate the horizontal stretching, in order to preserve local shapes. This makes it a bit better for these kinds of local measurements than one might think at first glance.
- bluenose69 4y ago
  Yes, `lm()` linearly interpolates. Thanks for pointing out my error.
  - codeflo 4y ago
    To be fair, I’m not sure this will actually make a significant difference. :) Even 200m wouldn’t be a huge error.
    - bluenose69 4y ago
      I tried a computation that I think may be correct. It gives 319 metres of northing separation between rhumb and creat-circle paths. In case it's of any interest, the R code is below. # Northing distance between rhumb and great-circle paths. library(oce) lat <- c(57+18.280/60, 56+54.152/60) lon <- -c(3+55.948/60, 4+56.664/60) # Great-circle path, first in longitude-latitude space, and then # (for Mercator projection) in easting-northing space. gc <- oce::geodGc(lon, lat, 0.01) # same results for 3rd arg any value < 0.01 gcXY <- lonlat2map(gc$longitude, gc$latitude, projection="+proj=merc") # Create function, f(easting), that computes northing value on # a rhumb line. We need this to make the two lines share easting # values. p <- oce::lonlat2map(lon, lat, "+proj=merc") m <- lm(y ~ x, data=p) X <- seq(p$x[1], p$x[2], length.out=length(gc$longitude)) Y <- predict(m, data.frame(x=X)) f <- approxfun(X, Y) # interpolating function # Compute the northing error. northingError <- f(gcXY$x) - gcXY$y maxNorthingError <- round(max(abs(northingError)), 1) message("maxNorthingError=", maxNorthingError, " [m]") plot(gcXY$x, northingError, type="l", xlab="Easting [m]", ylab="Northing error [m]") mtext(paste("max northing error ", maxNorthingError, " [m]"))

codeflo 4y ago

I’m not familiar with R, but the lm(…) line looks like it’s linearly interpolating the latitude and longitude. However, that’s not what straight lines in a Mercator projection are. You might be confusing it with the cylindrical projection. The Mercator projection stretches things away from the equator vertically to compensate the horizontal stretching, in order to preserve local shapes. This makes it a bit better for these kinds of local measurements than one might think at first glance.

bluenose69 4y ago
Yes, `lm()` linearly interpolates. Thanks for pointing out my error.
- codeflo 4y ago
  To be fair, I’m not sure this will actually make a significant difference. :) Even 200m wouldn’t be a huge error.
  - bluenose69 4y ago
    I tried a computation that I think may be correct. It gives 319 metres of northing separation between rhumb and creat-circle paths. In case it's of any interest, the R code is below. # Northing distance between rhumb and great-circle paths. library(oce) lat <- c(57+18.280/60, 56+54.152/60) lon <- -c(3+55.948/60, 4+56.664/60) # Great-circle path, first in longitude-latitude space, and then # (for Mercator projection) in easting-northing space. gc <- oce::geodGc(lon, lat, 0.01) # same results for 3rd arg any value < 0.01 gcXY <- lonlat2map(gc$longitude, gc$latitude, projection="+proj=merc") # Create function, f(easting), that computes northing value on # a rhumb line. We need this to make the two lines share easting # values. p <- oce::lonlat2map(lon, lat, "+proj=merc") m <- lm(y ~ x, data=p) X <- seq(p$x[1], p$x[2], length.out=length(gc$longitude)) Y <- predict(m, data.frame(x=X)) f <- approxfun(X, Y) # interpolating function # Compute the northing error. northingError <- f(gcXY$x) - gcXY$y maxNorthingError <- round(max(abs(northingError)), 1) message("maxNorthingError=", maxNorthingError, " [m]") plot(gcXY$x, northingError, type="l", xlab="Easting [m]", ylab="Northing error [m]") mtext(paste("max northing error ", maxNorthingError, " [m]"))