Projection d'objets sp dans R

J'ai un certain nombre de fichiers de formes dans différents CRS (principalement WGS84 lat / lon) que je voudrais transformer en une projection commune (probablement conique d'Albers Equal Area Conic, mais je pourrais demander de l'aide pour choisir une autre question une fois mon problème résolu -défini).

J'ai passé quelques mois à faire de la statistique spatiale en R, mais c'était il y a 5 ans. Pour ma vie, je ne me souviens pas comment transformer un spobjet (par exemple SpatialPolygonsDataFrame) d'une projection à une autre.

Exemple de code:

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

J'ai maintenant une SpatialPolygonsDataFrameinformation de projection appropriée, mais j'aimerais la transformer en projection souhaitée. Je me souviens qu'il existait une fonction quelque peu involontairement nommée pour cela, mais je ne me souviens plus de quoi il s'agit.

Notez que je ne veux pas simplement changer le CRS mais changer les coordonnées pour les faire correspondre ("reproject", "transformer", etc.).

modifier

En excluant AK / HI qui sont déplacés au Mexique pour ce fichier de formes:

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

Vous pouvez utiliser les spTransform()méthodes de rgdal - en utilisant votre exemple, vous pouvez transformer l'objet en NAD83 pour le Kansas (26978):

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)

Pour le sauvegarder dans la nouvelle projection:

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

EDIT : Ou, selon la suggestion de @ Spacedman (qui écrit un fichier .prj avec les informations sur le CRS):

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

Si vous ne savez pas par quel CRS projeter, reportez-vous au post suivant:

• Choisir la valeur correcte pour proj4string pour la lecture du fichier de formes dans R?

Et si on veut définir / assigner un CRS quand les données n'en ont pas, se référer à:

• Assigner CRS à un fichier de formes lorsqu'il n'en a pas, en R?

Depuis l'introduction du paquet sf (regardez les vignettes sf1 , sf2 , sf3 , sf4 et un guide de migration ici ), vous pouvez utiliser st_transform()pour reprojeter vos données vectorielles:

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 remplacera sp à l'avenir et a, en raison de sa simplicité et de sa rapidité, plusieurs avantages par rapport à sp.