Alternativen von spTransform?

Nehmen wir an, wir haben ein Shapefile mit einer bestimmten Projektion.

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

Wir haben auch die Flughäfen von Spanien als Punkte in einer anderen Projektion.

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

Wenn wir die Punkte auf dem Shapefile Spaniens platzieren möchten, müssen wir die Koordinaten so anordnen, dass sie gleich sind. Normalerweise wird es so gemacht:

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

Aber ist das auch so?

proj4string(a)<-proj4string(s)

Wenn ja, warum ist dies dann nicht die Standardmethode, da es einfacher ist und auf spTransform zurückgreifen muss?

Durch das Zuweisen eines anderen CRS wird die Projektion der zugrunde liegenden räumlichen Daten nicht geändert. Das CRS ist ein interner Teil des räumlichen Objekts, der R mitteilt, wie die räumlichen Koordinaten zu interpretieren sind.

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

Transformieren Sie nun die Form in eine benachbarte 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]

Die Koordinaten der Geometrie haben sich korrekt geändert, um sie an die benachbarte UTM-Zone anzupassen. Was ist, wenn wir das CRS einfach ohne Transformation neu zuweisen?

# 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

Wir wurden gewarnt ... wie sehen die Koordinaten aus?

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]

Die Zahlen (Koordinaten) sind die ursprünglichen aus der UTM-Zone 36, aber die Form wird jetzt der falschen UTM-Zone zugeordnet und an der falschen Stelle angezeigt.

Bearbeiten

Verwenden Sie genau die vom OP vorgeschlagene Methode:

# 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 ist das gleiche wie das ursprüngliche utm36-Objekt, aber proj4string ist jetzt utm37. (Interessanterweise gab es diesmal keine Warnung.) Hinweis: Das Ergebnis ist genau das gleiche wie das Zuweisen des CRS mit dem EPSG-Code. Zu testen:

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

Lassen Sie uns anhand einer tatsächlichen UTM37-Schicht der Regionen Tansanias überprüfen, wie die tatsächlichen Formen aussehen.

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

Sie reihen sich nicht aneinander, also funktioniert diese Methode nicht! Was ist mit dem transformierten Objekt?

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

Das transformierte Objekt wird an der richtigen Stelle angezeigt (dunkelgrüner Hintergrund).

Hinweis: Wenn Sie zwischen UTM-Zonen mit demselben Datum (hier WGS84) transformieren, sind die Unterschiede wie in meinem Beispiel dramatisch. Wenn Sie jedoch zwischen verschiedenen Datumsangaben transformieren, sind die Unterschiede möglicherweise viel subtiler. Im Folgenden finden Sie beispielsweise zwei Formen aus Sansibar, beide UTM 37, aber eine ist WGS84 und die andere ist ARC1960:

Nein, diese beiden Formen sind nicht gleich, und dies ist beabsichtigt. Der Grund dafür ist, dass Sie in einigen Fällen einem Objekt, das kein CRS hat, ein CRS zuweisen oder ein CRS ändern möchten, ohne die Koordinaten zu ändern, und in anderen Fällen möchten Sie Koordinaten von einem CRS in ein anderes CRS konvertieren oder transformieren. Das erste (Zuweisen, Ersetzen) wird erhalten durch

proj4string(a)<-proj4string(s)

und die zweite (konvertieren, transformieren) durch

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

Da Sie nicht der erste sind, der denkt, dass der erste Befehl tatsächlich das tut, was der zweite Befehl tut, spgibt er die folgende Warnung aus, wenn Sie das erste Formular versuchen, aber abereits ein anderes CRS haben:

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

(was nicht mehr ganz richtig ist, da spTransformeine Funktion spdarin Methoden aufruft rgdal).