r转换data.frame内部的坐标

Question

在r中转换坐标与Sptransform等很简单，但是有没有办法绕过空间对象并直接在数据框架中转换？只是在2列上运行转换方程？例如。从拉特隆到英国国家电网，作为新专栏：

# current way using a spatial object
require(raster)
require(rgdal)
# define BNG and latlon
BNG <- CRS("+init=epsg:27700")
LL <- CRS("+init=epsg:4326")
# dummy data
toconv <- data.frame(id=c("a","b","c"), lat=c(54.530776,54.551913,54.455268), lon=c(-2.6006958,-2.4084351,-2.4688599))
# promote to spatial points data frame and define CRS of points
coordinates(toconv) = ~lon + lat
crs(toconv) <- LL
# current LL coordinates as columns in the SPDF
toconv$Xlon <- coordinates(toconv)[,1]
toconv$Ylat <- coordinates(toconv)[,2]
# transform to BNG
conv <- spTransform(toconv, crs(BNG))
# rename the coords from original name to new wanted name
colnames(conv@coords) <- c("Xbng","Ybng")
# extract as data frame, new coords with new name are new columns. 
final <- as.data.frame(conv)

但是，我想直接从原始的虚拟数据('toconv'(直接转到最终输出('final'(而不会绕过，在一个函数中是否可以使用吗？(例如，包含Helmert转换或OSTN02转换的功能(

Answer 1

我在回答您的问题时搞砸了很多。

我知道您正在寻找一个直接的功能来在不同的预测之间转换坐标。

所需的最终输出是：

  id      Xlon     Ylat     Xbng     Ybng
1  a -2.600696 54.53078 361224.2 515221.4
2  b -2.408435 54.55191 373679.8 517484.1
3  c -2.468860 54.45527 369699.8 506754.6

我尝试使用两种软件包proj4和您在评论中命名的方法尝试了几种方法。不幸的是，结果与使用Pebesma博士创建的出色sp软件包获得的结果不同。

因此，我对问题的最终解决方案是创建一个名为 help_sam的函数，该函数使您直接更改data的坐标参考系统。帧构造为toconv。

BNG <- CRS("+init=epsg:27700")
LL <- CRS("+init=epsg:4326")
toconv <- data.frame(id=c("a","b","c"), lat=c(54.530776,54.551913,54.455268), lon=c(-2.6006958,-2.4084351,-2.4688599))

help_sam = function(data,
                    src.proj = CRS("+init=epsg:4326"),
                    dst.proj = CRS("+init=epsg:27700")) {
        require(sp)
        as.data.frame(
                spTransform(
                        SpatialPointsDataFrame(
                                coords = data.frame(Xbng = toconv$lon,
                                                    Ybng = toconv$lat),
                                data = data.frame(id = toconv$id,
                                                  Xlon = toconv$lon,
                                                  Ylat = toconv$lat),
                                proj4string = src.proj), dst.proj))
}
final <- help_sam(data = toconv)
print(final)
  id      Xlon     Ylat     Xbng     Ybng
1  a -2.600696 54.53078 361224.2 515221.4
2  b -2.408435 54.55191 373679.8 517484.1
3  c -2.468860 54.45527 369699.8 506754.6

如果要更改最终投影的CRS，则只需为函数help_sam()中的参数dst.proj设置其他espg值。

Answer 2

这是一个使用data.table，sf和rgdal的小函数。

## Packages
library(sf)
library(data.table)
library(rgdal)

## Data
# Example data from question
# Note: using a data.table instead of a data.frame here
#       so we can use the function below to add new projected columns
toconv <-
    data.table(
        id = c("a", "b", "c"),
        lat = c(54.530776, 54.551913, 54.455268),
        lon = c(-2.6006958, -2.4084351, -2.4688599)
    )

## Function
# Project coordinates inside a data.table within converting to a spatial object
# Uses rgdal::project a matrix of geographical positions to projected coordinates
# Expectated order: X coordinate, Y coordinate
project_coords <- function(DT, coords, projection, projcoords) {
    DT[, (projcoords) :=
            data.table::as.data.table(
                rgdal::project(
                    as.matrix(.SD, ncol = 2),
                    projection)
            ),
         .SDcols = coords][]
}

# Useful checks to add:
#   Are the columns named in coords numeric?
#   Is the projection a character?
#   Is it a valid projection?
#   Are the coord1inates in the right order?
## Usage
# Setup output CRS (using the sf function st_crs and returning the WKT representation)
projection <- st_crs(27700)$wkt

# Project coordinates
project_coords(
    DT = toconv,
    coords = c('lon', 'lat'),
    projection = projection,
    projcoords = c('Xbng', 'Ybng')
)
#>    id      lat       lon     Xbng     Ybng
#> 1:  a 54.53078 -2.600696 361131.2 515235.4
#> 2:  b 54.55191 -2.408435 373585.3 517497.9
#> 3:  c 54.45527 -2.468860 369605.8 506769.7

^{由Reprex软件包(v0.3.0(在2021-01-25创建}

还请注意以下问题：将纬度和经度转换为UTM