我有这个比利时地点坐标数据框:
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对this question的回答构建的,但这些似乎是另一扇门的关键。
如果它们可以被改编/用于解决方案(可能通过以某种方式循环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)
}
我认为它可能是Belgian Lambert 1972 reference system或接近它的东西,EPSG代码31370.我基本上通过搜索search for Belgian ones的前几个EPSG代码列表来做到这一点,假设x_y坐标在该CRS中然后转换为WGS84与纬度和经度进行比较。这是我的代码;你可以看到31370的均方根误差约为490m,比我搜索的其他误差小。 (你可以显然扩大搜索范围,我尝试了一些,但这是我发现的最接近的)。请注意使用possibly
来处理当没有找到EPSG代码的转换时代码会出错的事实,因为我找不到所有可用代码的简单列表。我也不知道预期的精度是多少,因为你说lat-lon点与未知crs点完全相同。根据地区的大小和合理的采样模式,470米可能是一个关键或指示我们仍然遥远...
library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
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"))
epsg_belgium <- c(3447, 3812, 31370, 31300, 21500, 4809, 4313, 4215, 25832, 8370, 6190, 5710, 25831)
guess_crs <- function(tbl, lon_lat, x_y, epsg_codes) {
if(! requireNamespace("lwgeom"))
stop("Package 'lwgeom' must be installed to use `guess_crs`")
wgs84 <- tbl %>%
st_as_sf(coords = lon_lat, crs = 4326)
compare_crs <- function(epsg) {
tbl %>%
st_as_sf(coords = x_y, crs = epsg) %>%
st_transform(4326) %>%
st_distance(wgs84, by_element = TRUE) %>%
as.numeric %>%
`^`(2) %>%
mean %>%
sqrt
}
poss_commpare <- possibly(compare_crs, otherwise = Inf)
tibble(
crs = epsg_codes,
rms_dist = map_dbl(epsg_codes, poss_commpare)
) %>%
arrange(rms_dist)
}
guess_crs(data, c("LON", "LAT"), c("x", "y"), epsg_belgium)
#> # A tibble: 13 x 2
#> crs rms_dist
#> <dbl> <dbl>
#> 1 31370 492.
#> 2 6190 492.
#> 3 21500 533.
#> 4 31300 1101.
#> 5 3447 1695.
#> 6 3812 706664.
#> 7 25832 5466924.
#> 8 25831 5478500.
#> 9 4809 9862261.
#> 10 4215 9882639.
#> 11 8370 Inf
#> 12 5710 Inf
#> 13 4313 NA
由reprex package创建于2019-02-08(v0.2.1)