Memproyeksikan objek sp dalam R

35

Saya memiliki sejumlah shapefile di CRS yang berbeda (kebanyakan WGS84 lat / lon) yang ingin saya ubah menjadi proyeksi umum (kemungkinan Albers Equal Area Conic, tetapi saya dapat meminta bantuan dalam memilih pertanyaan lain setelah masalah saya membaik -defined).

Saya menghabiskan beberapa bulan melakukan hal statistik spasial di R, tapi itu 5 tahun yang lalu. Untuk kehidupan saya, saya tidak ingat bagaimana mengubah suatu spobjek (misalnya SpatialPolygonsDataFrame) dari satu proyeksi ke proyeksi lainnya.

Kode contoh:

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

Sekarang saya memiliki SpatialPolygonsDataFrameinformasi proyeksi yang sesuai, tetapi saya ingin mengubahnya menjadi proyeksi yang diinginkan. Saya ingat ada fungsi yang agak tidak sengaja dinamai untuk ini, tapi saya tidak ingat apa itu.

Perhatikan bahwa saya tidak ingin hanya mengubah CRS tetapi untuk mengubah koordinat yang cocok ("proyeksi ulang", "transformasi", dll.).

Edit

Tidak termasuk AK / HI yang ditempatkan secara mengganggu di Meksiko untuk shapefile ini:

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
r coordinate-system Ari B. Friedman
sumber
Jawaban sebelumnya tentang memproyeksikan menggunakan paket proj4 di sini . Belum mencoba ini dengan SpatialPolygonsDataFrame.
Simbamangu
Sebenarnya sepertinya proj4 tidak bekerja dengan objek Spasial - tetapi lihat jawaban di bawah ini.
Simbamangu
2
Selalu ada Tampilan Tugas Spasial: cran.r-project.org/web/views/Spatial.html dan catatan saya tentang Data Tata Ruang [plug shameless]: maths.lancs.ac.uk/ ~ paylings
Teaching

Jawaban:

44

Anda dapat menggunakan spTransform()metode dalam rgdal - menggunakan contoh Anda, Anda dapat mengubah objek ke NAD83 untuk 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)

tidak diproyeksikan

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

diproyeksikan

Untuk menyimpannya dalam proyeksi baru:

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

EDIT : Atau, sesuai saran @ Spacedman (yang menulis file .prj dengan info CRS):

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

Jika seseorang tidak yakin dari CRS mana untuk memproyeksikan, lihat posting berikut:

Dan jika seseorang ingin mendefinisikan / menetapkan CRS ketika data tidak memilikinya, lihat:

Simbamangu
sumber
10
perhatikan bahwa writePolyShape TIDAK menulis file .prj! Anda harus menggunakan writeOGR dari rgdal (dan menggunakan readOGR untuk membaca shapefile) jika Anda ingin menulis dan membaca file .prj untuk mengatur CRS objek spasial Anda di R!
Spacedman
Jauh lebih baik (diedit sesuai) - terima kasih; belum sadar itu membuat file .prj! Omong kosong di halaman Anda, ngomong-ngomong.
Simbamangu
1
Sungguh aneh bagaimana proyeksi di Meksiko mempengaruhi penampilan inset Alaska dan Hawaii :-).
whuber
@whuber - hmm, ya ... seseorang mengedit posting saya yang tidak memiliki peta aktual di dalamnya yang menunjukkan insets yang agak tidak pantas ... tidak pernah merencanakannya sendiri untuk melihat bahwa mereka ada di sana.
Simbamangu
@Simbamangu Maaf, lupa bahwa file .shp ini agak tidak benar memasukkan insets ketika saya berusaha membantu dalam menambahkan grafik!
Ari B. Friedman
8

Sejak diperkenalkannya paket sf (lihat sketsa sf1 , sf2 , sf3 , sf4 dan panduan migrasi di sini ), Anda dapat menggunakan st_transform()untuk memproyeksi ulang data vektor Anda:

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 akan menggantikan sp di masa depan dan memiliki kesederhanaan dan kecepatan karena beberapa keunggulan dibandingkan sp.

andschar
sumber