3 ms·

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 Merc

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 3y ago
Yes, `lm()` linearly interpolates. Thanks for pointing out my error.
- codeflo 3y ago
  To be fair, I’m not sure this will actually make a significant difference. :) Even 200m wouldn’t be a huge error.
  - bluenose69 3y 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]"))

bluenose69 3y ago

Yes, `lm()` linearly interpolates. Thanks for pointing out my error.

codeflo 3y ago
To be fair, I’m not sure this will actually make a significant difference. :) Even 200m wouldn’t be a huge error.
- bluenose69 3y 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]"))