1

我有一个空间多边形对象和一个空间点对象。后者是从数据帧中的 xy latlong 数据向量（分别称为纬度和经度）创建的，而前者只是使用 rgdal 直接读入 R。我的代码如下：

boros <- readOGR(dsn = ".", "nybb")
rats <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv", header = TRUE)
coordinates(rats) <- ~longitude + latitude

此时，两个空间对象都没有投影。如果我按如下方式投影这些对象：

proj4string(boros) <- CRS("+proj=lcc")
proj4string(rats) <- CRS("+proj=lcc")

两个对象现在都被投影，并且都将成功地使用 plot() 函数映射，如下所示：

plot(boros)
plot(rats)

但是，当我尝试将它们绘制在一起时：

plot(boros)
plot(rats, add = TRUE)

我只得到第一个情节，没有老鼠对象叠加在 boros 上。然而，这是一个大问题，我没有收到错误消息，所以我无法确定这两个能够相互交谈的空间对象之间的断开连接在哪里。这两个命令运行顺利，没有错误或警告，但我只剩下一个情节。当我使用 proj4string() 检查每个对象的投影时，我会为每个对象返回相同的投影。

我现在花了很多很多时间在几天内尝试各种方法来创建两个空间对象，它们的 CRS 和投影匹配，以便可以使用 plot() 将它们映射在一起。顺便说一句，我采用的一种方法是在 ArcGIS 中为老鼠对象创建一个 shapefile，它可以很好地创建文件。但是我仍然无法在 R 中一起工作这两个空间对象。我为这两个对象尝试了许多不同的投影，两个对象上的 spTransform，似乎没有任何效果。非常感激任何的帮助。我还包含了一个包含我上面描述的 2 个数据文件的保管箱链接： https ://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384

4

1 回答 1

3

因此，正如一些评论指出的那样，问题在于您的数据和地图具有不同的投影。

看起来您的地图来自纽约市城市规划部。shapefile 绝对不在WGS84 (longlat) 中，但 CRS 不包含在文件中（顺便说一句，这很令人失望……）。然而，有一个元数据文件表明 shapefile 被投影为 EPSG 2263。

为了在 R 中使用它，我们需要一个投影字符串。在 R 中得到这个的惯用方法是：

library(rgdal)
EPSG <- make_EPSG()
NY   <- with(EPSG,EPSG[grepl("New York",note) & code==2263,]$prj4)
NY
# [1] "+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"

现在我们可以获取您的地图并将其重新投影到 WGS84 中，或者获取您的数据并将其重新投影到地图 CRS 中。

setwd("< directory with all your files >")
data   <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv")

# First approach: reproject map into long-lat
wgs.84       <- "+proj=longlat +datum=WGS84"
map          <- readOGR(dsn=".",layer="nybb",p4s=NY)
map.wgs84    <- spTransform(map,CRS(wgs.84))
map.wgs84.df <- fortify(map.wgs84)
library(ggplot2)
ggplot(map.wgs84.df, aes(x=long,y=lat,group=group))+
  geom_path()+
  geom_point(data=data, aes(x=longitude,y=latitude, group=NULL),
             colour="red", alpha=0.2, size=1)+
  ggtitle("Projection: WGS84 longlat")+
  coord_fixed()


# Second approach: reproject data
map.df <- fortify(map)
rats <- SpatialPoints(data[,c("longitude","latitude")],proj4string=CRS(wgs.84))
rats <- spTransform(rats,CRS(NY))
rats.df <- data.frame(coordinates(rats))
ggplot(map.df, aes(x=long,y=lat,group=group))+
  geom_path()+
  geom_point(data=rats.df, aes(x=longitude,y=latitude, group=NULL),
             colour="red", alpha=0.2, size=1)+
  ggtitle("Projection: NAD83.2263")+
  coord_fixed()

中央公园没有老鼠？

于 2014-04-04T00:16:30.800 回答