Объединение географических слоев с различными проекциями в R

Question

РЕДАКТИРОВАТЬ: я немного перефразировал заголовочный вопрос и скорректировал текст, чтобы ответить на комментарий с помощью @ DWin.

Объединение географических слоев, которые проецируются и не проецируются, может быть затруднительным.Похоже, что часто требуется некоторая трансформация, поскольку географические слои принадлежат разным продуктам и издателям.

Мне известно, что в R есть несколько инструментов для выполнения географических преобразований.Например:

1. Для объектов класса Spatial* в пакете sp может использоваться функция spTransform() в пакете rgdal;и
2. Для объектов класса Raster* в пакете raster можно использовать функцию projectRaster().

Вот конкретная задача, которую я хотел бы выполнитьвыполнить в R: преобразование в зону 15N сетки UTM (Datum: NAD83) слой многоугольников, описывающий озера в проекции зоны 15N сетки (UatM: NAD27) (в формате шейп-файла ESRI).

Answer 1

Полезная вещь - это база данных epsg, включенная в rgdal.

epsgs = make_EPSG()
subset(epsgs,grepl("15N",epsgs$note))

[etc]
      code
2703 26715                         # NAD27 / UTM zone 15N  [etc]
2851 26915                         # NAD83 / UTM zone 15N  [etc]
[etc]

Эти коды - то, что вам нужно в spTransform.Если ваши озера находятся в шейп-файле с этой проекцией NAD27, то:

require(maptools)
lakes = readShapeSpatial("lakes.shp")
proj4string(lakes)=CRS("+init=epsg:26715")

должно дать вам озера в том виде, в каком они есть (заметьте, я не думаю, что readShapeSpatial будет читать файл .prj с набором шейп-файлов, поэтомумы установили это здесь явно)

Теперь, чтобы преобразовать в базовую версию NT83 зоны UTM 15N:

lakes83 = spTransform(lakes,CRS("+init=epsg:26915"))

Растры немного сложнее, так как они обычно включают деформацию, так что в итоге вы получитерегулярная сетка в вашей системе координат проекции - вы не можете просто преобразовать координаты углов ...