将纬度/经度或 URM 坐标添加到 R 中的 shapefile

Question

我是 R 中的 GIS 新手，我正在尝试将 lat/long 或 UTM 坐标添加到 shapefile。我从这里下载了芝加哥市的边界（City_Boundaries.shp）：http://www.cityofchicago.org/city/en/depts/doit/supp_info/gis_data.html

我加载了 maptools 和 rgeos 库：

library(maptools)
library(rgeos)
library(rgdal)

我将数据带入 R 并尝试为 16T 区添加 UTM 代码：

city<- readShapeSpatial("C:/Users/Luke.Gerdes/Documents/GIS Files/Chicago/City_Boundary.shp")
proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84")

但是，生成的数据对我来说没有意义。当我查看“city”中的“coords”槽时，坐标值是，例如，X=1092925 和 Y=1944820。我使用外部工具 (http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html) 来查找 GPS 坐标。结果是：lat=17.511，long=-81.42。这是牙买加和洪都拉斯海岸之间的某个地方。在我看来，shapefile 坐标存在于它们自己的宇宙中；正如 shapefile 的名称所暗示的那样，坐标提供了对城市形状的准确描述，但这些坐标不会自动映射到地球上。

我的最终目标是确定是否在芝加哥发生了一些地理标记事件。我很乐意将这些按纬度/经度列出的点转换为 UTM，如下所示：

  SP <- SpatialPoints(cbind(-87.63044, 41.79625), proj4string=CRS("+proj=longlat +datum=WGS84"))
  SP <- spTransform(SP, CRS("+proj=utm +north +zone=16T +datum=WGS84"))

如果这样更有意义，我也愿意尝试以原始格式处理事件数据。最终，我需要帮助将芝加哥城市边界变成一个多边形，我可以对照我的事件数据进行检查。像这样的：

  gContains(city, SP)

如何将我的 shapefile 和事件数据放入一个通用（且准确）的参考框架中，以便我可以检查点是否在城市中？非常感谢您提供的任何帮助。

Answer 1

一般建议是使用rgdal::readOGR 读取 shapefile，如

library(rgdal)
city = readOGR(".", "City_Boundary")

这不仅会给city 提供适当的 CRS：

proj4string(city)
[1] "+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

但如果您想在不重新投影的情况下更改它，也会警告您：

proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84")
Warning message:
In `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
  A new CRS was assigned to an object with an existing CRS:
+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
without reprojecting.
For reprojection, use function spTransform in package rgdal

意思是这样做之后，city 是错误。您使用spTransform 重新投影city，如

city = spTransform(city, CRS("+proj=utm +north +zone=16T + datum=WGS84"))

您的rgeos::gContains 调用可能会失败，因为 utm 的两个 CRS 字符串不同；删除前者中+ datum=WGS84中的空格，你会没事的（TRUE）。