Proyectando objetos sp en R

31

Tengo varios shapefiles en diferentes CRS (principalmente WGS84 lat / lon) que me gustaría transformar en una proyección común (probablemente cónica Albers Equal Area, pero puedo pedir ayuda para elegir otra pregunta una vez que mi el problema se define mejor).

Pasé unos meses haciendo cosas de estadísticas espaciales en R, pero fue hace 5 años. Por mi vida, no puedo recordar cómo transformar un objeto sp (por ejemplo, SpatialPolygonsDataFrame ) de una proyección a otra.

Código de ejemplo:

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 la información de proyección adecuada, pero me gustaría transformarla en la proyección deseada. Recuerdo que hay una función con un nombre poco intuitivo para esto, pero no puedo recordar qué es.

Tenga en cuenta que no solo quiero cambiar el CRS sino cambiar las coordenadas para que coincidan ("reproyectar", "transformar", etc.).

Editar

Excluyendo AK / HI que se colocan de manera molesta en México para este shapefile:

library(taRifx.geo)
hrr.shp <- 
  subset(hrr.shp, !(grepl( "AK-" , [email protected]$HRRCITY ) |
                                     grepl( "HI-" , [email protected]$HRRCITY )) )
proj4string(hrr.shp) <- P4S.latlon
    
pregunta Ari B. Friedman 19.08.2012 - 02:53

2 respuestas

40

Puedes usar los métodos spTransform() en rgdal; usando tu ejemplo, puedes transformar el objeto a NAD83 para 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)

Paraguardarloenlanuevaproyección:

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

EDIT : O, según la sugerencia de @ Spacedman (que escribe un archivo .prj con la información de CRS):

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

Si no está seguro de a qué CRS proyectar, consulte la siguiente publicación:

Y si uno quiere definir / asignar un CRS cuando los datos no tienen uno, consulte:

respondido por el Simbamangu 19.08.2012 - 09:56
7

Desde la introducción del paquete sf (vea las viñetas sf1 , sf2 , sf3 , sf4 y una guía de migración aquí ) puede usar st_transform() para re-reproyectar sus datos vectoriales:

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 reemplazará sp en el futuro y, debido a su simplicidad y velocidad, tiene varias ventajas en comparación con sp.

    
respondido por el andrasz 24.04.2017 - 18:36

Lea otras preguntas en las etiquetas