如何在R中将英国网格参考转换为纬度和经度

3

我有一个英国英国国家网格参考系的向量:

x <- c("SK393744", "SK442746", "SK376747", "SK108191", "SP169914", "SP206935", "SK173105", "SJ993230", "SK448299", "SK112396")

我需要将这个向量转换为WGS84坐标系(纬度和经度)。

我该如何使用R进行转换?


http://www.movable-type.co.uk/scripts/latlong-gridref.html 提供了一个JavaScript实现,我可能会尝试进行转换。不确定它是否能直接与SpiderMonkey模块兼容。 - hrbrmstr
1个回答

1

试一试这些。如果它们有效,我将制作一个包含该JavaScript库中其他函数的一些更多功能的软件包(该库还有姊妹PHP和Java库,因此适合R也应该有一个)。

# takes numeric east/north generated from the os.grid.parse() function
# i shld have made it take the vector the os.grid.parse() returns but 
# we'll save that for a proper package version

os.grid.to.lat.lon <- function(E, N) {

  a <- 6377563.396
  b <- 6356256.909
  F0 <- 0.9996012717
  lat0 <- 49*pi/180
  lon0 <- -2*pi/180
  N0 <- -100000
  E0 <- 400000
  e2 <- 1 - (b^2)/(a^2)
  n <- (a-b)/(a+b)
  n2 <- n^2
  n3 <- n^3

  lat <- lat0
  M <- 0

  repeat {

    lat <- (N-N0-M)/(a*F0) + lat

    Ma <- (1 + n + (5/4)*n2 + (5/4)*n3) * (lat-lat0)
    Mb <- (3*n + 3*n*n + (21/8)*n3) * sin(lat-lat0) * cos(lat+lat0)
    Mc <- ((15/8)*n2 + (15/8)*n3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
    Md <- (35/24)*n3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))
    M <- b * F0 * (Ma - Mb + Mc - Md)

    if (N-N0-M < 0.00001) { break }

  }

  cosLat <- cos(lat)
  sinLat <- sin(lat)

  nu <- a*F0/sqrt(1-e2*sinLat*sinLat)
  rho <- a*F0*(1-e2)/((1-e2*sinLat*sinLat)^1.5)

  eta2 <- nu/rho-1

  tanLat <- tan(lat)
  tan2lat <- tanLat*tanLat
  tan4lat <- tan2lat*tan2lat
  tan6lat <- tan4lat*tan2lat

  secLat <- 1/cosLat
  nu3 <- nu*nu*nu
  nu5 <- nu3*nu*nu
  nu7 <- nu5*nu*nu

  VII <- tanLat/(2*rho*nu)
  VIII <- tanLat/(24*rho*nu3)*(5+3*tan2lat+eta2-9*tan2lat*eta2)
  IX <- tanLat/(720*rho*nu5)*(61+90*tan2lat+45*tan4lat)
  X <- secLat/nu
  XI <- secLat/(6*nu3)*(nu/rho+2*tan2lat)
  XII <- secLat/(120*nu5)*(5+28*tan2lat+24*tan4lat)
  XIIA <- secLat/(5040*nu7)*(61+662*tan2lat+1320*tan4lat+720*tan6lat)

  dE <- (E-E0)
  dE2 <- dE*dE
  dE3 <- dE2*dE
  dE4 <- dE2*dE2
  dE5 <- dE3*dE2
  dE6 <- dE4*dE2
  dE7 <- dE5*dE2

  lon <- lon0 + X*dE - XI*dE3 + XII*dE5 - XIIA*dE7
  lat <- lat - VII*dE2 + VIII*dE4 - IX*dE6

  lat <- lat * 180/pi
  lon <- lon * 180/pi

  return(c(lat, lon))

}


# takes a string OS reference and returns an E/N vector

os.grid.parse <- function(grid.ref) {

  grid.ref <- toupper(grid.ref)

  # get numeric values of letter references, mapping A->0, B->1, C->2, etc:
  l1 <- as.numeric(charToRaw(substr(grid.ref,1,1))) - 65
  l2 <- as.numeric(charToRaw(substr(grid.ref,2,2))) - 65

  # shuffle down letters after 'I' since 'I' is not used in grid:
  if (l1 > 7) l1 <- l1 - 1
  if (l2 > 7) l2 <- l2 - 1

  # convert grid letters into 100km-square indexes from false origin - grid square SV

  e <- ((l1-2) %% 5) * 5 + (l2 %% 5)
  n <- (19 - floor(l1/5) *5 ) - floor(l2/5)

  if (e<0 || e>6 || n<0 || n>12) { return(c(NA,NA)) }

  # skip grid letters to get numeric part of ref, stripping any spaces:

  ref.num <- gsub(" ", "", substr(grid.ref, 3, nchar(grid.ref)))
  ref.mid <- floor(nchar(ref.num) / 2)
  ref.len <- nchar(ref.num)

  if (ref.len >= 10) { return(c(NA,NA)) }

  e <- paste(e, substr(ref.num, 0, ref.mid), sep="", collapse="")
  n <- paste(n, substr(ref.num, ref.mid+1, ref.len), sep="", collapse="")

  nrep <- 5 - match(ref.len, c(0,2,4,6,8))

  e <- as.numeric(paste(e, "5", rep("0", nrep), sep="", collapse=""))
  n <- as.numeric(paste(n, "5", rep("0", nrep), sep="", collapse=""))

  return(c(e,n))

}

proj4包应该能够帮助你获取最终的WGS84数值(答案中的代码可以将你转换为OSGB36纬度/经度)。 - hrbrmstr
谢谢你的回答!我尝试过:os.grid.parse("SK393744") [1] 439350 374450 但是我期望的结果是:439300 374400... 你有什么想法吗? - Claudia
谢谢。我希望能得到一些验证,确认它有点偏差(我本来就觉得可能有问题,但我从不使用网格参考,所以很难判断)。我会将JavaScript实现与Java和PHP的代码进行比较,并进行调整和优化。这绝对需要成为一个包。 - hrbrmstr
还有一些不太对劲的地方。如果我使用网格参考“ND262549”,应该得到58.476061 -3.2670900。但是使用你的函数,我得到的是58.476881 -3.264716。在打包之前最好修复这些问题吧 ;) - Claudia
确实。我需要在这上面引用所有的库。敬请关注。联系方式:rudis dot net 的 bob,用于离线信息/评论/见解。对英国和其他地方表示没有冒犯之意,但这个网格事情真的很复杂 :-) - hrbrmstr

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接