У меня есть этот фрейм данных координат бельгийских локаций:
data
# # A tibble: 11 x 4
# LON LAT x y
# <dbl> <dbl> <dbl> <dbl>
# 1 3.618942 50.68165 96227.01 152551.2
# 2 3.473466 50.55899 86074.26 138702.0
# 3 3.442395 50.69979 84369.88 154860.5
# 4 3.293505 50.68766 73127.74 153413.9
# 5 3.352688 50.68009 77876.00 153115.7
# 6 3.567919 50.52372 93229.16 134975.7
# 7 3.333666 50.54388 76183.82 137373.9
# 8 3.394737 50.61322 80806.11 145230.0
# 9 3.410566 50.53073 82041.22 135743.9
# 10 3.613925 50.59610 96907.40 143057.9
# 11 3.502580 50.74147 88860.56 159399.0
Я знаю, x и y в метрах, но я не знаю, какое соглашение используется, я получил ихкопаясь в БД общедоступных бельгийских данных, но не смог выяснить это по имеющейся информации.
Долготы и широты, представленные в таблице, не являются точными значениями (x и y являются координатамицентры областей и LAT и LON являются средними координатами выборки жителей этих регионов), но они дают хорошее представление о том, каким должен быть перевод между ними.
Какя могу выяснить, в какую систему координат x и y кодируются?
Я посмотрел на пакет sp и собрал пару функций преобразования в и изобщая система (UTM), построенная вокруг ответа @ josh-Obrien на этот вопрос , но, похоже, это ключи для другой двери.
Пожалуйста, смотрите их ниже, если они могут быть адаптированы / использованыдля соль(может быть, как-то зацикливаясь на аргументах sp::CRS)?.
Я знаю, что библиотека rgdal также выполняет преобразования координат с синтаксисом, аналогичным sp.
данные
data <- structure(list(
LON = c(3.6189419546606, 3.47346614446389, 3.44239459327957,
3.29350462630471, 3.35268808777572, 3.56791893543689, 3.33366611318681,
3.39473714826007, 3.41056562146275, 3.61392544406354, 3.50258),
LAT = c(50.6816466868977, 50.5589876530483, 50.6997902260753,
50.6876572958438, 50.6800941411327, 50.523718459466, 50.5438833669109,
50.6132227223641, 50.5307279235646, 50.5960956577015, 50.7414748843137),
x = c(96227.0052771935, 86074.2589595734, 84369.8773101478,
73127.7357132523, 77875.9986049107, 93229.1592839806, 76183.8151614011,
80806.111537044, 82041.2236842105, 96907.4010078463, 88860.5615808823),
y = c(152551.212026743, 138702.046875, 154860.466229839, 153413.886429398,
153115.726084184, 134975.700053095, 137373.913804945, 145229.97987092,
135743.853978207, 143057.883184524, 159399.019607843)),
class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -11L),
.Names = c("LON", "LAT", "x", "y"))
функции преобразования
#' add UTM coordinates (x and y in km) from WGS84 coordinates (long and lat)
#'
#' @param data a data frame
#' @param out_names names of x and y columns in output data frame
#' @param in_names names of longlat columns in input, by default searches for
#' an unambiguous set of columns with names starting with "lat" and "lon"
#' (case insensitive)
#' @param zone UTM zone, for Belgium it's 31 (the default), see
#' http://www.dmap.co.uk/utmworld.htm
#'
#' @return a data frame with 2 more columns
#' @export
#'
#' @examples
#' xy <- data.frame(ID = 1:2, longitude = c(118, 119), latitude = c(10, 50))
# add_longlat(add_utm(xy))
add_utm <- function(data, out_names = c("x","y"), in_names = NULL, zone = 31){
nms <- names(data)
if(length(temp <- intersect(out_names,nms)))
stop(paste("data already contains",paste(temp,collapse =" and ")))
if(is.null(in_names)){
lon_col <- grep("^lon",nms, ignore.case = TRUE,value = TRUE)
lat_col <- grep("^lat",nms, ignore.case = TRUE,value = TRUE)
if((n <- length(lon_col)) != 1)
stop(sprintf(
"%s columns have names starting with 'lon' (case insensitive)", n))
if((n <- length(lat_col)) != 1)
stop(sprintf(
"%s columns have names starting with 'lon' (case insensitive)", n))
in_names <- c(lon_col, lat_col)
}
new_data <- data[in_names]
sp::coordinates(new_data) <- names(new_data)
sp::proj4string(new_data) <- sp::CRS("+proj=longlat +datum=WGS84") ## for example
res <- sp::spTransform(new_data, sp::CRS(sprintf("+proj=utm +zone=%s ellps=WGS84",zone)))
res <- setNames(as.data.frame(res), out_names)
cbind(data,res)
}
#' add WGS84 coordinates (long and lat) from UTM coordinates (x and y in km)
#'
#' @param data a data frame
#' @param out_names names of longitude and latitude columns in output data frame
#' @param in_names names of longlat columns in input, by default searches for
#' an unambiguous set of columns with names starting with "x" and "y"
#' (case insensitive)
#' @param zone UTM zone, for Belgium it's 31 (the default), see
#' http://www.dmap.co.uk/utmworld.htm
#'
#' @return a data frame with 2 more columns
#' @export
#'
#' @examples
#' xy <- data.frame(ID = 1:2, longitude = c(118, 119), latitude = c(10, 50))
# add_longlat(add_utm(xy))
add_longlat <- function(data, out_names = c("LON","LAT"), in_names = NULL, zone = 31){
nms <- names(data)
if(length(temp <- intersect(out_names,nms)))
stop(paste("data already contains",paste(temp,collapse =" and ")))
if(is.null(in_names)){
lon_col <- grep("^x",nms, ignore.case = TRUE,value = TRUE)
lat_col <- grep("^y",nms, ignore.case = TRUE,value = TRUE)
if((n <- length(lon_col)) != 1)
stop(sprintf(
"%s columns have names starting with 'x' (case insensitive)", n))
if((n <- length(lat_col)) != 1)
stop(sprintf(
"%s columns have names starting with 'y' (case insensitive)", n))
in_names <- c(lon_col, lat_col)
}
new_data <- data[in_names]
sp::coordinates(new_data) <- names(new_data)
sp::proj4string(new_data) <- sp::CRS(sprintf("+proj=utm +zone=%s ellps=WGS84",zone)) ## for example
res <- sp::spTransform(new_data, sp::CRS("+proj=longlat +datum=WGS84"))
res <- setNames(as.data.frame(res), out_names)
cbind(data,res)
}