Projetando objetos sp em R

35

Tenho vários shapefiles em diferentes CRSs (principalmente WGS84 lat / lon) que gostaria de transformar em uma projeção comum (provavelmente Albers Equal Area Conic, mas posso pedir ajuda para escolher outra pergunta assim que meu problema melhorar. -definiram).

Passei alguns meses fazendo estatísticas espaciais em R, mas foi há 5 anos. Pela minha vida, não me lembro como transformar um spobjeto (por exemplo SpatialPolygonsDataFrame) de uma projeção para outra.

Código de exemplo:

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

Agora tenho uma SpatialPolygonsDataFrameinformação de projeção apropriada, mas gostaria de transformá-la na projeção desejada. Lembro-me de haver uma função com um nome pouco intuitivo para isso, mas não consigo lembrar o que é.

Observe que eu não quero apenas alterar o CRS, mas alterar as coordenadas para corresponder ("reprojetar", "transformar" etc.).

Editar

Excluindo AK / HI, que são irritantemente colocados no 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
Ari B. Friedman
fonte
Resposta anterior sobre como projetar usando o pacote proj4 aqui . Ainda não tentei isso com SpatialPolygonsDataFrame.
Simbamangu 19/08/2012
Na verdade, parece que o proj4 não funciona com objetos espaciais - mas veja a resposta abaixo.
Simbamangu 19/08/2012
2
Sempre há a exibição da tarefa espacial: cran.r-project.org/web/views/Spatial.html e minhas anotações sobre dados espaciais [plug descarado]: maths.lancs.ac.uk/~rowlings/Teaching/UseR2012
Spacedman

Respostas:

44

Você pode usar os spTransform()métodos em rgdal - usando seu exemplo, você pode transformar o objeto em 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)

não projetado

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

projetado

Para salvá-lo na nova projeção:

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

EDIT : Ou, conforme a sugestão de @ Spacedman (que grava um arquivo .prj com as informações do CRS):

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

Se não tiver certeza de qual CRS projetar, consulte a seguinte postagem:

E se alguém quiser definir / atribuir um CRS quando os dados não tiverem um, consulte:

Simbamangu
fonte
10
observe que writePolyShape NÃO grava o arquivo .prj! Você deve usar writeOGR de rgdal (e readOGR para ler shapefiles) se desejar escrever e ler o arquivo .prj para definir o CRS de seus objetos espaciais no R!
Spacedman
Muito melhor (editado de acordo) - obrigado; não tinha percebido que cria o arquivo .prj! Cheatsheet impressionante em sua página, a propósito.
Simbamangu
11
É estranho como a projeção no México afeta a aparência das inserções no Alasca e no Havaí :-).
whuber
@ whuber - hmm, sim ... alguém editou minha postagem que não tinha os mapas reais mostrando aquelas inserções inadequadas ... nunca as plotei para ver que elas estavam lá.
Simbamangu 19/08/2012
@Simbamangu Desculpe, esqueci que esse arquivo .shp incluía inadequadamente as inserções quando tentei ser útil na adição de gráficos!
Ari B. Friedman
8

Desde a introdução do pacote-sf (ter um olhar para o vinhetas sf1 , sf2 , SF3 , SF4 e um guia de migração aqui ) você pode usar st_transform()para re-reprojetar seus dados vetoriais:

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

O sf substituirá o sp no futuro e, devido à sua simplicidade e velocidade, apresenta várias vantagens em relação ao sp.

andschar
fonte