Rzutowanie obiektów sp w R

35

Mam wiele plików kształtów w różnych CRS-ach (głównie WGS84 lat / lon), które chciałbym przekształcić we wspólną projekcję (prawdopodobnie stożek Albers Equal Area Conic, ale mogę poprosić o pomoc w wyborze innego pytania, gdy mój problem się poprawi) -definiowane).

Spędziłem kilka miesięcy, robiąc statystyki przestrzenne w R, ale to było 5 lat temu. Przez całe życie nie pamiętam, jak przekształcić spobiekt (np. SpatialPolygonsDataFrame) Z jednej projekcji w drugą.

Przykładowy kod:

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

Teraz mam SpatialPolygonsDataFrameodpowiednią informację o projekcji, ale chciałbym ją przekształcić w pożądaną projekcję. Przypominam sobie, że jest to funkcja o nieco nieinicjatywnej nazwie, ale nie pamiętam, co to jest.

Zauważ, że nie chcę po prostu zmieniać CRS, ale zmieniać współrzędne, aby dopasować („powtórzenie”, „transformacja” itp.).

Edytować

Z wyjątkiem AK / HI, które są denerwująco umieszczone w Meksyku dla tego pliku kształtu:

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
źródło
Poprzednia odpowiedź na temat projektowania przy użyciu pakietu proj4 tutaj . Jednak nie próbowałem tego z SpatialPolygonsDataFrame.
Simbamangu,
Właściwie wygląda na to, że proj4 nie działa z obiektami przestrzennymi - ale zobacz odpowiedź poniżej.
Simbamangu,
2
Zawsze znajduje się widok zadań przestrzennych: cran.r-project.org/web/views/Spatial.html i moje notatki na temat danych przestrzennych [bezwstydna wtyczka]: maths.lancs.ac.uk/~rowlings/Teaching/UseR2012
Spacedman

Odpowiedzi:

44

Możesz użyć tych spTransform()metod w rgdal - na twoim przykładzie możesz przekształcić obiekt w NAD83 dla 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)

nieproszony

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

przewidywane

Aby zapisać go w nowej projekcji:

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

EDYCJA : Lub, zgodnie z sugestią @ Spacedman (która zapisuje plik .prj z informacjami CRS):

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

Jeśli nie ma pewności, z którego CRS projektować, zapoznaj się z następującym postem:

A jeśli chcesz zdefiniować / przypisać CRS, gdy dane go nie mają, zapoznaj się z:

Simbamangu
źródło
10
zauważ, że writePolyShape NIE zapisuje pliku .prj! Powinieneś użyć writeOGR z rgdal (i readOGR do odczytu plików kształtów), jeśli chcesz napisać i odczytać plik .prj, aby ustawić CRS twoich obiektów przestrzennych w R!
Spacedman,
Znacznie lepiej (odpowiednio zredagowane) - dzięki; nie zdawałem sobie sprawy, że tworzy plik .prj! Nawiasem mówiąc, niesamowite ściągawki na twojej stronie.
Simbamangu,
1
Dziwne, jak projekcja w Meksyku wpływa na wygląd wypustek na Alasce i Hawajach :-).
whuber
@ whuber - hmm, tak ... ktoś edytował mój post, który nie zawierał rzeczywistych map pokazujących te raczej nieodpowiednie wstawki ... nigdy nie planowałem ich, aby się przekonać, że tam są.
Simbamangu,
@Simbamangu Niestety, zapomniałem, że ten plik .shp raczej niewłaściwie zawierał wstawki, gdy próbowałem pomóc w dodawaniu wykresów!
Ari B. Friedman,
8

Od czasu wprowadzenia pakietu sf (spójrz na winiety sf1 , sf2 , sf3 , sf4 i przewodnik po migracji tutaj ) możesz używać st_transform()do ponownej projekcji danych wektorowych:

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 zastąpi sp w przyszłości i ze względu na swoją prostotę i szybkość ma kilka zalet w porównaniu do sp.

andschar
źródło