[GIS] Error While Converting UTM data to Lat-Long Coordinate

coordinate systemrsputm

I have a a dataset of Vancouver,BC where Coordinate values are projected in UTM Zone 10 . Now I can convert it using spTransform by the following piece of code and looks like the it correctly transforms the coordinates ,as I can see if I plot them in a map .

cord_transform <- function(df){
  coordinates(df) <- ~X+Y 
  proj4string(df) <- "+proj=utm +zone=10 +ellps=GRS80 +units=m +no_defs"
  df.lat.long <- spTransform(df,CRSobj=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  df.lat.long


}

However I realized I dont need a SPdataframe ,rather my purpose serves if I can transform the UTM coordinates in lat -long coordinates in the original data.table (since the data is huge I want to use data.table ) but when I follow this answer,and write a similar code for the data.table

func <- function(x,y){
  proj4string <- "+proj=utm +zone=10 +ellps=GRS80 +units= m + no_defs"
  inp <- list("x"=x,y="y")
  print(inp)
  pj <- project(inp,proj4string,inverse = TRUE)
  print(pj)
  latlon <- list(lat = pj$y,lng = pj$x)
  latlon



}
crime.data.2016.jan[,c("lat","lng"):=func(X,Y)]

I get this error :

Error in project(inp, proj4string, inverse = TRUE) : 
  unknown unit conversion id

a subset of the data looks like this :

         NEIGHBOURHOOD        X       Y
 1:       Renfrew-Collingwood 497445.9 5454283
 2: Central Business District 491580.1 5458745
 3:                           492757.5 5458792
 4:              South Cambie 491371.1 5455550
 5:                  West End 490417.0 5459169
 6:               Shaughnessy 489926.6 5454845
 7:                   Marpole 490409.0 5450869
 8:            Mount Pleasant 492545.4 5456743
 9:                  Fairview 490767.8 5456820
10:                 Kitsilano 488624.5 5457182
11:                 Killarney 497326.3 5452640
12:          Hastings-Sunrise 496763.8 5458220
13:  Kensington-Cedar Cottage 494917.4 5455246
14:                Kerrisdale 488711.0 5452722
15:       Victoria-Fraserview 495207.1 5452375
16:                Strathcona 493345.6 5458515
17:         Dunbar-Southlands 486589.5 5454941
18:                  Oakridge 490997.4 5453021
19:                    Sunset 492937.0 5451823
20:             Arbutus Ridge 488145.0 5455054
21:        Grandview-Woodland 494940.8 5458213
22:                Riley Park 492607.3 5455102
23:           West Point Grey 485649.1 5456957
24:              Stanley Park 490531.0 5460906
25:                  Musqueam 485508.0 5453015
                NEIGHBOURHOOD        X       Y

Best Answer

Your string has wrong spaces:

proj4string <- "+proj=utm +zone=10 +ellps=GRS80 +units= m + no_defs"

should be

proj4string <- "+proj=utm +zone=10 +ellps=GRS80 +units=m +no_defs"

Related Solutions

[GIS] Huge differences between results for 7-param datum transform

First, Proj4 uses what EPSG calls the "Position Vector" version of the 7 parameter method. It's possible that GeoTrans and Leica GeoOffice use the other version which EPSG called "Coordinate Frame". Both methods are equivalent, but the rotation matrices are different and the signs of the angular parameters have to be changed.

Second, thank you for providing the screen shot of the transformation definition in Leica GeoOffice. Either the Proj4 definition or this definition is defined incorrectly. The defining parameters Ellipsoid A and Ellipsoid B should be switched. Currently, this transformation converts from WGS84 to Congo60. In Proj4, note that the option is +towgs84, so it's defined as FROM Congo60 TO WGS84. To change the direction in the definition, change the signs of ALL the parameters. Also check the help for GeoOffice and see whether the SF parameter wants the parts per million version or the already converted-to-a-scale-factor value.

I don't know about GeoTrans--do you mean the NGA software? Anyway, hopefully you can get GeoOffice and Proj4 to match.

Process Vector to Raster Faster with R – Speed Optimization Techniques

I tried to "parallelize" the function rasterize using the R package parallel in this way:

split the SpatialPolygonsDataFrame object in n parts
rasterize every part separately
merge all the parts into one raster

In my computer, the parallelized rasterize function took 2.75 times less than the no-parallelized rasterize function.

Note: the code below download a polygon shapefile (~26.2 MB) from the web. You can use any SpatialPolygonDataFrame object. This is only an example.

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Figure 2: Brazil Raster plot