Projizieren von sp-Objekten in R

Ich habe eine Reihe von Shapefiles in verschiedenen CRSs (meistens WGS84 lat / lon), die ich in eine gemeinsame Projektion umwandeln möchte (wahrscheinlich Albers Equal Area Conic), aber ich kann Sie um Hilfe bei der Auswahl einer anderen Frage bitten, sobald sich mein Problem bessert -definiert).

Ich habe ein paar Monate damit verbracht, Raumstatistiken in R zu erstellen, aber das war vor 5 Jahren. Ich kann mich nicht erinnern, wie ich ein spObjekt (z. B. SpatialPolygonsDataFrame) von einer Projektion in eine andere verwandeln kann .

Beispielcode:

P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
hrr.shp <- readShapePoly("HRR_Bdry"), verbose=TRUE, proj4string=P4S.latlon) 
# Shapefile available at 
#   http://www.dartmouthatlas.org/downloads/geography/hrr_bdry.zip 
#   but you must rename all the filenames to have the same 
#   capitalization for it to work in R

Jetzt habe ich eine SpatialPolygonsDataFramemit entsprechender Projektionsinformation, möchte diese aber in die gewünschte Projektion umwandeln. Ich erinnere mich, dass es dafür eine etwas unintuitiv benannte Funktion gibt, aber ich kann mich nicht erinnern, was es ist.

Beachten Sie, dass ich nicht nur das CRS ändern möchte, sondern die Koordinaten so ändern möchte, dass sie übereinstimmen ("Umprojektieren", "Transformieren" usw.).

Bearbeiten

Mit Ausnahme von AK / HI, die für dieses Shapefile ärgerlich in Mexiko platziert sind:

library(taRifx.geo)
hrr.shp <- 
  subset(hrr.shp, !(grepl( "AK-" , hrr.shp@data$HRRCITY ) |
                                     grepl( "HI-" , hrr.shp@data$HRRCITY )) )
proj4string(hrr.shp) <- P4S.latlon

Sie können die spTransform()Methoden in rgdal verwenden - anhand Ihres Beispiels können Sie das Objekt in NAD83 für Kansas (26978) umwandeln:

library(rgdal)
library(maptools)

P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
hrr.shp <- readShapePoly("HRR_Bdry", verbose=TRUE, proj4string=P4S.latlon)
plot(hrr.shp)

hrr.shp.2 <- spTransform(hrr.shp, CRS("+init=epsg:26978"))
plot(hrr.shp.2)

So speichern Sie es in der neuen Projektion:

writePolyShape(hrr.shp.2, "HRR_Bdry_NAD83")

BEARBEITEN : Oder laut dem Vorschlag von @ Spacedman (der eine .prj-Datei mit den CRS-Informationen schreibt):

writeOGR(hrr.shp.2, dsn = getwd(), layer = "HRR_Bdry_NAD83", driver="ESRI Shapefile")

Wenn Sie nicht sicher sind, von welchem ​​CRS Sie projizieren sollen, lesen Sie den folgenden Beitrag:

• Wählen Sie den richtigen Wert für proj4string für das Lesen von Shapefiles in R?

Und wenn man ein CRS definieren / zuweisen möchte, wenn Daten kein CRS haben, siehe:

• Zuweisen von CRS zu Shapefile in R?

Seit der Einführung des sf-Pakets (hat einen Blick auf den Vignetten SF1 , SF2 , SF3 , SF4 und einen Migrationsleitfaden hier ) Sie können st_transform()für den Wieder Reprojizierens Ihre Vektordaten:

require(sf)

hrr_sf = st_read('HRR_Bdry.shp', stringsAsFactors = FALSE,
    crs = 4326) # has +proj=longlat +datum=WGS84
plot(hrr_sf)

hrr_sf2 = st_transform(hrr_sf, "+init=epsg:26978") # 1st option sp::CRS() not working/ needed
hrr_sf2 = st_transform(hrr_sf, 26978) # 2nd option - EPSG code as an integer
plot(hrr_sf2)

# don't think about doing this:
hrr_sf3 = st_read('HRR_Bdry.shp', stringsAsFactors = FALSE,
    crs = 26978)

# Output layer
st_write(hrr_sf2, dsn = getwd(), layer = "HRR_Bdry_NAD83", driver = "ESRI Shapefile")

sf wird in zukunft sp ersetzen und hat aufgrund seiner einfachheit und geschwindigkeit imho mehrere vorteile gegenüber sp.