Alternatives of spTransform?

Question

Let's say that we have a shapefile with a certain projection.

s<-readOGR(dsn=".",layer="Spain")

We also have the airports of Spain as points in another projection.

a<-readOGR(dsn=".",layer="airports")

If we aim to place the points onto the Spain's shapefile we have to arrange the coordinates so they are the same. Usually, it is done like this:

 a<-spTransform(a,CRS(proj4string(s))

But is it the same with that?

proj4string(a)<-proj4string(s)

If yes, then why isn't the standard way of doing it since it is simpler and have to resort to the spTransform?

Answer 1

Assigning a different CRS does not change the projection of the underlying spatial data - the CRS is an internal part of the spatial object that tells R how to interpret the spatial coordinates.

library(rgdal)

# Load Tanzania in UTM 36
tz.36 <- readOGR(dsn = ".", layer = "tz_36")
summary(tz.36) # Show the bounding coordinates:

Coordinates:
         min     max
x   94000.58 1315978
y 8699697.87 9889701
Is projected: TRUE 
proj4string :
[+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
....

Now transform the shape to a neighbouring UTM zone:

tz.37 <- spTransform(tz.36, CRS("+init=epsg:32737"))
summary(tz.37)

Coordinates:
        min       max
x -576091.7  657248.1
y 8700995.0 9888925.5
Is projected: TRUE 
proj4string :
[+init=epsg:32737 +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]

The geometry's coordinates have changed - correctly - to fit neighbouring UTM zone. What if we simply reassign the CRS without transformation?

# Make a copy of the original object.
tz.37.b <- tz.36 
# Assign the CRS
proj4string(tz.37.b) <- CRS("+init=epsg:32737") 

Warning message:
In `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
  A new CRS was assigned to an object with an existing CRS:
+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
without reprojecting.
For reprojection, use function spTransform in package rgdal

We've been warned ... so what do the coordinates look like?

summary(tz.37.b)

    Coordinates:
         min     max
x   94000.58 1315978
y 8699697.87 9889701
Is projected: TRUE 
proj4string :
[+init=epsg:32737 +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]

The numbers (coordinates) are the original ones from UTM zone 36, but the shape will be mapped to the wrong UTM zone now, and appear in the wrong place.

Edit

Using the exact method suggested by the OP:

# Make a new copy of the original UTM 36 shape:
tz.37.c <- tz.36
# Now assign the proj4string using OP's suggestion:
proj4string(tz.37.c) <- proj4string(tz.37)
summary(tz.37.c)

Coordinates:
         min     max
x   94000.58 1315978
y 8699697.87 9889701
Is projected: TRUE 
proj4string :
[+init=epsg:32737 +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]

Extent is the same as original utm36 object, but proj4string is now utm37. (Interestingly, there was no warning this time.) Note: the result is EXACTLY the same as having assigned the CRS with the EPSG code. To test:

identical(tz.37.b, tz.37.c)
TRUE

Let's check what the actual shapes look like, using an actual UTM37 layer of the regions of Tanzania.

# Load the region shape
tz.regions.37 <- readOGR(".", "tz_regions_37_simp")
# Plot it
plot(tz.regions.37, lwd = 1, border = 'red')
# Now add the CRS-assigned (not transformed!) object:
plot(tz.37.c, add = T, border = 'blue', lwd = 2)

They don't line up, so this method isn't working! What about the transformed object?

plot(tz.37, add = T, border = 'darkgreen', lwd = 2)

The transformed object appears in the right place (dark green background).

N.B. If you are transforming between UTM zones but the same datum (WGS84 here) the differences will be dramatic, as in my example - but if you are transforming between different datums the differences may be much more subtle. For example, below is are two shapes from Zanzibar, both UTM 37, but one is WGS84 and the other is ARC1960:

Answer 2

No: these two forms are not the same, and this is by design. The reason is that in some cases you want to assign a CRS to an object which does not have a CRS or change a CRS without changing the coordinates, and in other cases you want to convert or transform coordinates from one CRS into another CRS. The first (assigning, replacing) is obtained by

proj4string(a)<-proj4string(s)

and the second (convert, transform) by

a<-spTransform(a,CRS(proj4string(s))

Becauses you are not the first to think that the first command actually does what the second command does, sp gives the following warning when you try the first form but a already has a different CRS:

Warning message:
In ReplProj4string(obj, CRS(value)) :
  A new CRS was assigned to an object with an existing CRS:
+init=epsg:28992 +proj=sterea [etc]
without reprojecting.
For reprojection, use function spTransform in package rgdal

(which is no longer entirely correct, since spTransform is a function in sp that calls methods in rgdal).