R 将 CSV 映射到 SHP

Question

我正在尝试在 shapefile 上绘制纬度和经度坐标。我被告知他们有相同的坐标系（NAD27）。研究后我尝试了很多方法，但没有一个是完美的。他们要么映射一个或另一个，但不能同时映射两者。下面是一个读取它们的示例，但我什至不知道如何开始绘图，因为我的 csv 来自 sp 包，而 shp 来自 maptools 包。

pts <- read.csv("locations.csv")
pts <- pts[,4:5]
names(pts) <- c("Latitude", "Longitude")
pts <- pts[complete.cases(pts),]
pts <- SpatialPoints(pts)
border <- maptools::readShapePoly("landman_one shape.shp")

任何帮助将不胜感激。

以下是不同的结构。

> str(pts)
Formal class 'SpatialPoints' [package "sp"] with 3 slots
  ..@ coords     : num [1:372, 1:2] 30.9 30.9 31 31 31 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:372] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "Latitude" "Longitude"
  ..@ bbox       : num [1:2, 1:2] 29.5 -105.2 36 -76.1
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "Latitude" "Longitude"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

> str(border)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1 obs. of  2 variables:
  .. ..$ Id       : int 0
  .. ..$ PERIM_GEO: num 998432
  .. ..- attr(*, "data_types")= chr [1:2] "N" "F"
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 922511 545445
  .. .. .. .. .. .. ..@ area   : num 3.45e+10
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:252, 1:2] 954182 954073 954006 953914 953828 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 922511 545445
  .. .. .. ..@ ID       : chr "0"
  .. .. .. ..@ area     : num 3.45e+10
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 819851 386743 1042230 669865
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

另外，我对预测不太熟悉，因此我们将不胜感激。另外，如果有办法用同一个包做到这一点，我会愿意接受的。

Answer 1

问题在于您的两个空间数据集都没有投影集并且使用不同的坐标系。您可以在@proj4string 行中看到这一点，它们当前都是NA。确定坐标系不同的一种快速方法是使用@bbox 插槽（“边界框”）。在pts 中，这些是经度/纬度，因为它们分别在 -180/+180 和 -90/+90 范围内。另一方面，border 的值采用完全不同的单位 (819851 386743 1042230 669865)。这些看起来都不在NAD27 坐标系中。

要解决这个问题，首先我会使用rgdal::readOGR() 或sf::read_sf() 读取您的border 形状文件。如果存在这些函数，则这些函数会加载到投影中，因此这应该会告诉您该数据在哪个投影中。

其次，要设置pts的投影，可以使用proj4string()：

proj4string(pts) = CRS("+init=epsg:4326")

一旦pts 和border 都有一个投影集，您可以使用spTransform() 将一个坐标系转换到另一个坐标系。因为我不知道 border 的 CRS/投影是什么，所以我将其转换为 WGS84（pts 的投影），但只需将它们交换为将 pts 转换为与 border 相同的投影一次你知道它的 CRS。

border = spTransform(border, CRS("+init=epsg:4326"))

一旦在同一个投影中，它们应该一起绘制。您可以确认他们的预测与：

proj4string(border) == proj4string(pts)