2

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)

Sam
  • 1,400
  • 13
  • 29
  • so you would like one line of code to do all this stuff instead of several lines? – Seymour Mar 28 '18 at 15:34
  • yes, i thought there would be a function to do this without the spatial object – Sam Mar 29 '18 at 07:18
  • i think the answer is `gdalUtils::gdaltransform` but i'm getting status 1 and NAs at the moment, i'll continue to explore – Sam Mar 29 '18 at 11:40

2 Answers2

1

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().

Seymour
  • 3,104
  • 2
  • 22
  • 46
  • Seymour, thanks. What I am wondering - and can change my qu - is if there is a specific function (or other) that directly works with numeric class columns and performs the mathematical transformation without faffing around with a spatial object – Sam Mar 29 '18 at 07:20
  • yours works thanks, and is just a condensed version of mine, but considering transformations are just equations, i thought there would be a function that did this to a column of numbers. Perhaps i'll just write one. – Sam Mar 29 '18 at 07:21
  • If this answer to this question please validate my answer – Seymour Mar 29 '18 at 13:47
  • it is not the answer i am looking for, sorry – Sam Mar 29 '18 at 13:50
0

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

Alec
  • 63
  • 1
  • 3