Sam
Sam

Reputation: 1482

R transforming coordinates inside the data.frame

Converting coordinates in R is straightforward with spTransform etc, but is there a way to bypass the Spatial object and convert directly in the dataframe? To just run the transformation equation on 2 columns? E.g. from latlon to British National Grid as new columns:

# 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)

However i want to go from the original dummy data ('toconv') straight to the final output ('final') without faffing around, is is possible in one function? (e.g. a function containing the Helmert Transformation or OSTN02 Transformation)

Upvotes: 2

Views: 4459

Answers (2)

Alec
Alec

Reputation: 65

Here's a small function that does this with data.table, sf and 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

Created on 2021-01-25 by the reprex package (v0.3.0)

Also note this related question: Converting latitude and longitude points to UTM

Upvotes: 0

Seymour
Seymour

Reputation: 3264

I messed around quite a lot for answering your question.

I understood that you are looking for a straightforward function to transform coordinates between different projections.

The desired final output is:

  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

I tried several approaches using both packages proj4 and the one you named in your comment. Unfortunately, the results were different than those obtained using the brilliant sp packages created by Dr. Pebesma.

Therefore, my final solution to your question is to create a function named help_sam that make straightforward for you to change the coordinate reference system of data.frame structured as 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

If you want to change the CRS of the final projection you just have to set a different espg value for parameter dst.proj in function help_sam().

Upvotes: 2

Related Questions