Невозможно построить два пространственных объекта с одинаковыми CRS в R

Question

У меня есть два пространственных объекта, один - объект пространственного многоугольника, а другой - файл.csv, который я превратил в объект пространственных точек. Первый - это официальный файл формы чилийского правительства для одной из его коммун, другой был создан путем геокодирования с помощью API HERE, уличных адресов той же коммуны.

Сначала я загрузил объект пространственного многоугольника readOGR от:

quilpue <- readOGR( dsn= getwd() , layer="quilpue-rgdal", 
                encoding = "UTF-8") 

Затем я загрузил файл.csv в R и преобразовал его в объект пространственных точек с coordinates()функция от sp пакет.

pointsCoords<- read.csv("../quilpueR/quilpueLayer.csv", header = TRUE)
coordinates(pointsCoords) <- ~Longitude+Latitude

Затем я проверил проекцию каждого объекта.

proj4string(quilpue) 
proj4string(pointsCoords)

"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" а также NAсоответственно.

Единственная проекция, которая работает для pointsCoords было CRS("+init=epsg:3857"), Поэтому я назначил эту проекцию quilpue

proj4string(pointsCoords) <- CRS("+init=epsg:3857")
quilpue_prj <- spTransform(quilpue, CRSobj = CRS(proj4string(pointsCoords)))

Тем не менее, когда я проверяю расширение обоих объектов с extent() от raster() пакет, они не перекрываются.

extent(quilpue_prj)
class       : Extent 
xmin        : -7957703 
xmax        : -7946463 
ymin        : -3907594 
ymax        : -3898059 

extent(pointsCoords)
class       : Extent 
xmin        : -71498550 
xmax        : -71334950 
ymin        : -33133030 
ymax        : -32769810 

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

plot(quilpue_prj)
plot(pointsCoords, add = TRUE) 

Чтобы проверить, была ли проблема с шейп-файлом или файлом.csv, я открыл оба файла на Maptitude другое программное обеспечение ГИС, и ему удается автоматически наложить их. Я хотел бы иметь возможность сделать то же самое в R.

Answer 1

Мне удается решить проблему, но я не понимаю, почему это сработало. После загрузки файла.csv и использования

coordinates(pointsCoords) <- ~Longitude+Latitude 

чтобы создать объект пространственных точек, я использовал projection() функция от rasterПакет для присвоения ему проекции:

projection(pointsCoords) = "+init=epsg:4326"

Затем я преобразовал проекцию объекта пространственного многоугольника quilpueсначала "+init=epsg:3857" а затем "+init=epsg:4326":

quilpue <- spTransform(quilpue, 
             CRSobj = CRS("+init=epsg:3857"))

quilpue <- spTransform(quilpue, 
                       CRSobj = CRS("+init=epsg:4326"))

С bbox() Я проверяю диапазон каждого пространственного объекта:

bbox(pointsCoords)

            min       max
Longitude -71498550 -71334950
Latitude  -33133030 -32769810

bbox(quilpue)

    min       max
x -71.48526 -71.38429
y -33.09254 -33.02075

И обратите внимание, они были очень похожи, и это pointsCoords содержался в quilpue, Единственным предупреждением было то, что quilpuecoords был "." после первых двух цифр, поэтому я использовал gsub добавить "." к coords в pointsCoords,

dfcoords <- as.data.frame(pointsCoords@coords)
dfcoords$Longitude <- as.numeric(gsub("([[:digit:]]{6,6})$", ".\\1", 
                      dfcoords$Longitude)) 
dfcoords$Latitude <- as.numeric(gsub("([[:digit:]]{6,6})$", ".\\1", 
                      dfcoords$Latitude)) 
coordinates(dfcoords) <-  ~Longitude+Latitude

И назначить измененный coordsк оригинальным.

pointsCoords@coords <- dfcoords@coords

Тогда я смог использовать over() и построить пространственные объекты.

df_over <- over(quilpue_prj, pointsCoords)
plot(quilpue)
plot(pointsCoords, add = TRUE)

Answer 2

Я думаю, что вы прошли через некоторые ненужные шаги, чтобы ответить на ваш вопрос.

Чтобы построить их вместе, они просто должны быть в одном CRS:

library(rgdal)    

quilpue <- spTransform(quilpue, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) #or any CRS you wish to use

pointsCoords <- spTransform(pointsCoords, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"")) #or any CRS you wish to use

plot(quilpue)
plot(pointsCoords, add = T)

Также полезно проверить размеры ваших функций, чтобы убедиться, что они совпадают. Иногда случается, что рассматриваемые функции находятся в одном и том же CRS, но из-за обработки или большого количества преобразований они заканчиваются в некоторой степени. Проверьте это с:

extent(quilpue)
extent(pointsCoords)