Проектирование sp объектов в R

35

У меня есть несколько шейп-файлов в разных CRS (в основном WGS84 lat / lon), которые я хотел бы преобразовать в общую проекцию (вероятно, коническую равную площадь Альберса, но я могу попросить помощи при выборе другого вопроса, как только моя проблема станет лучше) -определенный).

Я провел несколько месяцев, занимаясь пространственной статистикой в ​​R, но это было 5 лет назад. Что касается меня, я не могу вспомнить, как преобразовать spобъект (например SpatialPolygonsDataFrame) из одной проекции в другую.

Пример кода:

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

Теперь у меня есть SpatialPolygonsDataFrameсоответствующая информация о проекции, но я хотел бы преобразовать ее в желаемую проекцию. Я вспоминаю, что для этого была функция с некоторым неинтуитивным названием, но я не могу вспомнить, что это такое.

Обратите внимание, что я хочу не просто изменить CRS, а изменить координаты для соответствия («перепроектировать», «преобразовать» и т. Д.).

редактировать

За исключением AK / HI, которые досадно размещены в Мексике для этого шейп-файла:

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
Ари Б. Фридман
источник
Предыдущий ответ по проектированию с использованием пакета proj4 здесь . Хотя не пробовал это с SpatialPolygonsDataFrame.
Симбамангу
На самом деле похоже, что proj4 не работает с пространственными объектами - но смотрите ответ ниже.
Симбамангу
2
Всегда есть представление пространственных задач: cran.r-project.org/web/views/Spatial.html и мои заметки о пространственных данных [бесстыдный плагин]: maths.lancs.ac.uk/~rowlings/Teaching/UseR2012
Spacedman

Ответы:

44

Вы можете использовать spTransform()методы в rgdal - используя ваш пример, вы можете преобразовать объект в NAD83 для Канзаса (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)

проектируется

Чтобы сохранить его в новой проекции:

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

РЕДАКТИРОВАТЬ : Или, согласно предложению @ Spacedman (который записывает файл .prj с информацией CRS):

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

Если вы не уверены, с какого CRS проектировать, обратитесь к следующему сообщению:

И если кто-то хочет определить / назначить CRS, когда у данных его нет, обратитесь к:

Simbamangu
источник
10
обратите внимание, что writePolyShape НЕ записывает файл .prj! Вы должны использовать writeOGR из rgdal (и использовать readOGR для чтения шейп-файлов), если вы хотите записать и прочитать файл .prj, чтобы установить CRS ваших пространственных объектов в R!
Spacedman
Намного лучше (отредактировано соответственно) - спасибо; не понял, что создает файл .prj! Между прочим, потрясающая шпаргалка на вашей странице.
Симбамангу
1
Странно, как проекция в Мексике влияет на внешний вид насекомых Аляски и Гавайев :-).
whuber
@whuber - хм, да ... кто-то отредактировал мою публикацию, в которой не было настоящих карт, показывающих эти довольно неуместные вставки ... никогда не планировал их, чтобы увидеть, что они там.
Симбамангу
@Simbamangu Извините, я забыл, что этот файл .shp неуместно содержит вставки, когда я пытался помочь в добавлении графиков!
Ари Б. Фридман
8

С момента появления sf-пакета (посмотрите виньетки sf1 , sf2 , sf3 , sf4 и руководство по миграции здесь ) вы можете использовать st_transform()для перепроецирования ваших векторных данных:

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 заменит sp в будущем и благодаря своей простоте и скорости имеет несколько преимуществ по сравнению с sp.

andschar
источник