8

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?

6

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)

enter image description here

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)

enter image description here

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:

enter image description here

  • Try instead of proj4string(tz.37.b) <- CRS("+init=epsg:32737") to do as i did in the post pro4string(tz.37.b)<-proj4string(a shapefile with the proj u are after) and see if it works cause i did it like i mentioned in the description and it appeared in the correct position. – gsa Oct 17 '15 at 22:29
  • So if one shp doesnt have CRS what is the first step before we apply the spTransform?give it any crs we want and then apply the spTransform and all is good? – gsa Oct 18 '15 at 12:43
  • You should definitely figure out the correct CRS and apply it first, or spTransform won't have a basis on which to re-project the object! – Simbamangu Oct 18 '15 at 13:08
4

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.