43 votos

La proyección de sp objetos en R

Tengo una serie de archivos en diferentes Sir (en su mayoría WGS84 (lat/lon) que me gustaría transformar en común de proyección (probablemente proyección Cónica Equivalente de Albers, pero me pueden pedir ayuda en la elección en otra pregunta una vez que mi problema se presenta mejor definido).

Me pasé un par de meses haciendo espacial estadísticas cosas en R, pero fue hace 5 años. Para la vida de mí, no puedo recordar cómo transformar un sp objeto (por ejemplo SpatialPolygonsDataFrame) a partir de una proyección a otra.

Ejemplo de código:

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

Ahora tengo un SpatialPolygonsDataFrame con proyección adecuada de la información, pero me gustaría transformar a la proyección deseada. Recuerdo que hay un poco de unintuitively nombre de la función para esto, pero no puedo recordar lo que es.

Tenga en cuenta que no quiero solo para cambiar el CRS, pero a cambio de las coordenadas de partido ("reproyectar", "transformar", etc.).

Editar

Excluyendo AK/HI que son molestas y se coloca en México para este shapefile:

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

51voto

Usted puede utilizar el spTransform() métodos en rgdal - con su ejemplo, usted puede transformar el objeto a NAD83 de 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)

unprojected

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

projected

Para guardarlo en la nueva proyección:

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

EDITAR: O, como por @Spacedman la sugerencia (que escribe .prj archivo con la CRS info):

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

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X